velo2hfile.c

seismic software,very useful

#include "comva.h"
#include "par.h"

/*********************** self documentation *****************************/
string sdoc =
" 									\n"
" VELO2HFILE - VELO CARDS to HORIZON FILE (HFILE) CONVERSION 		\n"
" 									\n"
" velo2hfile <velocards oldhfile= [optional parameters] >newhfile	\n"
"\n"
"Required parameters:							\n"
"velocards              file name of VELO/DVDZ cards			\n"
"oldhfile=              file name of old hfile to update 		\n"
"newhfile=              file name of new hfile 				\n"
"Optional parameters: 							\n"
"vtopupdate=1,2,3,...   horizon numbers to update top velocities        \n"
"                       (default is all horizons)                       \n"
"vbotupdate=1,2,3,...   horizon numbers to update bottom velocities     \n"
"                       (default is all horizons)                       \n"
"ism=1                  smoothing option; to laterally smooth interval  \n"
"                         velocities at top and bottom, along interface,\n"
"                         before output  (applied to updated horizons   \n"
"                         (only)                                        \n"
"                         0=average all values at picks along interface and \n"
"                           this single value is assigned to all output \n"
"                           points along the interface \n"
"                         1=no smoothing \n"
"                         n=smooting values at picks along interface using \n"
"                           n data points \n"
"dcdp=from oldhfile     cdp spacing.  If not given in oldhfile, it mutst \n"
"                       be supplied	 				\n"
"xtol=0.                horizons extend, at both ends, a distance of xtol\n"
"                       to avoid discontinuity.               	\n"
"ztol=50.               vertical tolerance to choose VELO picks at horizon\n"
"                       positions 					\n"
"                       >0. --- will choose the pick within ztol from the \n"
"                               z position of the horizon at VELO location \n"
"                       =0. --- will interpolate velocity picks to small \n"
"                               z interval, then obtain velocity at horizon \n"
"                               depth					\n"
"\n"
" author: Zhiming Li		      		6/18/92			\n"
;
/**************** end self doc *******************************************/

/* define maximum number of horizons and maximum number of picks per horizon */
#define maxhs 128
#define maxpks 256

/* functions defined and used internally */

void vtopbot(int maxp,int maxcdp, 
     	int ncdp_velo,float *x_velo,float *z_velo,float *v_velo,int *nps_velo,
	int ncdp_dvdz,float *x_dvdz,float *zt_dvdz,float *zb_dvdz, 
	float *dvdz_dvdz,int *nps_dvdz,
	int nhs, float *xpicks, float *zpicks, int *npicks, float xtol, 
	float ztol, 
	float *veltops, float *velbots, float *velavgs, float *dvdzs,
	int *vtopupdate, int *vbotupdate, int *hnums);

void hfileout(FILE *hfilefp, float *xpicks, float *zpicks,
	float *veltop,float *velbot, int *dif,
	float *dens,float *qval, float *pois,
	int *npicks, int nhs, int ism, 
	float xmin, float xmax, float zmin, float zmax,
	int cdpxmin, int cdpxmax, float dcdp, 
	int *hnums, char hnames[maxhs][80],
	char *dunits, char *zattrib,
	int *vtopupdate, int *vbotupdate);

void findv(int ns, float *xs, float *vs, float x, float *v);  

void findbot(int maxp,int maxcdp,
	int nhs,float *xpicks,float *zpicks,int *npicks,float xtol,
	float x, int ihxz, int *ihbot);

void ilimit(int *i,int i1,int i2);

