qcvs3d.c

seismic software,very useful

/* quality control of velocity 3d cards */

#include "velo.h"
#include "usgrid.h"
#include "ghdr.h"
#include "gridhd.h"
#include "comva.h"
#include "par.h"



char *sdoc = 
"QCVS3D - quality control display of VS3D cards and volocity grids \n"
"\n"
"qcvs3d [parameters] xqc= yqc=					\n" 
"\n"
"Required parameters:						 	\n"
"xqc=              inline (trace) coordinate of VS3D to qc 			\n"
"yqc=              crossline (line) coordinate of VS3D to qc 		\n"
"\n"
"Optional parameters:							\n"
"ixy=0             x and y coordinate system code \n"
"                  =0 xqc and yqc are inline and crossline distance \n"
"                  =1 xqc is trace number and yqc is line number \n"
"vs3d1=            Name of first VS3D dataset to qc 			\n"
"vs3d2=            Name of second VS3D dataset to qc			\n"
"vs3d3=            Name of third VS3D dataset to qc			\n"
"dt=100            time interval in ms to interpolate input t-v pairs   \n"
"                  before converting to interval velocity 		\n"
"                  =0 will compute interval velocity at input times	\n"
"vstype1=0         Velocity type of VS3D card (0=rms 1=avg 2=interval)	\n"
"vstype2=0         Velocity type of VS3D card (0=rms 1=avg 2=interval)	\n"
"vstype3=0         Velocity type of VS3D card (0=rms 1=avg 2=interval)	\n"
"vgrid1=           Name of 1st 3D interval velocity grid stored as v(t,x,y) \n"
"                  i.e., stored as velocity vectors not time slices 	\n"
"vgtype1=2         Velocity type of 1st vgrid (0=rms 1=average 2=interval)\n"
"vgrid2=           Name of 2nd 3D interval velocity grid stored as v(t,x,y) \n"
"                  i.e., stored as velocity vectors not time slices 	\n"
"vgtype2=2         Velocity type of 2nd vgrid (0=rms 1=average 2=interval)\n"
"vgrid3=           Name of 3rd 3D interval velocity grid stored as v(t,x,y) \n"
"                  i.e., stored as velocity vectors not time slices 	\n"
"vgtype3=2         Velocity type of 3rd vgrid (0=rms 1=average 2=interval) \n"
"paper=0           Paper plot (1=yes 0=screen)				\n"
"vauxplot=0        Plot other velocity curve (0=no 1=rms 2=avg)		\n"
"viplot=1          Plot interval velocity curves (0=no 1=yes)	\n"
"nvfmax=4096       maximum number of velocity functions each vs3d data	\n"
"                  set may have						\n"
"ntvmax=256        maximum number of t-v pairs per velocity function in \n"
"                  vs3d datasets					\n"
"tmax=             maximum time to disply (default to data)		\n"
"title=            title of the plot							\n"
"Note:									\n"
" 1. see vgrid3d for the format of VS3D					\n"
" 2. velocity curve from vs3d1 will be ploted with red color	\n"
"    velocity curve from vs3d2 will be ploted with green color	\n"
"    velocity curve from vs3d3 will be ploted with dark blue color\n"
"    velocity curve from vgrid1 will be ploted with light blue color\n"
"    velocity curve from vgrid2 will be ploted with violet color\n"
"    velocity curve from vgrid3 will be ploted with yellow color\n"
" 3. by default, inline and crossline coordinates of the velocity \n"
"    grid are defined in o2,o3,d2,d3,n2,n3 of the velocity gridheader. \n"
"    however, when ixy=1, inline and crossline coordinates of the velocity \n"
"    grid are defined in ocdp2,dcdp2,oline3,dline3 of the velocity \n"
"    gridheader \n"
"\n"
"AUTHOR:	   Zhiming Li,       ,			9/14/93		\n"    
;

void times(float dt, float *t, int ntv, float *time, int *nt, float tmax);

void vrms2vint(float *t, float *v, int ntv, float *time, float *vi, int nt);

