kzmig.c.9.2.99

seismic software,very useful

void timeb_(int *nr, int *nz, float *dr, float *dz, float *fz, float *a,
	float *v0, float *t, float *p);
void resit_(float *t, int *nx, int *ny, int *nz, int *nz0, int *nr,
	float *dx, float *dy, float *dr, float *fx, float *fy, 
	float *x0, float *y0, float *tb);
void resit0_(float *t, int *nout, int *nz, int *nz0, int *nr, float *dr,
	float *x, float *y, float *x0, float *y0, float *tb);
void latint_(int *nx, int *ny, int *nz, float *x, float *y, float *t,
	int *n1, int *n2, int *nzt, int *n0, float *f1, float *f2, 
	float *d1, float *d2, float *tt);
void bisear_(int *nxin, int *nxout, float *xin, float *xout, int *indx);
void blinint_(int *nx, int *ny, int *nz, float *sx, float *sy, float *t0,
    float *t1, float *t2, float *t3, float *t);

float disp2l(float a, float b, float c,float x, float y);
void xy2abc(float x1, float y1, float x2, float y2,
        float *a, float *b, float *c);
float disaper(float aperx,float apery,float x1,float y1,float x2,float y2);

segytrace tra, tr;
segybhdr bh;
segychdr ch;


main(int argc, char **argv)
{
	int nlines, ncdppl;
   	int nt,nzo,nz,nz0,it,ix,iy;
    	FILE *infp,*outfp,*ttfp,*imgfp,*hdrfp,*hffp,*fldfp;
	float x1,x2,x3,y1,y2,y3,s1,s2,s3,l1,l2,l3;
	char *backupi="", *backupo="";
	char *diskimg="DISKIMAGE", *diskhdr="DISKHEADER";
	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,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 long lwork, ltmp;
	long long lpos, nseek;

	float tr1,tr2,tr3,ln1,ln2,ln3,fline,ftrace;
	String tracekey="tracl", linekey="tracr", trktype, lnktype;
	Value trkval, lnkval;
	int indxtrk, indxlnk;
	int intline, inttrace, ikey;


	int ntrdsk, ndskwt;

	int ncpu, traceout;
	char *envs;
	char *datetime;
	int on2trace, ntotal;
	int outputonly;

   	/* 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 = 1.;
	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("PARALLEL")) {
                sprintf(envs,"%s=%d","PARALLEL",ncpu);
                putenv(envs);
        }
	if (!getparint("traceout",&traceout)) traceout = 1;
	if (!getparint("on2trace",&on2trace)) on2trace = 0;

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

	ikey = 0;
	if ( getparstring("tracekey",&tracekey)
		&& getparstring("linekey",&linekey) ) {

        ikey = 1;
		trktype = hdtype(tracekey);
		lnktype = hdtype(linekey);
		indxtrk = getindex(tracekey);
		indxlnk = getindex(linekey);
		if(!getparfloat("tr1",&tr1)) err(" tr1 missing");
		if(!getparfloat("tr2",&tr2)) err(" tr2 missing");
		if(!getparfloat("tr3",&tr3)) err(" tr3 missing");
		if(!getparfloat("ln1",&ln1)) err(" ln1 missing");
		if(!getparfloat("ln2",&ln2)) err(" ln2 missing");
		if(!getparfloat("ln3",&ln3)) err(" ln3 missing");
	}

	/* read id headers */
    if(outputonly==0) {
		fgethdr(infp,&ch,&bh);
	}

	/* read in first trace for nt and dt */
	if(outputonly==0) {
   		if (!fgettr(infp,&tra)) err("can't get first trace");
	} else {
			getparstring("diskhdr",&diskhdr);
			hdrfp = efopen(diskhdr,"r");
			fseek2g(hdrfp,0,1);
			if(fread(&tra,sizeof(char),240,hdrfp)!=240)
				err(" read error on %s \n",diskhdr);
			efclose(hdrfp);
	}
   	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 = fopen(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 = 5000.;
   	if (!getparfloat("dvz",&dvz)) dvz = 0.6;

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

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