kmmerge.c

seismic software,very useful

#include "subc.h"
#include "ghdr.h"
#include "gridhd.h"
#include "grid.h"
#include "usgrid.h"
#include "comva.h"
#include "usu.h"
#include "su.h"
#include "segy.h"
#include "header.h"
#include "cwp.h"

char *sdoc[] = {
"KMMERGE - Merge results from different-input/same-output KZMIG/KTMIG runs",
" ",
"kmmerge <input [parameters] >output 				", 
" ",
"Required parameters:							",
"input=                 any output of one KZMIG/KTMIG run		",
"output=                merged seismic output data set name 		",
"x1=                    x coordinate of 1st corner of the 3D master grid",
"y1=                    y coordinate of 1st corner of the 3D master grid",
"s1=                    s coordinate of 1st corner of the 3D master grid",
"l1=                    l coordinate of 1st corner of the 3D master grid",
"cdp1=                  cdp number of 1st corner of the 3D master grid	",
"x2=                    x coordinate of 2nd corner of the 3D master grid",
"y2=                    y coordinate of 2nd corner of the 3D master grid",
"s2=                    s coordinate of 2nd corner of the 3D master grid",
"l2=                    l coordinate of 2nd corner of the 3D master grid",
"cdp2=                  cdp number of 2nd corner of the 3D master grid\n",
"x3=                    x coordinate of 3rd corner of the 3D master grid",
"y3=                    y coordinate of 3rd corner of the 3D master grid",
"s3=                    s coordinate of 3rd corner of the 3D master grid",
"l3=                    l coordinate of 3rd corner of the 3D master grid",
"cdp3=                  cdp number of 3rd corner of the 3D master grid	",
"dds=                   cdp spacing in inline direction                 ",
"ddl=                   Line spacing in crossline                       ",
" ",
"diskimg=               Disk name for storing imaged time sections 	",
"diskhdr=               Disk name for storing output trace headers  	",
"diskfld=               Disk name for storing output trace folds  	",
"                       (diskimg,diskhdr,diskfld must be repeated 	",
"                        n times for n inputs)				",
"kzorkt=                kzmig or ktmig merge (0=kzmig 1=ktmig)		",
"Optional parameters:			",
"when kzorkt=0 the following 5 parameters are used			",
"nzo=from_input         number of output samples per trace		",
"fzo=from_input         depth of first sample of output trace		",
"dzo=from_input         depth sampling interval of ouptut trace		",
"ttfile=ttfile          file name of input traveltime tables            ",
"or fzt= dzt= (when ttfile is missing) \n"
"fzt=                   depth of first sample of travel time ",
"dzt=                   depth sampling interval of travel time ",
"zpow=0.0               power of depth gain to be appled to the output  ",
"                       after migration   ",
"when ktorkt=1 the following 4 parameters are used			",
"ntau=from_input        number of output samples per trace		",
"tau0=from_input        time of first sample of output trace (in ms)	",
"dtau=from_input        time sampling interval of ouptut trace (in ms)	",
"tpowaf=0.0             power of time gain to be appled to the output    ",
"                       after migration ",
" ",
"flexbin=0              1=Input data to migration have been flexbined	",  
"                         the fold dataset will not be applied to output",
"                       0=sum the folds and apply to the summed image	",
"oofo=0.                output minimum offset (>0.)                     ",
"dofo=999999.           output offset increment                         ",
"nofo=1                 number of offsets			",
"ofo=0.                 output offsets  (specified as 200,250.,350,600,...) ",
"                       if specified, dofo and oofo will be ignored     ",
"sstart=s1              Minimum s (trace position within line) of       ",
"                       migration output                                ",
"                       (NOTE: this is coordinates, i.e., first trace   ",
"                       of a line at s1, 2nd at s1+ds, 3rd at s1+2*ds, ...)",
"ns=(s2-s1)/ds+1        Number of cdp (group) per line of migration output ",
"lstart=l1              Minimum l (line position) of migration output   ",
"                       (NOTE: this is coordinates, i.e., first line    ",
"                       l1, 2nd at l1+dl, 3rd at l1+2*dl, ...)          ",
"nl=(l3-l1)/dl+1        Number of lines (group) of migration output     ",
"ds=dds                 cdp group interval to output                    ",
"dl=ddl                 Line group  interval to output                  ",
"send=                  Ending s position of migration output           ",
"lend=                  Ending l position of migration output           ",
"ntrend=                Number of cdp bins to output from (sstart,lstart) ",
"                       to (send,lend)                                  ",
"                       (specify send, lend and ntrend will output only ",
"                       one slant line from (sstart,lstart) to          ",
"                       (send,lend) of total ntrend cdps; if not        ",
"                       specified, will output to ns by nl cdps)        ",
"nds=1                  Number of cdp locations per output cdp group    ",
"ndl=1                  Number of line locations per output line group  ",
"hzgrid=                horizon grid name to specify the starting depth ",
"                       of output depth migraiton. ",
"                       if given, delrt in the trace header will be 	",
"                       kept. otherwise reset to tau0 or fzo			",
"tracekey=              key word for trace number ",
"linekey=               key word for line number ",
"	",
" ",
"AUTHOR:		Z Li,        12/94   			",
NULL};

segytrace tra, tr;
segybhdr bh;
segychdr ch;