main (int argc, char **argv) 
{

	int *npicks,*hnums,ih,nhs,*difs,
		cdpxmini,cdpxmaxi,
		ism,ncdps_velo,ncdps_dvdz,maxp,maxcdp,
		*nps_velo,*nps_dvdz,*sortindex,*isort,
		*cdps_velo,*cdps_dvdz,jh,i,ihfile,
		vtopupdate[maxhs], vbotupdate[maxhs];
	float *xpicks,*zpicks,*veltops,*velbots,*velavgs,*dvdzs, 
		*dens,*qval,*pois,ztol,
		xmini,xmaxi,zmini,zmaxi,dcdpi,vmini,vmaxi,xtol,
		*z_velo,*x_velo,*v_velo,*sorts,*cdpsort,dcdp,
		*zt_dvdz,*zb_dvdz,*dvdz_dvdz,*x_dvdz;
	char hnames[maxhs][80], dunitsi[80], zattribi[80], vtypei[80],
		*oldhfile;
	FILE *infp=stdin, *hfoldfp, *outfp=stdout;

	/* initialize getpar */
	initargs(argc,argv);
	askdoc(1);

	/* get oldhfile name */
	if (!getparstring("oldhfile", &oldhfile)) 
		err("must specify oldhfile !");
	hfoldfp = efopen(oldhfile,"r"); 
	
	/* get horizon numbers to update velocities */
	for(ih=0;ih<maxhs;ih++) {
		vtopupdate[ih] = ih + 1;
		vbotupdate[ih] = ih + 1;
	}
	if (jh=getparint("vtopupdate",vtopupdate))
                for (ih=jh; ih<maxhs; ih++) vtopupdate[ih] = 0;

	if (jh=getparint("vbotupdate",vbotupdate))
                for (ih=jh; ih<maxhs; ih++) vbotupdate[ih] = 0;

	/* get length of smoothing operator */
	if(!getparint("ism",&ism)) ism=1;
	if(!getparfloat("xtol",&xtol)) xtol=0.;
	if(!getparfloat("ztol",&ztol)) ztol=50.;

	/* memory allocations */
	xpicks = (float *) malloc(maxhs*maxpks*sizeof(float));
	zpicks = (float *) malloc(maxhs*maxpks*sizeof(float));
	veltops = (float *) malloc(maxhs*maxpks*sizeof(float));
	velbots = (float *) malloc(maxhs*maxpks*sizeof(float));
	dens = (float *) malloc(maxhs*maxpks*sizeof(float));
	pois = (float *) malloc(maxhs*maxpks*sizeof(float));
	qval = (float *) malloc(maxhs*maxpks*sizeof(float));
	npicks = (int *) malloc(maxhs*sizeof(int));
	hnums = (int *) malloc(maxhs*sizeof(int));
	difs = (int *) malloc(maxhs*maxpks*sizeof(int));
	velavgs = (float *) malloc(maxhs*maxpks*sizeof(float));
	dvdzs = (float *) malloc(maxhs*maxpks*sizeof(float));

	for(ih=0;ih<maxhs;ih++) npicks[ih]=0;
	for(ih=0;ih<maxhs*maxpks;ih++) {
		veltops[ih]=0.;
		velbots[ih]=0.;
		difs[ih] = 0;
		velavgs[ih] = 0.;
		dvdzs[ih] = 0.;
	}

	/* read in old hfile*/
	nhs = 0;
	dcdpi = 0.;
	ifileread(hfoldfp,&nhs,xpicks,zpicks,veltops,velbots,difs,
		  dens,qval,pois,velavgs,dvdzs,
              	  (char *)hnames,hnums,npicks,maxhs,maxpks,
		  &xmini,&xmaxi,&zmini,&zmaxi,
		  &cdpxmini,&cdpxmaxi,&vmini,&vmaxi,
		  vtypei,dunitsi,zattribi,&dcdpi,&ihfile);
	if(!getparfloat("dcdp",&dcdp)) dcdp = 0.;
	if(dcdp!=0.) {
		if(dcdpi!=0. && dcdpi!=dcdp )
	warn("input dcdp not the same as dcdp from oldhfile, input dcdp used");
		dcdpi = dcdp;
	} 
	if(dcdpi==0.) err("dcdp must be specified \n");

	/* read in VELO cards and DVDZ cards dataset */
	ncdps_velo = 0;
	ncdps_dvdz = 0;
	maxp = 256;
	maxcdp= 128;

	cdps_velo = (int*) malloc(maxcdp*sizeof(int));
	x_velo = (float*) malloc(maxcdp*sizeof(float));
	z_velo = (float*) malloc(maxp*maxcdp*sizeof(float));
 	v_velo = (float*) malloc(maxp*maxcdp*sizeof(float));
	nps_velo = (int*) malloc(maxcdp*sizeof(int));
	sorts = (float*) malloc(maxp*maxcdp*sizeof(float));
	sortindex = (int*) malloc(maxcdp*sizeof(int));
	isort = (int*) malloc(maxcdp*sizeof(int));
	cdpsort = (float*) malloc(maxcdp*sizeof(float));

	veloread(infp,cdps_velo,z_velo,v_velo,
		  &ncdps_velo,nps_velo, maxp, maxcdp);
	/* sort cdps of VELO cards into ascending order */
	for(i=0;i<ncdps_velo;i++) {
		sortindex[i] = i;
		cdpsort[i] = cdps_velo[i];
	}
        if(ncdps_velo>0) qkisort(ncdps_velo,cdpsort,sortindex);
        for(i=0;i<ncdps_velo;i++) {
            	cdps_velo[i] = cdpsort[sortindex[i]];
                isort[i] = nps_velo[sortindex[i]];
        }
        for(i=0;i<ncdps_velo;i++) {
                nps_velo[i] = isort[i];
        }
        for(i=0;i<ncdps_velo;i++) {
                bcopy((char*)(v_velo+sortindex[i]*maxp),
                     (char*)(sorts+i*maxp),maxp*4);
        }
	for(i=0;i<ncdps_velo;i++) {
                bcopy((char*)(sorts+i*maxp),
                      (char*)(v_velo+i*maxp),maxp*4);
        }
        for(i=0;i<ncdps_velo;i++) {
                bcopy((char*)(z_velo+sortindex[i]*maxp),
                      (char*)(sorts+i*maxp),maxp*4);
        }
	for(i=0;i<ncdps_velo;i++) {
                bcopy((char*)(sorts+i*maxp),
                      (char*)(z_velo+i*maxp),maxp*4);
        }
	/* DVDZ cards */
	cdps_dvdz = (int*) malloc(maxcdp*sizeof(int));
	x_dvdz = (float*) malloc(maxcdp*sizeof(float));
	zt_dvdz = (float*) malloc(maxp*maxcdp*sizeof(float));
	zb_dvdz = (float*) malloc(maxp*maxcdp*sizeof(float));
	dvdz_dvdz= (float*) malloc(maxp*maxcdp*sizeof(float));
        nps_dvdz = (int*) malloc(maxcdp*sizeof(int));
	dvdzread(infp,cdps_dvdz,zt_dvdz,zb_dvdz,dvdz_dvdz,
		 &ncdps_dvdz,nps_dvdz,maxp,maxcdp);
	/* sort cdps of DVDZ cards into ascending order */
	for(i=0;i<ncdps_dvdz;i++) {
		sortindex[i] = i;
		cdpsort[i] = cdps_dvdz[i];
	}
    	if(ncdps_dvdz>0) qkisort(ncdps_dvdz,cdpsort,sortindex);
    	for(i=0;i<ncdps_dvdz;i++) {
        	cdps_dvdz[i] = cdpsort[sortindex[i]];
        	isort[i] = nps_dvdz[sortindex[i]];
    	}
    	for(i=0;i<ncdps_dvdz;i++) {
        	nps_dvdz[i] = isort[i];
	}
    	for(i=0;i<ncdps_dvdz;i++) {
        	bcopy((char*)(dvdz_dvdz+sortindex[i]*maxp),
		      (char*)(sorts+i*maxp),maxp*4);
	}
    	for(i=0;i<ncdps_dvdz;i++) {
       		bcopy((char*)(sorts+i*maxp),
		      (char*)(dvdz_dvdz+i*maxp),maxp*4);
	}
    	for(i=0;i<ncdps_dvdz;i++) {
        	bcopy((char*)(zt_dvdz+sortindex[i]*maxp),
		      (char*)(sorts+i*maxp),maxp*4);
	}
    	for(i=0;i<ncdps_dvdz;i++) {
        	bcopy((char*)(sorts+i*maxp),
		      (char*)(zt_dvdz+i*maxp),maxp*4);
	}
    	for(i=0;i<ncdps_dvdz;i++) {
       		bcopy((char*)(zb_dvdz+sortindex[i]*maxp),
		      (char*)(sorts+i*maxp),maxp*4);
	}
    	for(i=0;i<ncdps_dvdz;i++) {
       		bcopy((char*)(sorts+i*maxp),
		      (char*)(zb_dvdz+i*maxp),maxp*4);
	}
	free(isort);
	free(sorts);
	free(sortindex);
	free(cdpsort);

	/* convert cdp locations of VELO into x locations */
	for(i=0;i<ncdps_velo;i++) { 
		x_velo[i] = xmini + (cdps_velo[i]-cdpxmini)*dcdpi;
	}
	/* convert cdp locations of DVDZ into x locations */
	/* still store in cdps_dvdz */
	for(i=0;i<ncdps_dvdz;i++) 
		x_dvdz[i] = xmini + (cdps_dvdz[i]-cdpxmini)*dcdpi;
	free(cdps_velo);
	free(cdps_dvdz);
	

        vtopbot(maxp,maxcdp,ncdps_velo,x_velo,z_velo,
		v_velo,nps_velo,
		ncdps_dvdz,x_dvdz,zt_dvdz,zb_dvdz, 
		dvdz_dvdz,nps_dvdz,
		nhs,xpicks,zpicks,npicks,xtol,ztol, 
		veltops,velbots,velavgs,dvdzs,vtopupdate,vbotupdate,
		hnums);

	hfileout(outfp, xpicks, zpicks,
		veltops, velbots, difs,
		dens,qval,pois,
		npicks, nhs, ism, 
		xmini, xmaxi, zmini, zmaxi,
		cdpxmini, cdpxmaxi, dcdpi, 
		hnums, hnames,
		dunitsi, zattribi,
		vtopupdate,vbotupdate);

	return EXIT_SUCCESS;
}

