kzmig.c

seismic software,very useful

	}
   	nt = tra.ns;
   	dt = (float)tra.dt/1000000.;
	t0 = (float)tr.delrt/1000;

	if(outputonly==1) idhdrs(&ch,&bh,nt);

	ntotal = tra.tracl;

	ntl = nt * ifact;
	dtl = dt/ifact;
	tminl = t0;
	dldt = dtl/dt;

	/* get traveltime grid header info */
	if(!getparstring("ttfile",&ttfile))  ttfile="ttfile";
	ttfp = efopen(ttfile,"r");
	fseek2g(ttfp,0,1);
	ierr = fgetusghdr(ttfp,&ugh);
	if (ierr!=0) err("Grid parameters of %s required!\n",ttfile);

	nzt = ugh.n1;	fzt = ugh.o1;	dzt = ugh.d1;
	nst = ugh.n2;	fst = ugh.o2;	dst = ugh.d2;
	nlt = ugh.n3;	flt = ugh.o3;	dlt = ugh.d3;
	nxst = ugh.n4;	fxst = ugh.o4;	dxst = ugh.d4;
	nyst = ugh.n5;	fyst = ugh.o5;	dyst = ugh.d5;
	if (dyst==0.) dyst = 1.;
	if (dxst==0.) dxst = 1.;
	if (dst==0.) dst = 1.;
	if (dlt==0.) dlt = 1.;

	ezt = fzt+(nzt-1)*dzt;
	if(dlt<999990){
		est = fst+(nst-1)*dst;
		elt = flt+(nlt-1)*dlt;
	} else{
		est = send;
		elt = lend;
	}
	exst = fxst+(nxst-1)*dxst;
	eyst = fyst+(nyst-1)*dyst;
	if(!getparint("izw",&izw)) izw = (fzt<4*dzt)? 4-fzt/dzt: 0;

   	if (!getparfloat("v0",&v0)) v0 = 1500.;
   	if (!getparfloat("dvz",&dvz)) dvz = 0.2;

   	if (!getparfloat("offmax",&offmax)) offmax = 6000;

	tabx = offmax/dxst;
	if(tabx>nxst) tabx = nxst;
	taby = offmax/dyst;
	if(taby>nyst) taby = nyst;
	ntab = 0.5*tabx*taby+9.+3.*sqrt(tabx*tabx+taby*taby);
	if(ntab>(tabx+3)*nyst) ntab = (tabx+3)*nyst;
	if(ntab>(taby+3)*nxst) ntab = (taby+3)*nxst;
	if(ntab>nxst*nyst) ntab = nxst*nyst;
   
   	if (!getparint("nzo",&nzo)) nzo = 10*(nzt-1)+1;
   	if (!getparfloat("dzo",&dzo)) dzo = 0.1*dzt;
   	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	/*
	fprintf(stderr,"fzo=%g fzt=%g fzo+(nzo-1)*dzo=%g ezt=%g \n",
		fzo,fzt,fzo+(nzo-1)*dzo,ezt);
	*/

	if(fabs(fzo-fzt)<0.1*dzo) fzo = fzt;
	if(fabs(fzo+(nzo-1)*dzo-ezt)<0.1*dzo) ezt  = fzo+(nzo-1)*dzo;

	if(fzo<fzt || fzo+(nzo-1)*dzo>ezt) {
	  warn("migration depth fzo=%g fze=%g \n", fzo, fzo+(nzo-1)*dzo);
	  warn("traveltime depth fzt=%g ezt=%g \n", fzt, ezt);
	  err("migration depth is beyond the depth in traveltime table!\n");
	}
	if(getparstring("hzgrid",&hzgrid)) {
		ihzgrid = 1;
	} else {
		ihzgrid = 0;
	}

	if(ihzgrid==0) {
		nz0 = (int)((fzo-fzt)/dzt);
		fz = fzt+nz0*dzt;
		nz = 1+(int)((fzo+(nzo-1)*dzo-fzt)/dzt+0.99)-nz0;
	} else {
		fzo = fzt;
		nz0 = 0;
		fz = fzt;
		nz = 1+nzt;
	}
 
   	if (!getparint("offnear",&offnear)) offnear = 0;
   	if (!getparint("offfar",&offfar)) offfar = 0;
   	if (!getparint("flexbin",&flexbin)) flexbin = 0;

   	if (!getparfloat("oofo",&ofomin)) ofomin = 0.;
   	if (!getparfloat("dofo",&dofo)) dofo = 999999.;
   	if (!getparfloat("obtol",&obtol)) obtol = 0.5*dofo;
   	if (!getparint("nofo",&nofo)) nofo = 1;
	itmp = countparname("ofo");
	if(itmp>0 && itmp!=nofo) err("number of ofo not match with nofo");
	ofo = (float*) emalloc(nofo*sizeof(float));
	if(itmp>0) {
		getparfloat("ofo",ofo);
		if(nofo>1) {
			ofol = ofo[0] - 0.5*(ofo[1]-ofo[0]);
			ofor = ofo[nofo-1] + 0.5*(ofo[nofo-1]-ofo[nofo-2]);
		} else {
			ofol = ofo[0] - 0.5*dofo;
			ofor = ofo[nofo-1] + 0.5*dofo;
			ofomin = ofo[0];
		}
	} else {
		for(io=0;io<nofo;io++) ofo[io] = ofomin + io*dofo;
	}
	iofoin = itmp;
	
   	if (!getparfloat("f0",&f0)) f0 = 10.;
   	if (!getparfloat("df",&df)) df = 5.;
   	if (!getparint("nf",&nf)) nf = 12;
   	if (!getparfloat("ftaper",&ftaper)) ftaper = 5.;

	/* compute maximum wavenumbers to be used in migration */
   	if (!getparfloat("ksmax",&ksmax)) {
		if(ncdppl>1) {
			ksmax = 0.5/dds;
		} else { 
			ksmax = 0.5/0.001;
		}
	}
   	if (!getparfloat("klmax",&klmax)) {
		if(nlines>1) {
			klmax = 0.5/ddl;
		} else {
			klmax = 0.5/0.001;
		}
	}

   	if (!getparint("tapecntl",&tapecntl)) tapecntl = 0;
	
	nsave = ((ntl * 3 / 2)/2)*2;
	radix_(&nsave,&nfft);
	nsave = nfft * 2 + 30;

	/* update id headers and write to output */
	bh.hns = nzo;
	bh.fold = nofo;
	if(dzo>32.) {
		bh.hdt = dzo;
	} else {
		bh.hdt = dzo*1000;
	}
	if(nofo>1) {
		bh.tsort = 2;
	} else {
		bh.tsort = 4;
	}
 	if(traceout==1) fputhdr(outfp,&ch,&bh);


   	if (!getparint("mlimit",&mlimit)) mlimit = 256;
   	llimit = mlimit;
	llimit = llimit * 1024 * 1024;

	if (!getparint("gath",&gath)) gath = 0;

	getparstring("diskhdr",&diskhdr);
	getparstring("diskimg",&diskimg);
	getparstring("diskfld",&diskfld);
	if(getparstring("hisfile",&hisfile)) ihis=1;

	if(!getparstring("jpfile",&jpfile)) {
                jpfp = stderr;
        } else {
                jpfp = efopen(jpfile,"w");
        }


	if (!getparint("isave",&isave)) isave = 0;

	if (getparstring("backupi",&backupi)) {
		fprintf(jpfp," Backup from the previous run ");
                fprintf(jpfp," Backup %s , %s and %s from %s \n",
                        diskimg,diskhdr,diskfld,backupi);
		tar3fr(backupi,diskimg,diskhdr,diskfld);
	}

	ibacko = 0;
	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);
	}

	/* get hzgrid if specified */
	fzout = (float*) malloc(nsout*nlout*sizeof(float));
	if(ihzgrid==1) {
		hzfp = efopen(hzgrid,"r");
		ierr = fgetusghdr(hzfp,&hzgh);
		if (ierr!=0) err("Grid parameters of %s required!\n",hzgrid);
		nshz = hzgh.n1;
		nlhz = hzgh.n2;
		oshz = hzgh.o1;
		olhz = hzgh.o2;
		dshz = hzgh.d1;
		dlhz = hzgh.d2;
		hz = (float*) malloc(nshz*nlhz*sizeof(float));
		fseek(hzfp,0,0);
		efread(hz,sizeof(float),nshz*nlhz,hzfp);
		fclose(hzfp);
	} else {
		for(ix=0;ix<nlout*nsout;ix++) fzout[ix] = fzo;
	}

	/* 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;
	tmp = (lwork+1024*1024-1)/(1024*1024); 

	if(lwork>llimit) err("mlimit too small; need mlimit=%d \n ",(int)tmp);

	tmp = llimit - 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));	

	datetime=(char*) emalloc(28*sizeof(char));

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