kzmig.c

seismic software,very useful

" sstart+(ns-1)*ds (ssn), ssn+dds,      ...,       ssn+(nds-1)*dds      ",
"       the output crossline locations (l):				",
" lstart      (ll1),      ll1+ddl,      ...,       ll1+(ndl-1)*ddl,	",
" lstart+dl   (ll2),      ll2+ddl,      ...,       ll2+(ndl-1)*ddl 	",
" lstart+2*dl (ll3),      ll3+ddl,      ...,       ll3+(ndl-1)*ddl 	",
"    ...       								",
" lstart+(nl-1)*ds (lln), lln+ddl,      ...,       lln+(ndl-1)*ddl      ",
" 9. The primary memory requirement is 					",
"        ns*nds*nl*ndl*(nofo*nzo+nzt*ntab)*4 (byte)				", 
"    ntab = .5*(offmax/dxst)*(offmax/dyst)		", 
"         + 3.*sqrt( (offmax/dxst)**2 + (offmax/dyst)**2 )  ",
" 10. when nlines is less than 10, a 2.5-D rho filter is applied to 	",
"     compensate for 2D data input					",
" ",
NULL};
/* "AUTHOR:		Z. Liu and Z. Li        9/94   			", */

void migs(int nl, int ns, float **mig, int nz,
	int ktrace, int ntl, float dldt, float *trace, 
	float *tras, float dtl, float tminl, int nt, 
	float *ss, float *sl, float *gs, float *gl, 
	int *imutel, float *s, float *l,
	int *iofs, float apers, float aperl, float **fold,
	float *work1, float *work2, float *wsave, float *tracef, float *sqrtf,
	int tahd, int *ifcut, int *ltaper, float ksmax, float klmax,
	float f0, float df, int nf, float ftaper, int nsave, int nfft,
	float dz, float fz, float angmax,
	int nxs, int nys, float fxs, float fys, float dxs, float dys,
	int nzo, float *fzo, float dzo,
	int nr, float dr, float *tb, float *pb, float **ttab, int igh); 

void pszm3d_(float *trace,int *nt,float *t0,float *dt,float *sx,float *sy,
       float *gx,float *gy,int *imute,float *mig,float *x,float *y,
       int *nxy,int *nzo,float *dzo,float *fzo,float *aperx,float *apery,
       float *fold,float *f0,float *df,int *nf,float *ftaper,
       int *ifcut,int *ltaper,float *tracef,float *wsave,int *nsave,
       float *work1,float *work2,float *sqrtf,int *tahd,int *nfft,
       float *kxmax,float *kymax,float *angmax,float *ts,float *tg,
       float *tit,float *wt,float *ampt,int *nz,float *dz,float *fz,
       float *tb,float *pb,int *nr,float *dr);    

void radix_(int *nt, int *nfft);
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);

void bilint_(int *n1, int *nx, int *ny, float *x0, float *y0, 
	float *dx, float *dy, float *x, float *y, float *f, float *g); 

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; 
	long long llimit;
	float tpow, zpow, sstart, lstart, t0;
	float xs, ys, xg, yg, dt, dzo, fzo, *fzout, v0, dvz; 
	float *hz, sss, lll, hzz, hzmin, hzmax;
	char *hzgrid;
	FILE *hzfp;
	usghed hzgh;
	float oshz,olhz,dshz,dlhz;
	int nshz,nlhz;
	int ihzgrid=0, itr;
	float *ss, *sl, *gs, *gl, ms, ml, ofomin, dofo, obtol; 
	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 strace, jthtrace;
	int mintrace, maxtrace;
	int minline, maxline;
	int sline, jthline;


	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");
		if(!getparint("strace",&strace)) strace = 0;
		if(!getparint("jthtrace",&jthtrace)) jthtrace = 1;
		if(!getparint("sline",&sline)) sline = 0;
		if(!getparint("jthline",&jthline)) jthline = 1;
		if(!getparint("mintrace",&mintrace)) mintrace = -2000000000;
		if(!getparint("maxtrace",&maxtrace)) maxtrace = 2000000000;
		if(!getparint("minline",&minline)) minline = -2000000000;
		if(!getparint("maxline",&maxline)) maxline = 2000000000;
	}

	/* 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);