void ilimit(int *i, int i1, int i2) {
        if (*i<i1) {
                *i = i1;
        }else if(*i>i2) {
                *i = i2;
        }
}


void hfileout(FILE *hfilefp, float *xpicks, float *zpicks,
	float *veltop,float *velbot, int *dif,
	float *dens,float *qval, float *pois,
	int *npicks, int nhs, int ism, 
	float xmin, float xmax, float zmin, float zmax,
	int cdpxmin, int cdpxmax, float dcdp, 
	int *hnums, char hnames[maxhs][80],
	char *dunits, char *zattrib,
     	int *vtopupdate, int *vbotupdate) {


	int ic,ip,np,ih,jp,nhout,h0; 
	int p1,p2,i,ismx,ismxh;
	float vmax, vmin, s;
	float *w, *f;
	float vt, vb, g;

	int jh;



	/* lateral smooth velocities */
	if(ism==0) {
	/* take average */ 
		for(ih=0;ih<nhs;ih++) {
			vt = 0.;
			vb = 0.;
			g = 0.;
			np = npicks[ih];
			for(jh=0;jh<maxhs;jh++) {
				if(ih+1 == vtopupdate[jh]) {
					for(ip=0;ip<np;ip++) {
						vt += veltop[ip+ih*maxpks];
					}
					if(np>0) vt = vt/np; 
					for(ip=0;ip<np;ip++) 
						veltop[ip+ih*maxpks] = vt;
					break;
				}
			}
			for(jh=0;jh<maxhs;jh++) {
				if(ih+1 == vbotupdate[jh]) {
					for(ip=0;ip<np;ip++) {
						vb += velbot[ip+ih*maxpks];
					}
					if(np>0) vb = vb/np; 
					for(ip=0;ip<np;ip++) 
						velbot[ip+ih*maxpks] = vb;
					break;
				}
			}

		}
	} else if (ism>1) {
	/* smoothing */
		for(ih=0;ih<nhs;ih++) {
			h0 = ih * maxpks;
			np = npicks[ih];
			if(np<2) continue;
			ismx = ism;
			ismx = (ismx/2)*2+1;
			if(ismx>np) ismx=(np/2)*2+1;
			f = (float *) malloc(ismx*sizeof(float));
			w = (float *) malloc(np*sizeof(float));
			/* triangular smoothing coefficients */
			ismxh = (ismx+1)/2;
			s = 0.;
			for(i=0;i<ismx;i++) {
				f[i] = ismxh - abs(i+1-ismxh);
				s += f[i];
			}
			s = 1./s;
			for(i=0;i<ismx;i++) {
				f[i]=f[i]*s;
			}
			/* apply smoothing to veltop and velbot */
			ic = 1;
			for(jh=0;jh<maxhs;jh++) {
				if(ih+1 == vtopupdate[jh]) {
					smth1d_(veltop+h0,f,w,&np,&ismx,&ic);
					for(ip=0;ip<np;ip++) 
						veltop[ip+h0] = w[ip];
					break;
				}
			}
			for(jh=0;jh<maxhs;jh++) {
				if(ih+1 == vbotupdate[jh]) {
					smth1d_(velbot+h0,f,w,&np,&ismx,&ic);
					for(ip=0;ip<np;ip++) 
						velbot[ip+h0] = w[ip];
					break;
				}
			}

			free(w);
			free(f);
		}
	}

	/* find maximum and minimum velocities */
	vmin = veltop[0];
	vmax = veltop[0];
	for (ih=0;ih<nhs;ih++) {
		for(ip=0;ip<npicks[ih];ip++) {
			veltop[jp] += 0.5;
			velbot[jp] += 0.5;
			jp = ip + ih * maxpks;
			if(vmin>veltop[jp]) vmin = veltop[jp];
			if(vmin>velbot[jp]) vmin = velbot[jp];
			if(vmax<veltop[jp]) vmax = veltop[jp];
			if(vmax<velbot[jp]) vmax = velbot[jp];
		}
	}

	nhout = 0;
	for(ih=0;ih<nhs;ih++) {
		if(npicks[ih]>1) nhout = nhout + 1;
	}


	/* output picks */

	fprintf(hfilefp, 
	"HFILE   <--- FILE FORMAT DESCRIPTOR \n");
	fprintf(hfilefp,"\n");
	fprintf(hfilefp, "MODEL.X-MIN         %8.1f\n",xmin);
	fprintf(hfilefp, "MODEL.X-MAX         %8.1f\n",xmax);
	fprintf(hfilefp, "MODEL.Z-MIN         %8.1f \n",zmin);
	fprintf(hfilefp, "MODEL.Z-MAX         %8.1f \n",zmax);
	fprintf(hfilefp, "MODEL.CDP#-AT-X-MIN %6d\n",cdpxmin);
	fprintf(hfilefp, "MODEL.CDP#-AT-X-MAX %6d\n",cdpxmax);
	fprintf(hfilefp, "MODEL.MIN-VELOCITY  %8.1f\n",vmin); 
	fprintf(hfilefp, "MODEL.MAX-VELOCITY  %8.1f\n",vmax); 
	fprintf(hfilefp, "MODEL.#HORIZONS     %6d\n",nhout); 
	fprintf(hfilefp, "MODEL.VELOCITY-TYPE  INTERVAL\n"); 
	fprintf(hfilefp, "MODEL.DIST-UNITS     %s\n",dunits); 
	fprintf(hfilefp, "MODEL.Z-ATTRIBUTE    %s\n",zattrib); 
	fprintf(hfilefp, "MODEL.CDP-SPACING   %8.1f\n",dcdp); 
	fprintf(hfilefp, "\n"); 

	
	for(ih=0;ih<nhs;ih++) {

		np = npicks[ih];
		if(np<=1) continue;

		fprintf(hfilefp, "\n"); 
		fprintf(hfilefp, "\n"); 
		fprintf(hfilefp, "HORIZON-NUMBER %3d\n", hnums[ih]); 
		fprintf(hfilefp, "HORIZON-NAME     %s\n", hnames[ih]); 
		fprintf(hfilefp, "* \n"); 
		fprintf(hfilefp, 
"* X-VALUE   Z-VALUE   VELTOP   VELBOT  DIF  DENS  QVAL  POIS \n");
		fprintf(hfilefp, 
"* =======  ========  =======  =======  ===  ====  ====  ==== \n");

		h0 = ih*maxpks;
		
		for (ip=0;ip<npicks[ih];ip++) {