ktmig.c~

seismic software,very useful

	for(it=0;it<ntau;it++) {
		tmig[it] = tau0 + it*dtau;
		tm[it] = tmig[it]*0.5/dtl;
	}

	for(it=0;it<ntau;it++) tau[it] = tau0 + (it+1)*dtau;

	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.;
			}
		}
		for(it=0;it<ntau;it++) {
			tmp = tau0 + it*dtau;
			if(tmp!=0.0) {
				tugain[it] = pow(tmp,tpowaf);
			} else {
				tugain[it] = 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);
	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 ntau=%d ntl=%d nfft=%d\n",
		ntrpre,nt,ntau,ntl,nfft);
	fprintf(jpfp," dt=%f dtau=%f t0=%f tau0=%f \n",
		dt*1000.,dtau*1000.,t0*1000.,tau0*1000.);
	fprintf(jpfp," tpow=%g ofomin=%g dofo=%g nofo=%d isave=%d \n",
		tpow,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 \n",
		nlcore,mtrace,ntras,mlimit/(1024*1024));
	fprintf(jpfp," f0=%g df=%g nf=%d ftaper=%g ksmax=%g klmax=%g\n",
		f0,df,nf,ftaper,ksmax,klmax);
	fprintf(jpfp," t0v=%g dtv=%g ntv=%d \n",t0v,dtv,ntv);
	fprintf(jpfp," s0v=%g dsv=%g nsv=%d \n",s0v,dsv,nsv);
	fprintf(jpfp," l0v=%g dlv=%g nlv=%d \n",l0v,dlv,nlv);
	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," ntrdsk=%d tpowaf=%g \n",ntrdsk,tpowaf);
	fprintf(jpfp," ncpu=%d traceout=%d \n",ncpu,traceout);
	if(on2trace==1) fprintf(jpfp," total_number_of_input_trace=%d \n",ntotal);
	if(ikey==1) { 
		fprintf(jpfp," tr1=%g tr2=%g tr3=%g \n",tr1,tr2,tr3);
		fprintf(jpfp," ln1=%g ln2=%g ln3=%g \n",ln1,ln2,ln3);
		fprintf(jpfp," tracekey=%s linekey=%s \n",tracekey,linekey); 
		fprintf(jpfp," strace=%d jthtrace=%d sline=%d jthline=%d \n",
			strace,jthtrace,sline,jthline);
		fprintf(jpfp," mintrace=%d maxtrace=%d \n",mintrace,maxtrace);
		fprintf(jpfp," minline=%d maxline=%d \n",minline,maxline);
	}
						    


	fprintf(jpfp," \n");
	

	fprintf(jpfp," start migration ... \n");
	fflush(jpfp);

	t0v = t0v * 0.001;
	dtv = dtv * 0.001;

	m1000 = ntrpre/1000;
	ntrace = 0;

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

	i3d = 3;
	if(nlines<10) i3d=2;

	/* initialization of fft and taper for antialiasing filter */
	fhigint_(&f0,&df,&nf,&ftaper,&dtl,ifcut,ltaper,wsave,&nsave,&nfft);

	/* read in velocity file if ivdisk=0 */
	if(ivdisk==0) {
		fseek2g(velfp,0,0);
		efread(velos,sizeof(float),ntv*nsv*nlv,velfp);
	}

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

		/* 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; 
		/* set the mute time to be 10 samples shallower */
		/* just to be save */
		tra.mute = tra.mute - 10*tra.dt/1000;
		if(tra.mute<0) tra.mute = 0;

		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]);


		ms = (ss[ktrace] + gs[ktrace])/2.;
		ml = (sl[ktrace] + gl[ktrace])/2.;

	fprintf(stderr,"ml=%g ms=%g cdp=%d offset=%d \n",
		ml,ms,tra.cdp,tra.offset); 

		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;
		} 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(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;
	 	}

	fprintf(stderr,"il=%d is=%d ml=%g ms=%g cdp=%d offset=%d io=%d\n",
		il,is,ml,ms,tra.cdp,tra.offset,io); 

		if(io<0 || io>=nofo) continue;

		if( ms < sinmin || ms >sinmax) continue;
		if( ml < linmin || ml >linmax) continue;

		iofs[ktrace] = io + 1;

		/* apply rho filter */
		f2p5n_(tra.data,&nt,&i3d,&mm,&ll);
		
		/* 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;



		if(il>=0 && il<nlout && is>=0 && is<nsout) {
			if(ntrend==0) {
			   if(abs(tra.offset)<off[io+is*nofo+il*nofo*nsout]) {
				off[io+is*nofo+il*nofo*nsout] = abs(tra.offset);
				lpos=io+is*nofo+il*nofo*nsout;
				lpos=lpos*240;
				fseek2g(hdrfp,lpos,0);
				tmp = tra.mute*0.001;
				tmp = tmp*tmp-tra.offset*tra.offset/gmin/gmin; 
				if(tmp>=0.) {
					tmp = sqrt(tmp);
					tmp = tmp * 1000.;
					itmp = tra.mute;
					tra.mute = tmp;	
				}
				efwrite(&tra,sizeof(char),240,hdrfp);
				if(tmp>=0.) tra.mute = itmp;
			   } 
			} 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) {
			   	   if (abs(tra.offset)<off[io+isend*nofo]) {
				      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.dz = tmp;
				         lpos=io+isend*nofo;
				         lpos=lpos*240;
					 fseek2g(hdrfp,lpos,0);
				         tmp = tra.mute*0.001;
				         tmp = tmp*tmp-tra.offset*tra.offset/
						gmin/gmin;