kzmig.c

seismic software,very useful

		ltrace = 0;
	}

	/* 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;
    
	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);
*/

	    		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)*nlt*nst*nzt;
   			    	fseek2g(ttfp,nseek*sizeof(float),SEEK_SET);

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

				/* compute residual t */
				if(dlt<999990) 
				  /* not a slant line */ 
      	    			  resit_(tt,&nst,&nlt,&nzt,&izw,&nr,&dst,&dlt,
				    &dr,&fst,&flt,&xs,&ys,tb0);
				else   
				  /* a slant line output */ 
      	    			  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," 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;
			}

			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 aperature \n",jtrace);
			continue;
		}

 		/* update minimum and maximum (s,l) coordinates	*/
		if(ss[ktrace]<ssmin) ssmin = ss[ktrace];
		if(gs[ktrace]<ssmin) ssmin = gs[ktrace];
		if(sl[ktrace]<slmin) slmin = sl[ktrace];
		if(gl[ktrace]<slmin) slmin = gl[ktrace];
		if(ss[ktrace]>ssmax) ssmax = ss[ktrace];
		if(gs[ktrace]>ssmax) ssmax = gs[ktrace];
		if(sl[ktrace]>slmax) slmax = sl[ktrace];
		if(gl[ktrace]>slmax) slmax = gl[ktrace];

		/* find traveltime tables for this gather */
		nxs0 = (ssmin-fxst)/dxst;
		nys0 = (slmin-fyst)/dyst;
		nxs = 1+(int)((ssmax-fxst)/dxst+0.99)-nxs0;
		nys = 1+(int)((slmax-fyst)/dyst+0.99)-nys0;
 
 		if(nys*nxs>ntab){
			range = 0;
			nxs = nxsb;
			nys = nysb;
			nxs0 = nxs0b;
			nys0 = nys0b;
			continue;
		}else{
  			nxsb = nxs;
			nysb = nys;
			nxs0b = nxs0;
			nys0b = nys0;
		}


		iofs[ktrace] = io + 1;

 		/* apply tpow */
                if(itgain==1) for(it=0;it<nt;it++) tra.data[it] *= tgain[it];
		for(it=0;it<nt;it++) tras[it+ktrace*nt] = tra.data[it];

			
		itmp = (tra.mute-tra.delrt);
 		itmp = itmp*1000/(int)tra.dt*ifact + 1; 
		if(itmp<1) itmp = 1;
		if(itmp>ntl) itmp = ntl; 
		imutel[ktrace] = itmp;
 
 
/*
fprintf(stderr,"il=%d is=%d ml=%g ms=%g cdp=%d offset=%d io=%d tracl=%d\n",
		il,is,ml,ms,tra.cdp,tra.offset,io,tra.tracl); 
*/

		if(il>=0 && il<nlout && is>=0 && is<nsout) {
			if(ntrend==0) {
			   if(abs(abs(tra.offset)-ofo[io])
                                <abs(off[io+is*nofo+il*nofo*ns]-ofo[io])) {
				off[io+is*nofo+il*nofo*ns] = abs(tra.offset);
				lpos=io+is*nofo+il*nofo*ns;
				lpos=lpos*240;
				fseek2g(hdrfp,lpos,0);
 				efwrite(&tra,sizeof(char),240,hdrfp);
 			   } 
			} else {
				if(lend==lstart) {
					ps = ms;
					pl = lstart;
				} else {
					tmp = (send-sstart)/(lend-lstart);
					pl = ml + tmp*(lstart*tmp+ms-sstart); 
					pl = pl/(1.+tmp*tmp);
					ps = tmp*(pl-lstart) + sstart;
				}
				if(dsout!=0.) {
					tmp = (ps - sstart)/dsout + 0.5;
					isend = tmp; 
				} else { 
					tmp = (pl - lstart)/dlout + 0.5;
					isend = tmp;
				}
				if (isend>=0 && isend<ntrend) {
				      tmp = sqrt( (ms-s[isend])*(ms-s[isend])
				          + (ml-l[isend])*(ml-l[isend]) );
				      if (tmp<disend[io+isend*nofo]) { 
					 disend[io+isend*nofo] = tmp;
					 tra.ungpow  = tmp;
				         lpos=io+isend*nofo;
				         lpos=lpos*240;
					 fseek2g(hdrfp,lpos,0);
 				         efwrite(&tra,sizeof(char),240,hdrfp);
				      }
			        }  
			}
		}

 		ktrace += 1;

		if(jtrace/1000>m1000) { 
		fprintf(jpfp," processing done for %d input traces\n",jtrace);
			m1000 = jtrace/1000;
		fflush(jpfp);
		}
	    }while(done);

/*
fprintf(stderr,"jtrace=%i,range=%i ntrace=%d \n",jtrace,range,ntrace);
*/
	
 		if(range){
			if(ntrace==ntras) gottrace = 0;
			else if(!fgettr(infp,&tra))   
				gottrace = 0;
		} else {
			jtrace -= 1;
			ntrace -= 1;
		}

		if(gottrace) {
			if(gath==1){
	 			newheadx = tra.sx;
	 			newheady = tra.sy;
  			} else if(gath==2){
	 			newheadx = tra.gx;
	 			newheady = tra.gy;
 			}
			if(ktrace==mtrace || !range) ready = 1;
			if(newheadx!=oldheadx || newheady!=oldheady){
				oldheadx = newheadx;
				oldheady = newheady;
				ready = 1;
 			}
			if(ready){  
				range = 1;
   				ssmin = slmin = 99999999.0;
				ssmax = slmax = -99999999.0;
			}
			if(!ktrace) ready = 0;
		}
	
		if(ndskwt>=ntrdsk && ntrdsk>0) {
                /* update disk file */