kzmig.c

seismic software,very useful

	if (getparstring("backupo",&backupo)) ibacko = 1;
	if(ibacko==1) isave=1;
	if(ntrdsk>0) isave=1;

	fprintf(jpfp,"\n");
	fprintf(jpfp," -------- KZMIG PRINTOUT -------- \n");
	fprintf(jpfp,"\n");

	/* detect starting trace of migration */
	jtrace = 0;
	if( (hffp = fopen(hisfile,"r"))!=NULL ) {
		fclose(hffp);
                hffp = efopen(hisfile,"r+");
		buf = (char *) malloc(81*sizeof(char));
		do {
                                bzero(buf,81);
                                if(fgets(buf,81,hffp)==NULL) break;
                                if(strncmp(buf,
                                "Number of Traces Processed: ",28)==0) {
                                        sscanf(buf+28,"%d %d",&ntrpre,&jtrace);
                                }
                } while(feof(hffp)==0);
		free(buf);
		fprintf(jpfp," From history file: \n");
		fprintf(jpfp,"Number of Traces Processed: %d\n",ntrpre);
		if(jtrace>0) {
                fprintf(jpfp,"Number of Traces Processed for this input: %d\n",
                        jtrace);
		}
	} else {
		hffp = efopen(hisfile,"w");
		ntrpre = 0;
	}

	/*  compute the reference traveltimes and lateral slowness	*/
	dr = (dds<ddl)?dds:ddl;
	rmax = sqrt((exst-fst)*(exst-fst)+(eyst-flt)*(eyst-flt));
	if(rmax<sqrt((est-fxst)*(est-fxst)+(elt-fyst)*(elt-fyst)))
		rmax = sqrt((est-fxst)*(est-fxst)+(elt-fyst)*(elt-fyst));
	if(rmax>0.5*offmax+sqrt(apers*apers+aperl*aperl)+dxst+dyst)
		rmax = 0.5*offmax+sqrt(apers*apers+aperl*aperl)+dxst+dyst;

	nr = 1+(int)(rmax/dr+0.99);

