ktmig.c

seismic software,very useful

*/

		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; 
				         if (tmp>=0.) {
					    tmp = sqrt(tmp);
					    tmp = tmp * 1000.;
					    itmp = tra.mute;
					    tra.mute = tmp;	
				         }
				         efwrite(&tra,sizeof(char),240,hdrfp);
				      }
				   }
			        }  
			}
		}


		/* find velocity */
		cdpnow = tra.cdp;
		if(cdpnow!=cdppre || vread==0 ) { 
			cdppre = cdpnow;
			tmp = (ms-s0v)/dsv + .5;
			isv = tmp;
			if(isv<0) isv = 0;
			if(isv>nsv-1) isv = nsv - 1;
			tmp = (ml-l0v)/dlv + .5;
			ilv = tmp;
			if(ilv<0) ilv = 0;
			if(ilv>nlv-1) ilv = nlv - 1;

			lpos = (isv + ilv*nsv)*ntv;
			lpos = lpos * sizeof(float);

			if(ivdisk==1) {
				fseek2g(velfp,lpos,0);
				efread(vel,sizeof(float),ntv,velfp);
				if(etatype!=-1) {
					fseek2g(etafp,lpos,0); 
					efread(eta,sizeof(float),ntv,etafp);
				}
			} else {
				for(it=0;it<ntv;it++)
					vel[it] = velos[it+(isv+ilv*nsv)*ntv];
				if(etatype!=-1) {
					for(it=0;it<ntv;it++)
						eta[it] = etas[it+(isv+ilv*nsv)*ntv];
				}
			}
			vread = 1;
		}


		/* compute fovt2 = 4./(v*t)**2    */ 
		for(it=0;it<ntau;it++) {
			tmp = (tmig[it]-t0v)/dtv;
			itv = tmp;
			if(itv<0) {
				v = vel[0];
			} else if(itv>=ntv-1) {
				v = vel[ntv-1];
			} else {
				tmp = tmp - itv;
				v = (1.-tmp)*vel[itv]+tmp*vel[itv+1];
			}
			vr2[ktrace*ntau+it] = v*v;
			tmp = v*tmig[it];
			if(tmp==0.) tmp = v * dtau; 
			fovt2[ktrace*ntau+it] = 4./(tmp*tmp);
		}
		/* compute eta at output travel time */
		if(etatype==0) {
			for(it=0;it<ntau;it++) {
				tmp = (tmig[it]-t0v)/dtv;
				itv = tmp;
				if(itv<0) {
					v = eta[0];
				} else if(itv>=ntv-1) {
					v = eta[ntv-1];
				} else {
					tmp = tmp - itv;
					v = (1.-tmp)*eta[itv]+tmp*eta[itv+1];
				}
				etamig[ktrace*ntau+it] = v;
			}
		} else if(etatype==1) {
			tmp = 0.;
			etaeff[0] = eta[0];
			for(it=1;it<ntv;it++) {
				vr4 = vel[it]*vel[it]*vel[it]*vel[it];
				tmp = tmp + vr4*(1.+8.*eta[it])*dtv;
				etaeff[it] = 0.125*(tmp/(it*dtv+t0v)/vr4-1.);
			}
			for(it=0;it<ntau;it++) {
				tmp = (tmig[it]-t0v)/dtv;
				itv = tmp;
				if(itv<0) {
					v = etaeff[0];
				} else if(itv>=ntv-1) {
					v = etaeff[ntv-1];
				} else {
					tmp = tmp - itv;
					v = (1.-tmp)*etaeff[itv]+tmp*etaeff[itv+1];
				}
				etamig[ktrace*ntau+it] = v;
			}
		}
		ktrace += 1;
			
		if(ktrace==mtrace || ktrace==ntras) {
			migs(nlout,nsout,nofo,mig,ntau,
				imgfp,ktrace,ntl,dldt,trace, 
				tras,dtl,tminl,nt,
				ss,sl,gs,gl, 
				imutel,fovt2,s,l,
				tm,iofs,apers,aperl,fold,
				w1,work1,work2,wsave,tracef,
				ifcut,ltaper,ksmax,klmax,
				f0,df,nf,ftaper,nsave,nfft,
				dtau,tau0,angmax,
				indxw,nindxw,incdxw,
				resamp,ires,ntres,
				ss2,scs,scg,ncpu,
				vr2,vi2,vq4,tau,etamig,etatype);
			itrace = itrace + ktrace;
			ndskwt = ndskwt + ktrace;
			ktrace = 0;
		}


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

		if(ndskwt>=ntrdsk && ntrdsk>0) {
		/* update disk file */
			fldfp = efopen(diskfld,"r+");
			fseek2g(fldfp,0,1);
			efwrite(fold,sizeof(float),nlout*nsout*nofo,fldfp);
			efclose(fldfp);
			imgfp = efopen(diskimg,"r+");
			fseek2g(imgfp,0,1);
			fwrite(mig,sizeof(float),nofo*nsout*nlout*ntau,imgfp);
			efclose(imgfp);
		fprintf(hffp,"Number of Traces Processed: %d %d \n",
				jtrace,jtrace-ntrpre+ltrace);
			fflush(hffp);
			ndskwt=0;
		}
	} while (ntrace<ntras && fgettr(infp,&tra));

	if(ktrace>0) {
		migs(nlout,nsout,nofo,mig,ntau,
			imgfp,ktrace,ntl,dldt,trace, 
			tras,dtl,tminl,nt, 
			ss,sl,gs,gl, 
			imutel,fovt2,s,l,
			tm,iofs,apers,aperl,fold,
			w1,work1,work2,wsave,tracef,
			ifcut,ltaper,ksmax,klmax,
			f0,df,nf,ftaper,nsave,nfft, 
			dtau,tau0,angmax, 
			indxw,nindxw,incdxw,
			resamp,ires,ntres,
			ss2,scs,scg,ncpu,
			vr2,vi2,vq4,tau,etamig,etatype);
		itrace = itrace + ktrace;
		ktrace = 0;
	fprintf(jpfp," processing done for total %d input traces \n",jtrace);
	}

	fprintf(jpfp," %d input traces processed in this run \n",ntrace);
	fflush(jpfp);
	ntrpre = ntrpre + ntrace;


	efclose(hdrfp);
	hdrfp = efopen(diskhdr,"r");
	fseek2g(hdrfp,0,1);

	if(isave==1 && ntrace>0) {
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		efwrite(fold,sizeof(float),nlout*nsout*nofo,fldfp);
		efclose(fldfp);
	} 
	if(isave==0) unlink(diskfld);

	if(ftell(btfp)==0) { 
		efclose(btfp);
		unlink("BAD_TRACE_FILE");
	} else {
		efclose(btfp);
	}
	if(nlcore<nlout || (ipre==1&&ntrdsk==0)) efclose(imgfp);

	if(isave==1 && nlcore==nlout)  {
		imgfp = efopen(diskimg,"r+");
		fseek2g(imgfp,0,1);
		fwrite(mig,sizeof(float),nofo*nsout*nlout*ntau,imgfp);
		efclose(imgfp);
	}

	if(isave==1) fprintf(jpfp," backup to disk done \n");

	fprintf(hffp,"Number of Traces Processed: %d 0 \n",ntrpre);
	efclose(hffp);
	if(ihis==0) unlink(hisfile);


	/* output */
	if(nlcore<nlout) fseek2g(imgfp,0,0);
	fseek2g(hdrfp,0,0);

	if( (ntrace == ntras || feof(infp)!=0) && traceout==1) {
	     fprintf(jpfp," start output ... \n");

	     for(iy=0;iy<nlout;iy++) {
		for(ix=0;ix<nsout;ix++) {
			for(io=0;io<nofo;io++) {
				itmp = (io*nsout*nlout+iy*nsout+ix);
				itmp2 = itmp*ntau;
				bcopy(mig+itmp2,tra.data,ntau*sizeof(float));
				if(fold[itmp]>1.) {
					scale = 1./fold[itmp];
					for(it=0;it<ntau;it++) 
						tra.data[it]=tra.data[it]*scale;
				}

				/* update headers */
				bzero(&tra,240);
				efread(&tra,sizeof(char),240,hdrfp);
				tra.offset = ofo[io];
				tra.tracl = ix + 1;
				tra.tracr = iy + 1;
				tra.cdpt = io + 1;
				tra.dz = 0.;
				ofs = ofo[io];
				if(tra.trid==1) {
					mx = (tra.gx+tra.sx)/2;
                        		my = (tra.gy+tra.sy)/2;
					tmp = tra.gx-tra.sx;
					dd = tra.gy-tra.sy;
					tmp = tmp*tmp + dd*dd;
                   	dd = sqrt(tmp);
                    if(dd>0.) {
                   		tra.sx = mx+ofs*(tra.sx-mx)/dd;
                   		tra.gx = mx+ofs*(tra.gx-mx)/dd;
                   		tra.sy = my+ofs*(tra.sy-my)/dd;
                   		tra.gy = my+ofs*(tra.gy-my)/dd;
                   	} else {
                        tra.sx = mx-ofs/2;
                   		tra.gx = mx+ofs/2;
                       	tra.sy = my;
                       	tra.gy = my;
                    }
				} else {
					ms = s[iy*nsout+ix];
					ml = l[iy*nsout+ix];
					sl2xy(s1,l1,x1,y1,s2,l2,x2,y2,s3,l3,x3,y3,ms,ml,&mx,&my);
					tra.sx = mx-ofs/2;
					tra.gx = mx+ofs/2;
					tra.sy = my;
					tra.gy = my;
				}
 
				tra.delrt = tau0 * 1000.;
				tra.ns = ntau;
				tra.dt = dtau * 1000000.;
				if(tra.trid!=1) {
					tmp = (tau0 + (ntau-1)*dtau)*1000;
					tra.mute = tmp;
				}

				if(itgain==1) {
					for(it=0;it<ntau;it++) 
						tra.data[it]=tra.data[it]
							*tugain[it];

				}

				if(ntrend>0) {
					tmp = (l[iy*nsout+ix] - l1)/ddl + 0.5;
					iiy = tmp;
					tmp = (s[iy*nsout+ix] - s1)/dds + 0.5;
					iix = tmp;
					if(cdpnum==0) {
						tra.cdp = iiy*ncdppl+iix+cdp1;
					} else {
						tra.cdp = iix*nlines+iiy+cdp1;
					}
				}
				if(ikey==1) {
					ms = s[iy*nsout+ix];
					ml = l[iy*nsout+ix];
					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);
				}

				/* output */
				fputtr(outfp,&tra);
			}
		}
	    }
	    fprintf(jpfp," output done \n");
	    fflush(jpfp);
	} 
	
	efclose(outfp);
	

	efclose(hdrfp);
	if(isave==0) unlink(diskhdr);
	if(isave==0) unlink(diskimg);


	free(mig);
	if(ivdisk==0) free(velos);
	if(ivdisk==0) free(etas);

	/* backup to tape */
	if(ibacko==1) {
                fprintf(jpfp," Backup %s , %s and %s to %s \n",
                        diskimg,diskhdr,diskfld,backupo);
                tar3to(backupo,diskimg,diskhdr,diskfld);
        }

	return 0;
}



