ktmig.c

seismic software,very useful

void radix_(int *nt,int *nfft);

void bisear_(int *n,int *nnew,float *x,float *xnew,int *indx);

void f2p5n_(float *data,int *nt,int *i3d,int *m, int *l);

void vr22vi2_(float *vr2,float *vi2,int *ntau);
void vi22vq4_(float *vi2,float *vq4,int *ntau);
void vt2s2_(float *vr2,float *vq4,float *s2,float *tau,int *ntau);

void pstm3d_(float *trace,int *nt,float *t0,float *dt,float *sx,float *sy,
	float *gx,float *gy,int *imute,
	float *fovt2,float *mig,int *nxy,int *ntm,float *x,float *y,
	float *tm,float *aperx,float *apery,float *fold,float *w1,
	float *f0,float *df,int *nf,float *ftaper,
	int *ifcut,int *ltaper,float *tracef,
	float *wsave,int *nsave,float *work1,float *work2,int *nfft,
	float *kxmax,float *kymax,float *dtau,float *tau0,float *angmax,
	int *indxw,int *nindxw,int *incdx,int *np,
	float *resamp,int *ires,int *ntres,float *s2,float *scs,float *scg);


segytrace tra, tr;
segybhdr bh;
segychdr ch;


main(int argc, char **argv)
{
	int nlines, ncdppl;
   	int nt,ntau,it,ix,iy;
    	FILE *infp,*outfp,*velfp,*imgfp,*hdrfp,*hffp,*fldfp;
	float x1,x2,x3,y1,y2,y3,s1,s2,s3,l1,l2,l3,sm,lm;
	char *backupi="", *backupo="";
	char *diskimg="DISKIMAGE", *diskhdr="DISKHEADER";
	char *velfile,*datain,*dataout;
	char *diskfld="DISKFOLD", *hisfile="KTMIG.HISFILE";
	int ibacki, ibacko,isave, ihis=0, ipre=0;
	float ds, dl, tmp;
	int mtrace, ntrpre=0, ktrace, itrace, jtrace, ntrace, nofo;
	int ltrace; 
	int *imutel, cdppre, cdpnow, m1000, ifact, it0;
	int vread=0, ivdisk=1;
	float *velos;
	float *mig, *fold, *tm, *tmig, *tgain, *s, *l, *tugain, *trace, *tras;
	float *fovt2;
	float *ofo;
	int *off, itgain, nl, ns, ntras, ntv, is, il, io, nsv, nlv;
	int nsread, nlread;
	int mlimit, nlcore, nlpos, nlmig, *iofs, itv; 
	float tpow, sstart, lstart, t0, dtv, t0v, dsv, dlv, s0v, l0v;
	float tpowaf;
	float xs, ys, xg, yg, dt, dtau, tau0, v, *vel; 
	float *ss, *sl, *gs, *gl, ms, ml, ofomin, dofo, obtol; 
	char *buf;
	int ilv, isv, itmp, itmp2, ntl, itr, ir0, i;
	float tminl, dtl, dldt, apers, aperl;
	float *w1;
	float f0, df, ftaper;
	int nf, nfft, *ifcut, *ltaper, nsave;
	int mm,ll;
	float *work1, *work2, *tracef, *wsave;
	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, ilend;
	float send, lend, dsout, dlout, dend, *disend;
	float ps, pl;
	int iofoin=0;
	float ofol, ofor;
	int one=1; 
	float *resamp;
	int ires, ntres; 
	int *indxw, nindxw, incdxw;
	int ntrdsk, ndskwt;

	float *vr2, *vi2, *vq4, *tau, *ss2, *scs, *scg;
	int i3d, ncpu, traceout;
	char *envs;
	int on2trace, ntotal;
	
	
	ghed gh;
	int n1,n2,n3,n4,n5,dtype,ierr,orient,gtype;
	float o1,o2,o3,o4,o5,d1,d2,d3,d4,d5,ocdp2,oline3,dcdp2,dline3;
	float gmin=0.,gmax,scale;
	
	long long lwork, ltmp;
	long long lpos;

	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 ispre, isnow, ilpre, ilnow;

   	/* get parameters */
   	initargs(argc,argv);
   	askdoc(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.;
    	if (!getparstring("velfile",&velfile)) err("must specify velfile");

	/* 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 (!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("tpowaf",&tpowaf)) tpowaf = 0.;
	itgain = 0;
	if(tpow!=0. || tpowaf!=0.) itgain = 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;

	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 */
    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)tra.delrt/1000;

	ntotal = tra.tracl;

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

	/* get velociry grid header info */
	velfp = efopen(velfile,"r");
	fseek2g(velfp,0,1);
	ierr = fgetghdr(velfp,&gh);
	if (ierr==0) fromghdr(&gh,&scale,&dtype,&n1,&n2,&n3,
        		&n4,&n5,&d1,&d2,&d3,&d4,&d5,&o1,&o2,&o3,&o4,&o5,
        		&dcdp2,&dline3,&ocdp2,&oline3,
        		&gmin,&gmax,&orient,&gtype);

	if(gmin==0. && ierr==0) 
		warn(" gmin=0 from velocity header file");

	if(gmin==0.) {
    		if (!getparfloat("vmin",&gmin)) err(" vmin must be specified");
	}

	/* optional parameters */
    	if (!getparint("ntv",&ntv)) {
		if(ierr==0) { ntv = n1; } else { err("Must specify ntv "); }
	}
    	if (!getparfloat("dtv",&dtv)) {
		if(ierr==0) { dtv = d1; } else { err("Must specify dtv "); }
	}
    	if (!getparfloat("t0v",&t0v)) { 
		if(ierr==0) { t0v = o1; } else { err("Must specify t0v "); }
	}
    	if (!getparint("nsv",&nsv)) {
		if(ierr==0) { nsv = n2; } else { err("Must specify nsv "); }
	}
    	if (!getparfloat("dsv",&dsv)) {
		if(ierr==0) { dsv = d2; } else { err("Must specify dsv "); }
	}
    	if (!getparfloat("s0v",&s0v)) { 
		if(ierr==0) { s0v = o2; } else { err("Must specify s0v "); }
	}
    	if (!getparint("nlv",&nlv)) {
		if(ierr==0) { nlv = n3; } else { err("Must specify nlv "); }
	}
    	if (!getparfloat("dlv",&dlv)) {
		if(ierr==0) { dlv = d3; } else { err("Must specify dlv "); }
	}
    	if (!getparfloat("l0v",&l0v)) { 
		if(ierr==0) { l0v = o3; } else { err("Must specify l0v "); }
	}

	if(dlv==0.)  dlv = dl;
	if(dsv==0.)  dsv = ds;
	 

   	if (!getparint("ntau",&ntau)) ntau = nt;
   	if (!getparfloat("dtau",&dtau)) {