kzmig.c

seismic software,very useful

	char *ttfile,*datain,*dataout;
	char *diskfld="DISKFOLD", *hisfile="KZMIG.HISFILE";
	int ibacko,isave, ihis=0, ipre=0;
	float ds, dl, tmp;
	int mtrace, ntrpre=0, ktrace, itrace, jtrace, ntrace, nofo;
	int ltrace;
	int *imutel, m1000, ifact, it0;
	float **mig, **fold, *tgain, *zgain,*s, *l, *trace, *tras;
 	float *ofo;
	int *off, itgain, izgain, nl, ns, ntras, is, il, io;
	int nsread, nlread;
	int mlimit, nlcore, *iofs; 
	float tpow, zpow, sstart, lstart, t0;
	float xs, ys, xg, yg, dt, dzo, fzo, v0, dvz; 
	float *ss, *sl, *gs, *gl, ms, ml, ofomin, dofo; 
	float *st, *lt;
	char *buf;
	int itmp, itmp2, ntl;
	float tminl, dtl, dldt, apers, aperl;
 	float f0, df, ftaper;
	int nf, nfft, *ifcut, *ltaper, nsave, tahd;
	float *work1, *work2, *tracef, *wsave, *sqrtf;
	float ddl, dds, ksmax, klmax;
	int ndl, nds, ntmp, tapecntl;
	float rmsmar, ampmax, amprms;
	char *badtracefile="BAD_TRACE_FILE";
	FILE *btfp;
	float dd, mx, my, ofs;
	int ids, idl; 
	int cdp1, cdp2, cdp3;
	int cdppds, cdppdds, cdppdl, cdppddl, cdpnum;
	int nss, nll, iix, iiy, is0, il0;
	float sinmin, sinmax, linmin, linmax; 
	float angmax; 
	char *jpfile;
        FILE *jpfp;
	int ntrend, nsout, nlout, isend;
	float send, lend, dsout, dlout, dend, *disend;
	float ps, pl;
	int iofoin=0;
	float ofol, ofor;
	int offnear, offfar, flexbin;
	int one=1; 
	int range,gath,oldheadx,oldheady,newheadx,newheady,gottrace,ready,done,
	    nr,nseek,ntab,itab,nxs0,nys0,nxs,nys,iz,ixs,iys,izi,izw,
	    nxsb,nysb,nxs0b,nys0b;
	float **ttab,*tb0,*pb0,*tb,*pb,*tt;
	float offmax,ssmin,ssmax,slmin,slmax,rmax,dr,zi,fz,fxs,fys,scale,fold0;	
	
	usghed ugh;
	int nzt,nst,nlt,nxst,nyst,ierr;
	float fzt,fst,flt,fxst,fyst,exst,eyst,dzt,dst,dlt,dxst,dyst,ezt,est,elt;
	float aa,bb,cc,dismax,tabx,taby;
	float ltmin, ltmax, stmin, stmax;
 	
	long lwork, ltmp;
	long lpos;

	int ntrdsk, ndskwt;

	int ncpu;
	char *envs;

    	/* get parameters */
    	initargs(argc,argv);
	requestdoc(1);

	/* open input and output data sets */
    	if (!getparstring("datain", &datain)) {
		infp = stdin;
	} else {
    		infp = efopen(datain,"r");
	}
    	if (!getparstring("dataout", &dataout)) {
		outfp = stdout;
	} else {
    		outfp = efopen(dataout,"w");
	}
	fseek2g(infp,0,1);
        fseek2g(outfp,0,1);

	/* required parameters */
    	if (!getparfloat("x1",&x1)) err("must specify x1");
    	if (!getparfloat("y1",&y1)) err("must specify y1");
    	if (!getparfloat("s1",&s1)) err("must specify s1");
    	if (!getparfloat("l1",&l1)) err("must specify l1");
    	if (!getparint("cdp1",&cdp1)) err("must specify cdp1");
    	if (!getparfloat("x2",&x2)) err("must specify x2");
    	if (!getparfloat("y2",&y2)) err("must specify y2");
    	if (!getparfloat("s2",&s2)) err("must specify s2");
    	if (!getparfloat("l2",&l2)) err("must specify l2");
    	if (!getparint("cdp2",&cdp2)) err("must specify cdp2");
    	if (!getparfloat("x3",&x3)) err("must specify x3");
    	if (!getparfloat("y3",&y3)) err("must specify y3");
    	if (!getparfloat("s3",&s3)) err("must specify s3");
    	if (!getparfloat("l3",&l3)) err("must specify l3");
    	if (!getparint("cdp3",&cdp3)) err("must specify cdp3");
    	if (!getparfloat("dds",&dds)) err("Must specify dds!\n");
    	if (!getparfloat("ddl",&ddl)) ddl = 0.;
 
	/* optional parameters */
    	if (!getparfloat("sstart",&sstart)) sstart = s1;
    	if (!getparfloat("lstart",&lstart)) lstart = l1;
    	if (!getparfloat("ds",&ds)) ds = dds;
    	if (!getparfloat("dl",&dl)) dl = ddl;
    	if (!getparint("ndl",&ndl)) ndl = 1; if(ndl<1) ndl=1;
    	if (!getparint("nds",&nds)) nds = 1; if(nds<1) nds=1;
    	if (!getparfloat("angmax",&angmax)) angmax = 60.;

	if(ds<dds) err(" dds larger than ds ");
	if(dl<ddl) err(" ddl larger than dl ");

	if(dds>0.) {
		tmp = (s2-s1)/dds + 1.5;
		ncdppl = tmp;
	} else {
		ncdppl = 1;
 	}

	if(ddl!=0.) {
		tmp = (l3-l1)/ddl + 1.5;
		nlines = tmp;
	} else {
		nlines = 1;
	}
 	if(nlines<10) {
		tahd = 1;
		fold0 = 10.;
	}
	else {
		tahd = 0;
		fold0 = 100.;
	}
    	if (!getparint("nl",&nlread)) nlread = nlines;
    	if (!getparint("ns",&nsread)) nsread = ncdppl;
	nl = nlread;
	ns = nsread;
    	if (dl==0. && nl>1) err("Must specify nonzero dl for migration");
    	if (ds==0. && ns>1) err("Must specify nonzeor ds for migration");

	if(dds==0.) { dds = 1.0; ds = dds; }
	if(ddl==0.) { ddl = 1.0; dl = ddl; }

	if(ncdppl>1) {
		tmp = (cdp2 - cdp1)/(ncdppl-1.) + .5;
		cdppdds = tmp;
		tmp = (cdp2 - cdp1)/(ncdppl-1.) + .5;
		cdppds = tmp;
	} else {
		cdppdds = nlines;
		cdppds = nlines;
	}
	if(cdppdds<1) cdppdds = 1;
	if(cdppds<1) cdppds = 1;
	if(nlines>1) {
		tmp = (cdp3 - cdp1)/(nlines-1.) + .5;
		cdppddl = tmp;
		tmp = (cdp3 - cdp1)/(nlines-1.) + .5;
		cdppdl = tmp;
	} else {
		cdppddl = ncdppl;
		cdppdl = ncdppl;
	}
	if(cdppddl<1) cdppddl = 1;
	if(cdppdl<1) cdppdl = 1;

	if(cdppdds!=1 && cdppddl!=1 ) {
		warn("cdppdds=%d cdppddl=%d ",cdppdds,cdppddl);
		err("check cdp1, cdp2, cdp3 ");
	}
	if(cdppdds==1) {
		cdpnum = 0;
	} else {
		cdpnum = 1;
	}


	tmp = ds/dds+0.5;
	nss  = tmp;
	tmp = dl/ddl+0.5;
	nll  = tmp;
	tmp = (sstart-s1)/dds+0.5;
	is0 = tmp;
	tmp = (lstart-l1)/ddl+0.5;
	il0 = tmp;


    	if (!getparint("mtrace",&mtrace)) mtrace = 1000;
    	if (!getparint("ntras",&ntras)) ntras = 2000000000;
	

    	if (!getparfloat("tpow",&tpow)) tpow = -1.;
    	if (!getparfloat("zpow",&zpow)) zpow = 0.;
	itgain = 0;
	izgain = 0;
	if(tpow!=0.) itgain = 1;
	if(zpow!=0.) izgain = 1;

    	if (!getparint("ifact",&ifact)) ifact = 1;
    	if (!getparfloat("apers",&apers)) apers = fabs(ns*ds)/4.;
    	if (!getparfloat("aperl",&aperl)) aperl = fabs(nl*dl)/4.;

    	if (!getparfloat("rmsmar",&rmsmar)) rmsmar = 50.;
	getparstring("badtracefile",&badtracefile);
	if( (btfp = fopen(badtracefile,"r"))!=NULL ) {
		efclose(btfp);	
		btfp = efopen(badtracefile,"a");
	} else {
		btfp = efopen(badtracefile,"w");
	}

    	if ( getparfloat("send",&send) 
    	  && getparfloat("lend",&lend) 
    	  && getparint("ntrend",&ntrend) ) {
		if(ntrend<1) err("ntrend must be greater than 0 ");
		nsout = ntrend;
		nlout = 1;
	} else {
		ntrend = 0;
		nsout = ns;
		nlout = nl;
	} 

	if(!getparint("ntrdsk",&ntrdsk)) ntrdsk = 0;
        ndskwt = 0;

	if (!getparint("ncpu",&ncpu)) ncpu = 1;
	envs = (char*) emalloc(80*sizeof(char));
        if(!getenv("MP_SET_NUMTHREADS")) {
                sprintf(envs,"%s=%d","MP_SET_NUMTHREADS",ncpu);
                putenv(envs);
        }

	/* read id headers */
        fgethdr(infp,&ch,&bh);

	/* read in first trace for nt and dt */
    	if (!fgettr(infp,&tra))  err("can't get first trace");
    	nt = tra.ns;
    	dt = (float)tra.dt/1000000.;
	t0 = (float)tr.delrt/1000;

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

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

	tabx = offmax/dxst;
	taby = offmax/dyst;
	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;
	if(fzo<fzt || fzo+(nzo-1)*dzo>ezt) 
	  err("migration depth is beyond the depth in traveltime table!\n");
	nz0 = (int)((fzo-fzt)/dzt);
	fz = fzt+nz0*dzt;
	nz = 1+(int)((fzo+(nzo-1)*dzo-fzt)/dzt+0.99)-nz0;
 
    	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 (!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 = 10.;
    	if (!getparint("nf",&nf)) nf = 6;
    	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(nofo>1) {
		bh.tsort = 2;
	} else {
		bh.tsort = 4;
	}
 	fputhdr(outfp,&ch,&bh);


    	if (!getparint("mlimit",&mlimit)) mlimit = 256;
    	mlimit = mlimit * 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;