/*
	fprintf(stderr," before timeb dds=%g ddl=%g \n",dds,ddl);
	fprintf(stderr," before timeb rmax=%d dr=%g nr=%d nzt=%d \n",
	rmax,dr,nr,nzt);
*/

 	tb0 = alloc1float(nzt*nr);
 	pb0 = alloc1float(nzt*nr);


 	timeb_(&nr,&nzt,&dr,&dzt,&fzt,&dvz,&v0,tb0,pb0);
	free1float(pb0);
 	pb = alloc1float(nz*nr);
 	tb = alloc1float(nz*nr);
 	timeb_(&nr,&nz,&dr,&dzt,&fz,&dvz,&v0,tb,pb);


	/* compute memory requirement for migration*/
 
	ntmp = nl*ndl;
	nl = ntmp;
	ntmp = ns*nds;
	ns = ntmp;

	nlcore = nl;

	if(ntrend==0) {
		nsout = ns;
		nlout = nl;
		dsout = ds;
		dlout = dl;
	} else {
		nsout = ntrend;
		nlout = 1;
		nlcore = 1;
		dsout = (send-sstart)/(nsout-1);
		dlout = (lend-lstart)/(nsout-1);
		dend = sqrt ( (send-sstart)*(send-sstart)
		      + (lend-lstart)*(lend-lstart) ) / (nsout-1);
	}

	/* check grid paramters if migration output line is slant	*/
	if(dlt>999990) {
		if(fabs(fst-sstart)<0.1) fst = sstart;
		if(fabs(flt-lstart)<0.1) flt = lstart;
		tmp = sqrt((send-sstart)*(send-sstart)
			+(lend-lstart)*(lend-lstart));
		tmp = tmp/dst+1.5;
		itmp = tmp;
		if(nlt!=1 || fst!=sstart || flt!=lstart || itmp!=nst)
err("grid parameters in traveltime table are inconsistent with output!\n");
 	}
		
	lwork = nsout * nlcore;
	lwork = lwork*nzo*nofo*sizeof(float);

	ltmp = nz*nsout*nlout*nofo+nt+ntl+nt*mtrace;
	ltmp = ltmp+2*nsout*nlout+mtrace*4;
	ltmp = ltmp  + nfft*(nf+2) + ntl*nf + nsave;
	ltmp = ltmp * sizeof(float);
	ltmp = ltmp + (nsout*nlout*nofo+mtrace+nf*2)*sizeof(int);
	ltmp += (nz*nsout*nlout*ntab+nzt*nst*nlt)*sizeof(float);

	ltmp += (nzt*nr+2*nz*nr)*sizeof(float);
	ltmp += 2*nz*nsout*nlout*sizeof(float);
	lwork = lwork + ltmp;
	if(lwork>mlimit) err("mlimit too small; need mlimit=%d \n ",
		(lwork+1024*1024-1)/(1024*1024));

	tmp = mlimit - lwork;
	tmp = tmp / (nz*nsout*nlout*sizeof(float));
	lwork = lwork - nz*nsout*nlout*ntab*sizeof(float); 
	ntab =  ntab + tmp;
	if(ntab>nxst*nyst)  ntab = nxst*nyst; 
	lwork = lwork + nz*nsout*nlout*ntab*sizeof(float); 
	
	fprintf(jpfp," Total Memory used (in Byte) = %g \n", (float)lwork);
	fprintf(jpfp,"\n");

	/* allocation of memory */
	mig = ealloc2float(nsout*nlcore*nzo,nofo);	
	fold = ealloc2float(nsout*nlcore*nz,nofo);	
	off = (int*) emalloc(nsout*nlout*nofo*sizeof(int));	
	tgain = (float*) emalloc(nt*sizeof(float));
	zgain = (float*) emalloc(nzo*sizeof(float));
 	trace = (float*) emalloc(ntl*sizeof(float));	
	tras = (float*) emalloc(nt*mtrace*sizeof(float));	
	imutel = (int*) emalloc(mtrace*sizeof(int));
	s = (float*) emalloc(nsout*nlout*sizeof(float));
	l = (float*) emalloc(nsout*nlout*sizeof(float));
	st = (float*) emalloc(nst*sizeof(float));
	lt = (float*) emalloc(nst*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));
 	work1 = (float*) emalloc(nfft*nf*sizeof(float));
 	sqrtf = (float*) emalloc(nfft*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));
 	ttab = ealloc2float(nz*nsout*nlout,ntab);
	tt = (float*) emalloc(nzt*nst*nlt*sizeof(float));	

	itrace = 0;

	bzero(fold[0],nofo*nsout*nlout*nz*sizeof(float));
	bzero(mig[0],nofo*nsout*nlcore*nzo*sizeof(float));
	if(ntrend>0) {
		for(io=0;io<nsout*nofo;io++) disend[io] =  99999999999.;
		dismax = disaper(apers,aperl,sstart,lstart,send,lend);
		xy2abc(sstart,lstart,send,lend,&aa,&bb,&cc);
	}

	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 {
		for(is=0;is<ntrend;is++) {
			s[is] = sstart+is*(send-sstart)/(ntrend-1.);
			l[is] = lstart+is*(lend-lstart)/(ntrend-1.);
		}
		for(is=0;is<nst;is++) {
			st[is] = sstart+is*(send-sstart)/(nst-1.);
			lt[is] = lstart+is*(lend-lstart)/(nst-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];
	}

/*	check migration output range	*/
	stmin = fst; if (est<fst) stmin = est;
        ltmin = flt; if (elt<flt) ltmin = elt;
        stmax = est; if (est<fst) stmax = fst;
        ltmax = elt; if (elt<flt) ltmax = flt;
        if(stmin>sinmin || stmax<sinmax-0.01
        || ltmin>linmin || ltmax<linmax-0.01) {
                warn("stmin=%g > sinmin=%g ? \n",stmin,sinmin);
                warn("stmax=%g < sinmax-0.01=%g ? \n",stmax,sinmax-0.01);
                warn("ltmin=%g > linmin=%g ? \n",ltmin,linmin);
                warn("ltmax=%g < linmax-0.01=%g ? \n",ltmax,linmax-0.01);
                err(" migration output is out of time tables!\n");
        }

	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 = nzo*nofo*nlout*nsout;
		efread(mig[0],sizeof(float),ltmp,imgfp);
		efclose(imgfp);
		fprintf(jpfp," diskimg file opened  \n");

		fprintf(jpfp," open diskfld file... \n");
		fldfp = efopen(diskfld,"r+");
		fseek2g(fldfp,0,1);
		efread(fold[0],sizeof(float),nofo*nsout*nlout*nz,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.ungpow;
		}
		ipre = 1;
	} else if(isave==1) {
		imgfp = efopen(diskimg,"w+r");
		fseek2g(imgfp,0,1);
		bzero(mig[0],nofo*nsout*nlout*nzo*sizeof(float));
		efwrite(mig[0],sizeof(float),nofo*nsout*nlout*nzo,imgfp);
 		fprintf(jpfp," initialize diskimg file... \n");
 		fprintf(jpfp," initialization of diskimg file done \n");
		efclose(imgfp);

		fprintf(jpfp," initialize diskfld file... \n");
		fldfp = efopen(diskfld,"w+r");
		fseek2g(fldfp,0,1);
		bzero(fold[0],nofo*nsout*nlout*nz*sizeof(float));
		efwrite(fold[0],sizeof(float),nofo*nsout*nlout*nz,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.ungpow = 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;
	m1000 = 0;
	oldheadx = newheadx = 0;
	oldheady = newheady = 0;
	ssmin = slmin = 99999999.0;
	ssmax = slmax = -99999999.0;
 
	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.;
			}
		}
 	} else {
		for(it=0;it<nt;it++) tgain[it] = 1.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);
	if(ntrend>0) fprintf(jpfp," dismax=%g aa=%g bb=%g cc=%g \n",
			dismax,aa,bb,cc);
	fprintf(jpfp," fzo=%g nzo=%d dzo=%g\n",fzo,nzo,dzo);
	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 ntl=%d dt=%g t0=%g nfft=%d\n",
		ntrpre,nt,ntl,dt,t0,nfft);
 	fprintf(jpfp," tpow=%g zpow=%g ofomin=%g dofo=%g nofo=%d isave=%d \n",
		tpow,zpow,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 ntrdsk=%d \n",
		nlcore,mtrace,ntras,mlimit/(1024*1024),ntrdsk);
	fprintf(jpfp," f0=%g df=%g nf=%d ftaper=%g ksmax=%g klmax=%g\n",
		f0,df,nf,ftaper,ksmax,klmax);
	fprintf(jpfp," offmax=%g ntab=%d \n",offmax,ntab);
 	fprintf(jpfp," fzt=%g dzt=%g nzt=%d \n",fzt,dzt,nzt);
	fprintf(jpfp," fst=%g dst=%g nst=%d \n",fst,dst,nst);
	fprintf(jpfp," flt=%g dlt=%g nlt=%d \n",flt,dlt,nlt);
	fprintf(jpfp," fxst=%g dxst=%g nxst=%d \n",fxst,dxst,nxst);
	fprintf(jpfp," fyst=%g dyst=%g nyst=%d \n",fyst,dyst,nyst);
	fprintf(jpfp," v0=%g dvz=%g gath=%d \n",v0,dvz,gath);
	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," \n");
	

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

 	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) {
		jtrace = ntrpre - jtrace;