main(int argc, char **argv)
{
	float x1,x2,x3,y1,y2,y3,s1,s2,s3,l1,l2,l3;
	float dds,ddl,ds,dl;
	float send, lend;
	int cdp1, cdp2, cdp3;
	int ntrend,nsout,nlout,ntmp,il,is,flexbin;
	float *s, *l, *migs, *mig, *fold, *folds; 
	float zi, scale;
	string *diskimg, *diskfld, *diskhdr;
	char *ttfile;
	FILE **difp, **dffp, **dhfp, *ttfp, *infp=stdin, *outfp=stdout;
	int ni, i, lpos,izi,iz,it;
	float mx,my,dd,ofs;
	float ms, ml;
	int nofo, ns, nl, nzo;
	int nz, ntau, kzorkt,nz0,ierr,io;
	float fz,fzo,dzo,tau0,dtau,fzt,dzt;
	float ofomin, dofo, sstart,lstart,tmp,fold0; 
	int ndl,nds,ncdppl,nlines,nlread,nsread;
	int cdppds, cdppdds, cdppdl, cdppddl, cdpnum;
        int iix, iiy;
	float *zgain, zpow, tpowaf;
	String tracekey, linekey;
	int ikey=0;

	float *ofo;
	usghed ugh;
	
	int itmp,itmp1;

	char *hzgrid;
	int ihz=1;
	int noinput=1;
	long long isize;

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



   	fseek2g(outfp,0,1);
	fseek2g(infp,0,1);
	fseek2g(infp,0,SEEK_END);
	isize= ftell64(infp);

	if(isize>3600) {
		fseek2g(infp,0,0);
		/* read id headers */
    	fgethdr(infp,&ch,&bh);
		fgettr(infp,&tr);
		noinput = 0;
	}


	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)) err("Must specify ddl!\n"); 

		if(!getparstring("hzgrid",&hzgrid)) ihz=0;


	ni = countparname("diskimg");
	if(ni<=1) err(" number of inputs %d ust be greater than 1",ni);
	if(countparname("diskfld")!=ni) err("check diskfld");
	if(countparname("diskhdr")!=ni) err("check diskhdr");

	diskimg = (string *) emalloc(ni*sizeof(string)); 
	diskfld = (string *) emalloc(ni*sizeof(string)); 
	diskhdr = (string *) emalloc(ni*sizeof(string)); 
	difp = (FILE **)emalloc(ni*sizeof(FILE *));
	dffp = (FILE **)emalloc(ni*sizeof(FILE *));
	dhfp = (FILE **)emalloc(ni*sizeof(FILE *));

	for(i=0;i<ni;i++) {
		if(!getnparstring(i+1,"diskimg",&diskimg[i]))
                        err("cannot get diskimg file for %i -th input \n",i+1);
		if(!getnparstring(i+1,"diskfld",&diskfld[i]))
                        err("cannot get diskfld file for %i -th input \n",i+1);
		if(!getnparstring(i+1,"diskhdr",&diskhdr[i]))
                        err("cannot get diskhdr file for %i -th input \n",i+1);

		difp[i]=efopen(diskimg[i],"r");
		dffp[i]=efopen(diskfld[i],"r");
		dhfp[i]=efopen(diskhdr[i],"r");
	}

	fprintf(stderr," %d sets of disk files open ... \n", ni);

	if(!getparint("kzorkt",&kzorkt)) err(" kzorkt missing");
	if(kzorkt==0) {
		if(noinput==0) {
			if(!getparint("nzo",&nzo)) nzo = tr.ns;
			if(!getparfloat("dzo",&dzo)) dzo = tr.dz;
			if(!getparfloat("fzo",&fzo)) fzo = tr.fz;
		} else {
			if(!getparint("nzo",&nzo)) err(" nzo missing");
			if(!getparfloat("dzo",&dzo)) err(" dzo missing");
			if(!getparfloat("fzo",&fzo)) err(" fzo missing");
		}
		if(!getparfloat("zpow",&zpow)) zpow = 0.0;
		if(!getparstring("ttfile",&ttfile))  err(" ttfile missing");
		if((ttfp = fopen(ttfile,"r"))!=NULL) {
			fseek2g(ttfp,0,1);
			ierr = fgetusghdr(ttfp,&ugh);
        	if (ierr!=0) err(" Grid parameters of %s required!\n",ttfile);
			fzt = ugh.o1;   dzt = ugh.d1;
		} else {
			if(!getparfloat("fzt",&fzt)) err(" missing fzt ");
			if(!getparfloat("dzt",&dzt)) err(" missing dzt ");
		}
		nz0 = (int)((fzo-fzt)/dzt);
		fz = fzt+nz0*dzt;
        nz = 1+(int)((fzo+(nzo-1)*dzo-fzt)/dzt+0.99)-nz0;

	} else {
		if(noinput==0) {
			if(!getparint("ntau",&ntau)) ntau = tr.ns;
        	if(!getparfloat("dtau",&dtau)) {
				dtau = tr.dt * 0.000001;
			} else {
				dtau = dtau * 0.001;
			}
        	if(!getparfloat("tau0",&tau0)) {
				tau0 = tr.delrt * 0.001;
			} else {
				tau0 = tau0 * 0.001;
			}
		} else {
			if(!getparint("ntau",&ntau)) err(" ntau missing");
        	if(!getparfloat("dtau",&dtau)) err(" dtau missing");
			dtau = dtau * 0.001;
        	if(!getparfloat("tau0",&tau0)) err(" tau0 missing");
			tau0 = tau0 * 0.001;
		}
		if(!getparfloat("tpowaf",&tpowaf)) tpowaf = 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);
        } else {
                for(io=0;io<nofo;io++) ofo[io] = ofomin + io*dofo;
        }
	
	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(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) {
                fold0 = 10.;
        }
        else {
                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; }