void vavg2vint(float *t, float *v, int ntv, float *time, float *vi, int nt);

void vint2vrms(float *t, float *vi, float *vrms, int nt);

void vint2vavg(float *t, float *vi, float *va, int nt);

void vrmsout(FILE *datafp, float *ts, float *vs, float *work, int ns);

void vavgout(FILE *datafp, float *ts, float *vs, float *work, int ns);

void findvf(char *in, float xqc, float yqc, float *t, float *v, int *np, 
	float *xx, float *yy, int ntvmax, int nvfmax); 

main(int argc, char **argv)
{
   	FILE *vfp1, *vfp2, *vfp3;
	FILE *datafp;

	char *in1, *in2, *in3, *vgrid1, *vgrid2, *vgrid3;
	char cname[80];
	int i1, i2, i3, iv1, iv2, iv3, it;
	float xqc, yqc, dt; 
	int vstype1, vstype2, vstype3, vgtype1, vgtype2, vgtype3, nt;
	int paper, vauxplot, ia, ii, viplot;
	int ixy;

   	float *time, *work;
	float *t1s, *t2s, *t3s, *t4s, *t5s, *t6s;
	float *v1s, *v2s, *v3s, *v4s, *v5s, *v6s;
	int n1s, n2s, n3s, n4s, n5s, n6s;

	float *t1, *v1, *t2, *v2, *t3, *v3, *t4, *v4, *t5, *v5, *t6, *v6; 
	float x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6;
	int np1=0, np2=0, np3=0, np4=0, np5=0, np6=0;
	int ntmax=4096;
	float tmp, tmax;
	int ix, iy, iseek;
	char cmd[1024];
	int ntvmax,nvfmax;
	char *title;
	

	float x4min,x4max,y4min,y4max;
	float x5min,x5max,y5min,y5max;
	float x6min,x6max,y6min,y6max;

	usghed usgh;


	
	initargs(argc,argv);
        askdoc(1);

	/* required parameters */
	if (!getparfloat("xqc",&xqc)) err(" xqc missing ");
	if (!getparfloat("yqc",&yqc)) err(" yqc missing ");

	/* optional parameters */

	if (!getparint("ixy",&ixy)) ixy = 0;

	if (!getparstring("vs3d1",&in1)) {
		i1 = 0;
	} else {
		i1 = 1;
	}

	if (!getparstring("vs3d2",&in2)) {
		i2 = 0; 
	} else {
		i2 = 1;
	}

	if (!getparstring("vs3d3",&in3)) {
		i3 = 0; 
	} else {
		i3 = 1;
	}

	if (!getparstring("vgrid1",&vgrid1)) {
		iv1 = 0; 
	} else {
		iv1 = 1;
		vfp1 = efopen(vgrid1,"r");
	}
	if (!getparstring("vgrid2",&vgrid2)) {
		iv2 = 0; 
	} else {
		iv2 = 1;
		vfp2 = efopen(vgrid2,"r");
	}
	if (!getparstring("vgrid3",&vgrid3)) {
		iv3 = 0; 
	} else {
		iv3 = 1;
		vfp3 = efopen(vgrid3,"r");
	}

	if (!getparfloat("dt",&dt)) dt = 100.;
	if (!getparint("vgtype1",&vgtype1)) vgtype1 = 2;
	if (!getparint("vgtype2",&vgtype2)) vgtype2 = 2;
	if (!getparint("vgtype3",&vgtype3)) vgtype3 = 2;
	if (!getparint("vstype1",&vstype1)) vstype1 = 0;
	if (!getparint("vstype2",&vstype2)) vstype2 = 0;
	if (!getparint("vstype3",&vstype3)) vstype3 = 0;
	if (!getparint("paper",&paper)) paper = 0;
	if (!getparint("vauxplot",&vauxplot)) vauxplot = 0;
	if (!getparint("viplot",&viplot)) viplot = 1;
	ia = vauxplot;
	if(ia==2) ia = 1;
	ii = viplot;
	if (!getparint("nvfmax",&nvfmax)) nvfmax = 4096;
	if (!getparint("ntvmax",&ntvmax)) ntvmax = 256;
	if (!getparfloat("tmax",&tmax)) tmax = 100000.;
	if (!getparstring("title",&title)) {
		title = (char*) malloc(30*sizeof(char));
		strcpy(title,"velocity qc ");
	}		

	time = (float*) malloc(ntmax*sizeof(float));
	work = (float*) malloc(4*ntmax*sizeof(float));
	t1 = (float*) malloc(ntmax*sizeof(float));
	t2 = (float*) malloc(ntmax*sizeof(float));
	t3 = (float*) malloc(ntmax*sizeof(float));
	t4 = (float*) malloc(ntmax*sizeof(float));
	t5 = (float*) malloc(ntmax*sizeof(float));
	t6 = (float*) malloc(ntmax*sizeof(float));
	v1 = (float*) malloc(ntmax*sizeof(float));
	v2 = (float*) malloc(ntmax*sizeof(float));
	v3 = (float*) malloc(ntmax*sizeof(float));
	v4 = (float*) malloc(ntmax*sizeof(float));
	v5 = (float*) malloc(ntmax*sizeof(float));
	v6 = (float*) malloc(ntmax*sizeof(float));

	/* read in t, v at closest location of data in1 */

	n1s = 0;
	if(i1==1) {
		findvf(in1, xqc, yqc, t1, v1, &np1, &x1, &y1, ntvmax, nvfmax); 
		times(dt, t1, np1, time, &nt, tmax);
		t1s = (float*) malloc(nt*sizeof(float));
		v1s = (float*) malloc(nt*sizeof(float));
		bcopy(time,t1s,nt*sizeof(float));
		bcopy(v1,v1s,nt*sizeof(float));
		if(vstype1==0) {
			vrms2vint(t1, v1, np1, t1s, v1s, nt);
		} else if(vstype1==1) {
			vavg2vint(t1, v1, np1, t1s, v1s, nt);
		}
		n1s = nt;
	}

	n2s = 0;
	if(i2==1) {
		findvf(in2, xqc, yqc, t2, v2, &np2, &x2, &y2, ntvmax, nvfmax); 
		times(dt, t2, np2, time, &nt, tmax);
		t2s = (float*) malloc(nt*sizeof(float));
		v2s = (float*) malloc(nt*sizeof(float));
		bcopy(time,t2s,nt*sizeof(float));
		bcopy(v2,v2s,nt*sizeof(float));
		if(vstype2==0) {
			vrms2vint(t2, v2, np2, t2s, v2s, nt);
		} else if (vstype2==1) {
			vavg2vint(t2, v2, np2, t2s, v2s, nt);
		}
		n2s = nt;
	}

	n3s = 0;
	if(i3==1) { 
		findvf(in3, xqc, yqc, t3, v3, &np3, &x3, &y3, ntvmax, nvfmax); 
		times(dt, t3, np3, time, &nt, tmax);
		t3s = (float*) malloc(nt*sizeof(float));
		v3s = (float*) malloc(nt*sizeof(float));
		bcopy(time,t3s,nt*sizeof(float));
		bcopy(v3,v3s,nt*sizeof(float));
		if(vstype3==0) {
			vrms2vint(t3, v3, np3, t3s, v3s, nt);
		} else if (vstype3==1) {
			vavg2vint(t3, v3, np3, t3s, v3s, nt);
		}
		n3s = nt;
	}

	/* read in volocity function from vgrid */
	n4s = 0;
	if(iv1==1) {
		fgetusghdr(vfp1, &usgh);
		if(ixy==0) {
			tmp = (xqc - usgh.o2)/usgh.d2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.o3)/usgh.d3 + .5;
			iy = tmp;
		} else {
			if(usgh.dcdp2==0 || usgh.dline3==0)
				err(" check $s gridheader ",vgrid1);
			tmp = (xqc - usgh.ocdp2)/usgh.dcdp2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.oline3)/usgh.dline3 + .5;
			iy = tmp;
		}
		if(ix<0) ix = 0;
		if(ix>=usgh.n2) ix = usgh.n2 - 1;
		if(iy<0) iy = 0;
		if(iy>=usgh.n3) iy = usgh.n3 - 1;

		if(ixy==0) {
			x4 = usgh.o2 + ix*usgh.d2;
			y4 = usgh.o3 + iy*usgh.d3;
		} else {
			x4 = usgh.ocdp2 + ix*usgh.dcdp2;
			y4 = usgh.oline3 + iy*usgh.dline3;
		}

		fprintf(stderr,"Velociy grid %s: x=%g y=%g t0=%g dt=%g nt=%d\n",
				vgrid1,x4,y4,usgh.o1,usgh.d1,usgh.n1);

		iseek = (iy*usgh.n2+ix)*usgh.n1*usgh.dtype;
		efseek(vfp1,iseek,0);
		efread(v4,sizeof(float),usgh.n1,vfp1);
		for(it=0;it<usgh.n1;it++) t4[it] = usgh.o1 + it*usgh.d1;
		np4 = usgh.n1; 
		if(vgtype1==2) {
			if(tmax < 100000.) {
				tmp = (tmax - usgh.o1)/usgh.d1;
				nt = tmp;
			} else {
				nt = np4;
			}
			t4s = (float*) malloc(nt*sizeof(float));
			v4s = (float*) malloc(nt*sizeof(float));
			bcopy(t4,t4s,nt*sizeof(float));
			bcopy(v4,v4s,nt*sizeof(float));
		} else if(vgtype1==0) {
			times(dt, t4, np4, time, &nt, tmax);
			t4s = (float*) malloc(nt*sizeof(float));
			v4s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t4s,nt*sizeof(float));
			vrms2vint(t4, v4, np4, t4s, v4s, nt);
		} else if(vgtype1=1) {
			times(dt, t4, np4, time, &nt, tmax);
			t4s = (float*) malloc(nt*sizeof(float));
			v4s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t4s,nt*sizeof(float));
			vavg2vint(t4, v4, np4, t4s, v4s, nt);
		}
		n4s = nt;
		if(ixy==0) {
			x4min = usgh.o2;
			x4max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y4min = usgh.o3;
			y4max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		} else {
			x4min = usgh.o2;
			x4max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y4min = usgh.o3;
			y4max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		}
	}

	n5s = 0;
	if(iv2==1) {
		fgetusghdr(vfp2, &usgh);

		if(ixy==0) {
			tmp = (xqc - usgh.o2)/usgh.d2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.o3)/usgh.d3 + .5;
			iy = tmp;
		} else {
			if(usgh.dcdp2==0 || usgh.dline3==0)
				err(" check $s gridheader ",vgrid2);
			tmp = (xqc - usgh.ocdp2)/usgh.dcdp2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.oline3)/usgh.dline3 + .5;
			iy = tmp;
		}
		if(ix<0) ix = 0;
		if(ix>=usgh.n2) ix = usgh.n2 - 1;

		if(ixy==0) {
			x5 = usgh.o2 + ix*usgh.d2;
			y5 = usgh.o3 + iy*usgh.d3;
		} else {
			x5 = usgh.ocdp2 + ix*usgh.dcdp2;
			y5 = usgh.oline3 + iy*usgh.dline3;
		}

		fprintf(stderr,"Velociy grid %s: x=%g y=%g t0=%g dt=%g nt=%d\n",
				vgrid2,x5,y5,usgh.o1,usgh.d1,usgh.n1);
		iseek = (iy*usgh.n2+ix)*usgh.n1*usgh.dtype;
		efseek(vfp2,iseek,0);
		efread(v5,sizeof(float),usgh.n1,vfp2);
		for(it=0;it<usgh.n1;it++) t5[it] = usgh.o1 + it*usgh.d1;
		np5 = usgh.n1; 
		if(vgtype2==2) {
			if(tmax < 100000.) {
				tmp = (tmax - usgh.o1)/usgh.d1;
				nt = tmp;
			} else {
				nt = np5;
			}
			t5s = (float*) malloc(nt*sizeof(float));
			v5s = (float*) malloc(nt*sizeof(float));
			bcopy(t5,t5s,nt*sizeof(float));
			bcopy(v5,v5s,nt*sizeof(float));
		} else if(vgtype2==0) {
			times(dt, t5, np5, time, &nt, tmax);
			t5s = (float*) malloc(nt*sizeof(float));
			v5s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t5s,nt*sizeof(float));
			vrms2vint(t5, v5, np5, t5s, v5s, nt);
		} else if(vgtype2==1) {
			times(dt, t5, np5, time, &nt, tmax);
			t5s = (float*) malloc(nt*sizeof(float));
			v5s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t5s,nt*sizeof(float));
			vavg2vint(t5, v5, np5, t5s, v5s, nt);
		}
		n5s = nt;

		if(ixy==0) {
			x5min = usgh.o2;
			x5max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y5min = usgh.o3;
			y5max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		} else {
			x5min = usgh.o2;
			x5max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y5min = usgh.o3;
			y5max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		}

	}

	n6s = 0;
	if(iv3==1) {
		fgetusghdr(vfp3, &usgh);

		if(ixy==0) {
			tmp = (xqc - usgh.o2)/usgh.d2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.o3)/usgh.d3 + .5;
			iy = tmp;
		} else {
			if(usgh.dcdp2==0 || usgh.dline3==0)
				err(" check $s gridheader ",vgrid3);
			tmp = (xqc - usgh.ocdp2)/usgh.dcdp2 + .5;
			ix = tmp;
			tmp = (yqc - usgh.oline3)/usgh.dline3 + .5;
			iy = tmp;
		}

		if(ix<0) ix = 0;
		if(ix>=usgh.n2) ix = usgh.n2 - 1;
		if(iy<0) iy = 0;
		if(iy>=usgh.n3) iy = usgh.n3 - 1;

		if(ixy==0) {
			x6 = usgh.o2 + ix*usgh.d2;
			y6 = usgh.o3 + iy*usgh.d3;
		} else {
			x6 = usgh.ocdp2 + ix*usgh.dcdp2;
			y6 = usgh.oline3 + iy*usgh.dline3;
		}

		fprintf(stderr,"Velociy grid %s: x=%g y=%g t0=%g dt=%g nt=%d\n",
				vgrid3,x6,y6,usgh.o1,usgh.d1,usgh.n1);

		iseek = (iy*usgh.n2+ix)*usgh.n1*usgh.dtype;
		efseek(vfp3,iseek,0);
		efread(v6,sizeof(float),usgh.n1,vfp3);
		for(it=0;it<usgh.n1;it++) t6[it] = usgh.o1 + it*usgh.d1;
		np6 = usgh.n1; 
		if(vgtype3==2) {
			if(tmax < 100000.) {
				tmp = (tmax - usgh.o1)/usgh.d1;
				nt = tmp;
			} else {
				nt = np6;
			}
			t6s = (float*) malloc(nt*sizeof(float));
			v6s = (float*) malloc(nt*sizeof(float));
			bcopy(t6,t6s,nt*sizeof(float));
			bcopy(v6,v6s,nt*sizeof(float));
		} else if(vgtype3==0) {
			times(dt, t6, np6, time, &nt, tmax);
			t6s = (float*) malloc(nt*sizeof(float));
			v6s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t6s,nt*sizeof(float));
			vrms2vint(t6, v6, np6, t6s, v6s, nt);
		} else if(vgtype3==1) {
			times(dt, t6, np6, time, &nt, tmax);
			t6s = (float*) malloc(nt*sizeof(float));
			v6s = (float*) malloc(nt*sizeof(float));
			bcopy(time,t6s,nt*sizeof(float));
			vavg2vint(t6, v6, np6, t6s, v6s, nt);
		}
		n6s = nt;

		if(ixy==0) {