kzmig.c

seismic software,very useful

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

		/* trace select */
        if(ikey==1) {
			gethval(&tra,indxtrk,&trkval);
			inttrace = vtoi(trktype,trkval);
			if( (inttrace-strace)%jthtrace!=0 
				|| inttrace<mintrace || inttrace>maxtrace) continue;
			gethval(&tra,indxlnk,&lnkval);
			intline = vtoi(lnktype,lnkval);
			if( (intline-sline)%jthline!=0 
				|| intline<minline || intline>maxline ) 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-1) {
			tmp = tra.offset;
			if(iofoin==0 || nofo==1 ) {
				if(tmp<0) tmp = -tmp;
			}
			tmp = tmp - ofo[io];
			if(tmp<0.) tmp =  - tmp;
			if(tmp>obtol*1.001) io=-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);
			*/
			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(ikey==1) {
				ms = s[il*nsout+is];
				ml = l[il*nsout+is];
				fline = (ml - l1)*(ln3-ln1)/(l3-l1) + ln1 + 0.5;
				ftrace = (ms - s1)*(tr2-tr1)/(s2-s1) + tr1 + 0.5;
				intline = fline;
				inttrace = ftrace;
				itov(trktype,&trkval,inttrace);
				itov(lnktype,&lnkval,intline);
				puthval(&tra,indxtrk,&trkval);
				puthval(&tra,indxlnk,&lnkval);
			}
			tmp = fzout[is+il*nsout]+0.5;
			itmp = tmp;
			tra.delrt = (short)itmp;
			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);
 				fwrite(&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);
 				     	if(fwrite(&tra,sizeof(char),240,hdrfp)!=240)
							err(" error write %s \n",diskhdr);
				      }
			        }  
			}
		}

 		ktrace += 1;

		if(jtrace/1000>m1000) { 
		/*
		gettime(datetime);
		fprintf(jpfp," processing done for %d input traces at %s\n",
			jtrace,datetime);
		*/
		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 */
                fldfp = efopen(diskfld,"r+");
                fseek2g(fldfp,0,1);
                fwrite(fold[0],sizeof(float),
				nz*nlout*nsout*nofo,fldfp);
                efclose(fldfp);
                imgfp = efopen(diskimg,"r+");
                fseek2g(imgfp,0,1);
                fwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
                efclose(imgfp);
                fprintf(hffp,"Number of Traces Processed: %d %d \n",