ktmig.c

seismic software,very useful

	ltmp = ltmp + ntau + nfft*(nf+1) + ntl*nf + nsave;
	ltmp = ltmp * sizeof(float);
	ltmp = ltmp + (nsout*nlout*nofo+mtrace+nf*2)*sizeof(int);
	ltmp = ltmp + ntres*nf*sizeof(float);
	ltmp = ltmp + nlout*nindxw*sizeof(int);
	lwork = lwork + ltmp;


	if(lwork>mlimit) err("mlimit too small ");
	
	fprintf(jpfp," Total Memory used (in Byte) = %g \n", (float)lwork);
	fprintf(jpfp,"\n");

	/* allocation of memory */
	mig = (float*) emalloc(nsout*nlcore*nofo*ntau*sizeof(float));	
	fold = (float*) emalloc(nsout*nlout*nofo*sizeof(float));	
	off = (int*) emalloc(nsout*nlout*nofo*sizeof(int));	
	tgain = (float*) emalloc(nt*sizeof(float));
	tugain = (float*) emalloc(ntau*sizeof(float));
	trace = (float*) emalloc(ntl*sizeof(float));	
	tras = (float*) emalloc(nt*mtrace*sizeof(float));	
	imutel = (int*) emalloc(mtrace*sizeof(int));
	tm = (float*) emalloc(ntau*sizeof(float));
	tmig = (float*) emalloc(ntau*sizeof(float));
	fovt2 = (float*) emalloc(ntau*mtrace*sizeof(float));
	s = (float*) emalloc(nsout*nlout*sizeof(float));
	l = (float*) emalloc(nsout*nlout*sizeof(float));
	vel = (float*) emalloc(ntv*sizeof(float));
	ss = (float*) emalloc(mtrace*sizeof(float));
	sl = (float*) emalloc(mtrace*sizeof(float));
	gs = (float*) emalloc(mtrace*sizeof(float));
	gl = (float*) emalloc(mtrace*sizeof(float));
	iofs = (int*) emalloc(mtrace*sizeof(int));
	w1 = (float*) emalloc(ntau*sizeof(float));
	work1 = (float*) emalloc(nfft*nf*sizeof(float));
	tracef = (float*) emalloc(ntl*nf*sizeof(float));
	work2 = (float*) emalloc(nfft*sizeof(float));
	wsave = (float*) emalloc(nsave*sizeof(float));
	ifcut = (int*) emalloc(nf*sizeof(int));
	ltaper = (int*) emalloc(nf*sizeof(int));
	disend = (float*) emalloc(nsout*nofo*sizeof(float));
	resamp = (float*) emalloc(ntres*nf*sizeof(float));
	indxw = (int*) emalloc(nindxw*sizeof(int));

	ss2 = (float*) emalloc(ntau*sizeof(float));
	scg = (float*) emalloc(ntau*sizeof(float));
	scs = (float*) emalloc(ntau*sizeof(float));
	vr2 = (float*) emalloc(ntau*mtrace*sizeof(float));
	vi2 = (float*) emalloc(ntau*sizeof(float));
	tau = (float*) emalloc(ntau*sizeof(float));
	vq4 = (float*) emalloc(ntau*sizeof(float));

	itrace = 0;

	bzero(fold,nofo*nsout*nlout*sizeof(float));
	bzero(mig,nofo*nsout*nlcore*ntau*sizeof(float));
	if(ntrend>0) {
		for(io=0;io<nsout*nofo;io++) disend[io] =  99999999999.;
	}

	if(ntrend==0) { 
		for(il=0;il<nlout;il++) {
			for(is=0;is<nsout;is++) {
				s[is+il*nsout] = sstart+is/nds*ds+is%nds*dds;
			}
		}
		for(il=0;il<nlout;il++) {
			for(is=0;is<nsout;is++) {
				l[is+il*nsout] = lstart+il/ndl*dl+il%ndl*ddl;
			}
		}
	} else {
		fprintf(jpfp,
		" Output s and l Locations of Target Line\n");
		for(is=0;is<ntrend;is++) {
			s[is] = sstart+is*(send-sstart)/(ntrend-1.);
			l[is] = lstart+is*(lend-lstart)/(ntrend-1.);
			fprintf(jpfp,
				"   s=%g l=%g trace=%d \n",s[is],l[is],is+1); 
		}
	}

	sinmin = s[0];
	sinmax = s[0];
	linmin = l[0];
	linmax = l[0];

	for(is=0;is<nsout*nlout;is++) {
		if(sinmin>s[is]) sinmin = s[is];
		if(sinmax<s[is]) sinmax = s[is];
		if(linmin>l[is]) linmin = l[is];
		if(linmax<l[is]) linmax = l[is];
	}
		
	sinmin = sinmin - apers;
	sinmax = sinmax + apers;
	linmin = linmin - aperl;
	linmax = linmax + aperl;

	/* allocate disk space for image */
	if(ntrpre>0) {

		fprintf(jpfp," open diskimg file... \n");
		imgfp = efopen(diskimg,"r+");
		fseek2g(imgfp,0,1);
		ltmp = ntau*nofo*nlout*nsout;
		efread(mig,sizeof(float),ltmp,imgfp);
		fprintf(jpfp," diskimg file opened  \n");

		fprintf(jpfp," open diskfld file... \n");
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		efread(fold,sizeof(float),nofo*nsout*nlout,fldfp);
		efclose(fldfp);
		fprintf(jpfp," diskfld file opened  \n");

		fprintf(jpfp," open diskhdr file... \n");
		hdrfp = efopen(diskhdr,"r+");
		fseek2g(hdrfp,0,1);
		fprintf(jpfp," diskhdr file opened  \n");
		
		for(ix=0;ix<nsout*nlout*nofo;ix++) {
			efread(&tr,sizeof(char),240,hdrfp);
			off[ix] = tr.offset;
			if(ntrend>0 ) disend[ix] = tr.dz;
		}
		ipre = 1;
	} else if(isave==1) {
		imgfp = efopen(diskimg,"w+r");
		fseek2g(imgfp,0,1);
		bzero(tm,ntau*sizeof(float));
		fprintf(jpfp," initialize diskimg file... \n");
		for(ix=0;ix<nsout*nlout*nofo;ix++) 
			efwrite(tm,sizeof(float),ntau,imgfp);
		fprintf(jpfp," initialization of diskimg file done \n");
		if(nlcore==nlout) efclose(imgfp);

		fprintf(jpfp," initialize diskfld file... \n");
		fldfp = efopen(diskfld,"w+r");
		fseek2g(fldfp,0,1);
		bzero(fold,nofo*nsout*nlout*sizeof(float));
		efwrite(fold,sizeof(float),nofo*nsout*nlout,fldfp);
		efclose(fldfp);
		fprintf(jpfp," initialization of diskfld file done  \n");
	}

	if(ntrpre==0) {
		fprintf(jpfp," initializing diskhdr file... \n");
		hdrfp = efopen(diskhdr,"w");
		bzero(&tr,240);
		tr.offset = 99999999;
		tr.ns = tra.ns;
		tr.dt = tra.dt;
		tr.trid = 2;
		if(ntrend>0) tr.dz = 99999999999.;
		for(iy=0;iy<nlout;iy++) {
			tr.tracr = iy + 1;
			for(ix=0;ix<nsout;ix++) {
				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) {
					tr.cdp = iiy*ncdppl + iix + cdp1;
				} else {
					tr.cdp = iix*nlines + iiy + cdp1;
				}
				tr.tracl = ix + 1;
				for(io=0;io<nofo;io++) {
					tr.cdpt  = io+1;
					efwrite(&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;
	cdppre = tra.cdp - 1;
	m1000 = 0;
	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);
	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);


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


		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;
		} else {
			if(tmp<ofol) {
				io = -1;
			} else if(tmp>ofor) {
				io = nofo;