kzmig.c.9.2.99

seismic software,very useful

					puthval(&tr,indxtrk,&trkval);
					puthval(&tr,indxlnk,&lnkval);
				}
				if(dzo<=32.) { 
					tr.dt = dzo*1000+0.5;
				} else {
					tr.dt = dzo;
				}
				tr.dt = dzo*1000.+0.5;
				for(io=0;io<nofo;io++) {
					tr.cdpt  = io+1;
					ofs = ofo[io];
					tr.sx = mx-ofs/2;
					tr.gx = mx+ofs/2;
					tr.sy = my;
					tr.gy = my;
					fwrite(&tr,sizeof(char),240,hdrfp);
				}
			}
		}
		fprintf(jpfp," initializing diskhdr file done \n");
		for(ix=0;ix<nofo*nsout*nlout;ix++) off[ix] = 99999999;
		fflush(hdrfp);
		
	}


	ktrace = 0;
	m1000 = 0;
	oldheadx = newheadx = 0;
	oldheady = newheady = 0;
	ssmin = slmin = 99999999.0;
	ssmax = slmax = -99999999.0;
 
	if(itgain==1) {
		for(it=0;it<nt;it++) {
			tmp = t0 + it*dt;
			if(tmp!=0.0) {
				tgain[it] = pow(tmp,tpow);
			} else {
				tgain[it] = 0.;
			}
		}
 	} else {
		for(it=0;it<nt;it++) tgain[it] = 1.0;
	}
		
	fprintf(jpfp," \n");
	fprintf(jpfp," Migration Parameters \n");
	fprintf(jpfp," -------------------- \n");
	fprintf(jpfp," x1=%g y1=%g s1=%g l1=%g cdp1=%d\n",x1,y1,s1,l1,cdp1);
	fprintf(jpfp," x2=%g y2=%g s2=%g l2=%g cdp2=%d\n",x2,y2,s2,l2,cdp2);
	fprintf(jpfp," x3=%g y3=%g s3=%g l3=%g cdp3=%d\n",x3,y3,s3,l3,cdp3);
	fprintf(jpfp," sstart=%g lstart=%g ns=%d nl=%d ds=%g dl=%g \n",
		sstart,lstart,nsread,nlread,ds,dl);
	if(ntrend>0) fprintf(jpfp," send=%g lend=%g ntrend=%d dend=%g \n",
			send,lend,ntrend,dend);
	if(ntrend>0) fprintf(jpfp," dismax=%g aa=%g bb=%g cc=%g \n",
			dismax,aa,bb,cc);
	fprintf(jpfp," fzo=%g nzo=%d dzo=%g\n",fzo,nzo,dzo);
	fprintf(jpfp," nds=%d ndl=%d dds=%g ddl=%g \n",
		nds,ndl,dds,ddl);
	fprintf(jpfp," nss=%d nll=%d is0=%d il0=%d \n",
		nss,nll,is0,il0);
	fprintf(jpfp," cdppdds=%d cdppds=%d cdppddl=%d cdppdl=%d cdpnum=%d\n",
		cdppdds, cdppds, cdppddl, cdppdl, cdpnum);
	fprintf(jpfp," ntrpre=%d nt=%d ntl=%d dt=%g t0=%g nfft=%d\n",
		ntrpre,nt,ntl,dt,t0,nfft);
 	fprintf(jpfp," tpow=%g zpow=%g ofomin=%g dofo=%g nofo=%d isave=%d \n",
		tpow,zpow,ofomin,dofo,nofo,isave);
	fprintf(jpfp," apers=%g aperl=%g angmax=%g tapecntl=%d\n",
		apers,aperl,angmax,tapecntl);
	fprintf(jpfp," nlcore=%d mtrace=%d ntras=%d mlimit=%d ntrdsk=%d \n",
		nlcore,mtrace,ntras,mlimit/(1024*1024),ntrdsk);
	fprintf(jpfp," f0=%g df=%g nf=%d ftaper=%g ksmax=%g klmax=%g\n",
		f0,df,nf,ftaper,ksmax,klmax);
	fprintf(jpfp," offmax=%g ntab=%d nr=%d \n",offmax,ntab,nr);
 	fprintf(jpfp," fzt=%g dzt=%g nzt=%d \n",fzt,dzt,nzt);
	fprintf(jpfp," fst=%g dst=%g nst=%d \n",fst,dst,nst);
	fprintf(jpfp," flt=%g dlt=%g nlt=%d \n",flt,dlt,nlt);
	fprintf(jpfp," fxst=%g dxst=%g nxst=%d \n",fxst,dxst,nxst);
	fprintf(jpfp," fyst=%g dyst=%g nyst=%d \n",fyst,dyst,nyst);
	fprintf(jpfp," v0=%g dvz=%g gath=%d \n",v0,dvz,gath);
	fprintf(jpfp," diskimg=%s\n",diskimg);
	fprintf(jpfp," diskfld=%s\n",diskfld);
	fprintf(jpfp," diskhdr=%s\n",diskhdr);
	fprintf(jpfp," backupi=%s\n",backupi);
	fprintf(jpfp," backupo=%s\n",backupo);
	fprintf(jpfp," hisfile=%s\n",hisfile);
	fprintf(jpfp," ncpu=%d traceout=%d\n",ncpu,traceout);
	if(on2trace==1) fprintf(jpfp," total_number_of_input_trace=%d \n",ntotal);
	fprintf(jpfp," \n");