void migs(int nl, int ns, int nofo, float *mig, int ntau,
	FILE *imgfp, int ktrace, int ntl, float dldt, float *trace, 
	float *tras, float dtl, float tminl, int nt, 
	float *ss, float *sl, float *gs, float *gl, 
	int *imutel, float *fovt2, float *s, float *l,
	float *tm, int *iofs, float apers, float aperl, float *fold,
	float *w1, float *work1, float *work2, float *wsave, float *tracef,
	int *ifcut, int *ltaper, float ksmax, float klmax,
	float f0, float df, int nf, float ftaper, int nsave, int nfft,
	float dtau, float tau0, float angmax, 
	int *indxw,int nindxw, int incdxw,
	float *resamp, int ires, int ntres,
	float *s2, float *scs, float *scg, int ncpu, 
	float  *vr2, float *vi2, float *vq4, float *tau,
	float *etamig, int etatype) {


	int il, itr, i, ir0, it, ii;
	float tmp;

	long lwork, ltmp;
	int nsl;

	nsl = ns*nl;
	/* linearly interpolate input trace */
	for(itr=0;itr<ktrace;itr++) {
		ir0 = itr*nt;
		if(ntl!=nt) {
               		for(it=0;it<ntl;it++) {
                   		tmp = it*dldt;
                       		i = (int)tmp;
                       		tmp = tmp - i;
                   		if(i>=0 && i<nt-1) {
                       			trace[it] = 
						(1.-tmp)*tras[i+ir0]+
						tmp*tras[i+1+ir0];
                       		}
			}	
		} else {
               		for(it=0;it<ntl;it++)
                     		trace[it] = tras[it+ir0];
		}
	
		/* 3d prestack time migration */
		ii = iofs[itr]-1;
		vr22vi2_(vr2+itr*ntau,vi2,&ntau);
		vi22vq4_(vi2,vq4,&ntau);
		vt2s2_(vr2+itr*ntau,vq4,s2,tau,&ntau);

		if(etatype==-1) {
			pstm3d_(trace,&ntl,&tminl,&dtl,&ss[itr],&sl[itr],
			&gs[itr],&gl[itr],&imutel[itr],
			fovt2+itr*ntau,mig+ii*ntau*nsl,&nsl,&ntau,s,l,
			tm,&apers,&aperl,fold+nsl*ii,w1,
			&f0,&df,&nf,&ftaper,
			ifcut,ltaper,tracef,
			wsave,&nsave,work1,work2,&nfft,
			&ksmax,&klmax,&dtau,&tau0,&angmax,
			indxw,&nindxw,&incdxw,&ncpu,
			resamp,&ires,&ntres,s2,scs,scg);
		} else {
			pstmeta_(trace,&ntl,&tminl,&dtl,&ss[itr],&sl[itr],
			&gs[itr],&gl[itr],&imutel[itr],
			fovt2+itr*ntau,mig+ii*ntau*nsl,&nsl,&ntau,s,l,
			tm,&apers,&aperl,fold+nsl*ii,w1,
			&f0,&df,&nf,&ftaper,
			ifcut,ltaper,tracef,
			wsave,&nsave,work1,work2,&nfft,
			&ksmax,&klmax,&dtau,&tau0,&angmax,
			indxw,&nindxw,&incdxw,&ncpu,
			resamp,&ires,&ntres,s2,scs,scg,etamig+itr*ntau);
		}
	}

}