ktmig.c

seismic software,very useful

	
	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;

   	/* 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 (!getparint("etatype",&etatype)) etatype=0;
    	if (!getparstring("etafile",&etafile)) etatype=-1;

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

	/* 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(etatype!=-1) { 
		etafp = efopen(etafile,"r");
		fseek2g(etafp,0,1);
	}

	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)) {
		dtau = dt;
	} else {
		dtau = dtau * 0.001;
	}
   	if (!getparfloat("tau0",&tau0)) {
		tau0 = 0.;
	} else {
		tau0 = tau0 * 0.001;
	}
		
   	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 = 7;
    	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 = ntau;
	bh.hdt = dtau * 1000000; 
	if(traceout==1) fputhdr(outfp,&ch,&bh);


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

	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;

	ibacki = 0;
	if (getparstring("backupi",&backupi)) {
		ibacki = 1;
		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;

	fprintf(jpfp,"\n");
	fprintf(jpfp," -------- KTMIG 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 memory requirement for migration*/
 
	ntmp = nl*ndl;
	nl = ntmp;
	ntmp = ns*nds;
	ns = ntmp;

	nlcore = nl;

	incdxw = 10;
	nindxw = ntau / incdxw + 1;
	ires = 4;
	ntres = nt*ires + 1;
	

	if(ntrend==0) {
		nsout = ns;
		nlout = nl;
	} 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);
	}

	lwork = nsout * nlcore;
	lwork = lwork*ntau*nofo*sizeof(float);
	ltmp = nsout*nlout*nofo+nt + ntau + ntl + nt*mtrace;
	ltmp = ltmp + 2*ntau + 2*ntau*mtrace + 2*nsout*nlout + ntv + mtrace*4;
	ltmp = ltmp + ntau + nfft*(nf+nsout*nlout) + ntl*nf + nsave;
	ltmp = ltmp + nsout*nofo + ntres*nf + ntau;
	ltmp = ltmp + ntau*nsout*nlout*2 + ntau*mtrace + ntau*3;
	ltmp = ltmp * sizeof(float);
	ltmp = ltmp+(nsout*nlout*nofo+mtrace*2+nf*2+nindxw*nsout*nlout)*sizeof(int);
	lwork = lwork + ltmp;
	ltmp = ntv*nsv*nlv*sizeof(float);
	lwork = lwork + ltmp;

	ivdisk = 0;
	if(lwork>mlimit) {
		warn("Velocity Disk I/O Used to Reduce Memory \n");
		lwork = lwork - ltmp;
		ivdisk = 1;
	}

	if(lwork>mlimit) err("mlimit too small ");
	
	fprintf(jpfp," Total Memory used (in Byte) = %g \n", (float)lwork);
	fprintf(jpfp,"\n");

	/* allocation of memory */
	mig = (float*) emalloc(nsout*nlcore*nofo*ntau*sizeof(float));	
	fold = (float*) emalloc(nsout*nlout*nofo*sizeof(float));	
	off = (int*) emalloc(nsout*nlout*nofo*sizeof(int));	
	tgain = (float*) emalloc(nt*sizeof(float));
	tugain = (float*) emalloc(ntau*sizeof(float));
	trace = (float*) emalloc(ntl*sizeof(float));	
	tras = (float*) emalloc(nt*mtrace*sizeof(float));	
	imutel = (int*) emalloc(mtrace*sizeof(int));
	tm = (float*) emalloc(ntau*sizeof(float));
	tmig = (float*) emalloc(ntau*sizeof(float));
	fovt2 = (float*) emalloc(ntau*mtrace*sizeof(float));
	if(etatype!=-1) {
		etamig = (float*) emalloc(ntau*mtrace*sizeof(float));