/*
	gettime(datetime);
 	fprintf(jpfp," start migration at %s ... \n",datetime);
*/
 	fprintf(jpfp," start migration at ... \n");
 	fflush(jpfp); 

 	m1000 = ntrpre/1000;
	ntrace = 0;

	if(tapecntl==1) {
		if(jtrace>0) {
        	lpos = nt*sizeof(float)+240;
            lpos = lpos * jtrace + 3600;
            fseek2g(infp,lpos,0);
            if(outputonly==0) fgettr(infp,&tra);
        }
		ltrace = jtrace;
		jtrace = ntrpre;
	} else if(tapecntl==-1) {
		ltrace = 0;
		jtrace = ntrpre;
		lpos = nt*sizeof(float)+240;
		lpos = lpos * ntrpre + 3600;
		fseek2g(infp,lpos,0);
		if(outputonly==0) fgettr(infp,&tra);
	} else if(tapecntl==0) {
		ltrace = jtrace;
		jtrace = ntrpre - jtrace;
	} else if(tapecntl==2) {
		ltrace = jtrace;
		jtrace = ntrpre;
	}

	/* set old header for the first trace */
	if(gath==1){
	 	oldheadx = tra.sx;
 	 	oldheady = tra.sy;
 	} else if(gath==2){
	 	oldheadx = tra.gx;
 	 	oldheady = tra.gy;
 	} 
	range = gottrace = 1;
	ready = done = 0;
	nxsb = nysb = nxs0b = nys0b = 0;

	if (outputonly==1) goto no_mig;
    
	do {
 		
		/* if gather has changed or no more input traces */
		if ( ready || (!gottrace && ktrace>0) ) {
  			fxs = fxst+nxs0*dxst;
			fys = fyst+nys0*dyst;

/*
	fprintf(stderr,"fxs=%g nxs=%d fys=%g nys=%d \n",fxs,nxs,fys,nys);
				fprintf(jpfp,"read ttable time=%f \n",walltime());
*/

	    	for(iys=0,ys=fys; iys<nys; ++iys,ys+=dyst) 
  		    for(ixs=0,xs=fxs; ixs<nxs; ++ixs,xs+=dxst){
 				itab = iys*nxs+ixs;
			
 		    	nseek = nxst*(nys0+iys)+nxs0+ixs;
				nseek = nseek*nlt*nst*nzt;
   		    	fseek2g(ttfp,nseek*sizeof(float),SEEK_SET);

 				fread(tt,sizeof(float),nst*nlt*nzt,ttfp);

				/* compute residual t */
				if(dlt<999990) 
      	   			resit_(tt,&nst,&nlt,&nzt,&izw,&nr,&dst,&dlt,
				    	&dr,&fst,&flt,&xs,&ys,tb0);
				else   
      	    		resit0_(tt,&nst,&nzt,&izw,&nr,&dr,
			    		st,lt,&xs,&ys,tb0);
  
				/* intrepolation t to output lateral positions */ 
   				latint_(&nsout,&nlout,&nz,s,l,ttab[itab],
			 		&nst,&nlt,&nzt,&nz0,&fst,&flt,&dst,&dlt,tt);
	      	}
 
			/* migration of ktrace */
			migs(nlout,nsout,mig,nz,
				ktrace,ntl,dldt,trace, 
				tras,dtl,tminl,nt, 
				ss,sl,gs,gl, 
				imutel,s,l,
				iofs,apers,aperl,fold,
				work1,work2,wsave,tracef,sqrtf,
				tahd,ifcut,ltaper,ksmax,klmax,
				f0,df,nf,ftaper,nsave,nfft, 
				dzt,fz,angmax, 
				nxs,nys,fxs,fys,dxst,dyst,
				nzo,fzo,dzo,
				nr,dr,tb,pb,ttab,gath);

			itrace = itrace + ktrace;
			ndskwt = ndskwt + ktrace;
 			ktrace = 0;
/*	fprintf(jpfp," processing done for total %d input traces \n",jtrace);*/
 	fflush(jpfp); 

                        /* if no more traces, break */
                        if (!gottrace) break;

			ready = 0;
		}

		if(!gottrace) break;

	    do {
		jtrace = jtrace + 1;
		if(jtrace<=ntrpre) continue;
		ntrace = ntrace + 1;

		/* dead trace */ 
		if(tra.trid!=1) {
		fprintf(btfp," trid not 1, at trace=%d cdp=%d cdpt=%d\n",
			jtrace,tra.cdp,tra.cdpt); 
			fflush(btfp);
			continue;
		}

		/* compute rms and max of trace amplitudes */
		itmp = (tra.mute-tra.delrt);
		itmp = itmp*1000/(int)tra.dt; 

		ampmax = 0.;
		amprms = 0.;
		if(itmp<0) itmp = 0;
		it0 = 0;

		for(it=itmp;it<nt;it++) {
			tmp = fabs(tra.data[it]);
			if(it0==0 && tmp>0.) it0 = it;
			amprms = amprms + tmp*tmp;
			if(tmp>ampmax) ampmax = tmp;
 		}
		itmp = it0;
		
		it0 = tra.delrt + itmp * (int)tra.dt/1000; 
		if(it0>tra.mute) tra.mute = it0; 
		
		if(nt>itmp) amprms = sqrt(amprms/(nt-itmp));
		if(amprms==0.) {
		fprintf(btfp," zero trace, at trace=%d cdp=%d cdpt=%d\n",
			jtrace,tra.cdp,tra.cdpt); 
			fflush(btfp);
			continue;
		} else {
			if (ampmax>amprms*rmsmar) {
fprintf(btfp," bad trace with ampmax=%g  rms=%g, at trace=%d cdp=%d cdpt=%d\n",
		ampmax,amprms,jtrace,tra.cdp,tra.cdpt); 
				fflush(btfp);
				continue;
			}
		}


		/* source and group x, y locations */
		xs = tra.sx;
		ys = tra.sy;
		xg = tra.gx;
		yg = tra.gy;

		if(tra.scalco>1) {
			xs = xs *tra.scalco;
			ys = ys *tra.scalco;
			xg = xg *tra.scalco;
			yg = yg *tra.scalco;
		} else if(tra.scalco<0) {
			xs = xs / (-tra.scalco);
			ys = ys / (-tra.scalco);
			xg = xg / (-tra.scalco);
			yg = yg / (-tra.scalco);
		}

 		/* source and group s, l locations */
		xy2sl(x1,y1,s1,l1,x2,y2,s2,l2,x3,y3,s3,l3,xs,ys,
			&ss[ktrace],&sl[ktrace]);

		xy2sl(x1,y1,s1,l1,x2,y2,s2,l2,x3,y3,s3,l3,xg,yg,
			&gs[ktrace],&gl[ktrace]);
 
 		if( ss[ktrace]<fxst || gs[ktrace]<fxst ||
		    ss[ktrace]>exst || gs[ktrace]>exst ||
		    sl[ktrace]<fyst || gl[ktrace]<fyst ||
		    sl[ktrace]>eyst || gl[ktrace]>eyst ) {

/*
	    	fprintf(jpfp," ss=%g gs=%g fxst=%g exst=%g \n",
				ss[ktrace],gs[ktrace],fxst,exst);
	    	fprintf(jpfp," sl=%g gl=%g fyst=%g eyst=%g \n",
				sl[ktrace],gl[ktrace],fyst,eyst);
	    	fprintf(jpfp," trace %d is out of time tables\n",jtrace);
*/

  			continue;
		}
 		if( (ss[ktrace]-gs[ktrace])*(ss[ktrace]-gs[ktrace])+ 
 		    (sl[ktrace]-gl[ktrace])*(sl[ktrace]-gl[ktrace])>
 			offmax*offmax) {
			/*
 		    fprintf(btfp," trace %d has offset too large \n",jtrace);
			*/
			continue;
		}
		ms = (ss[ktrace] + gs[ktrace])/2.;
		ml = (sl[ktrace] + gl[ktrace])/2.;

		if(ntrend==0) { 

			/*
			itmp = tra.cdp - cdp1;
			if(cdpnum==0) {
				itmp = itmp/ncdppl;
				iiy = itmp;
				iix = tra.cdp - cdp1 - iiy*ncdppl;
			} else {
				itmp = itmp/nlines;
				iix = itmp;
				iiy = tra.cdp - cdp1 - iix*nlines;
			}
			*/

			/* replaced by the following */
			tmp = (ml - l1)/ddl + 0.5;
			iiy = tmp;
			tmp = (ms - s1)/dds + 0.5;
			iix = tmp;

			is = (iix-is0)/nss;
			ids = iix - is*nss - is0;

			if(ids>=0 && ids<nds) {
				is = is * nds + ids; 
			} else {
				is = -1;
			}
		
			il = (iiy-il0)/nll;
			idl = iiy - il*nll - il0;

			if(idl>=0 && idl<ndl) {
				il = il * ndl + idl;
			} else {
				il = -1;
			} 

		} else {
			is = 0;
			il = 0;
		}
		
 			
		tmp = tra.offset;
		if(iofoin==0 || nofo==1 ) {
			if(tmp<0.) tmp = - tmp;
			tmp = (tmp-ofomin)/dofo+.5;
			io = tmp;
			if(tmp<0.) io = -1;
		} else {
			if(tmp<ofol) {
				io = -1;
			} else if(tmp>ofor) {
				io = nofo;
			} else {
				bisear_(&nofo,&one,ofo,&tmp,&io);
				io = io - 1;
				if(io<nofo-1) {
					if(abs(tmp-ofo[io])>abs(tmp-ofo[io+1]))
						io = io + 1;
				}
			}
		}

		if(offnear==1 && io<0 ) io = 0;
		if(offfar==1 && io>=nofo) io = nofo-1;
 
		if(io<0 || io>=nofo) {
/*
 		    fprintf(btfp," trace %d has offset io=%d \n",jtrace,io);
*/
			continue;
		}

		if( ms < sinmin || ms >sinmax) continue;
		if( ml < linmin || ml >linmax) continue;
/*
	fprintf(stderr,"ms=%g ml=%g jtrace=%d dis=%g\n",
		ms,ml,jtrace,disp2l(aa,bb,cc,ms,ml));
*/
		if( ntrend>0 && disp2l(aa,bb,cc,ms,ml)>dismax ) {
			/*
 		    fprintf(btfp," trace %d outside aperture \n",jtrace);