qcvs3d.c

seismic software,very useful

			x6min = usgh.o2;
			x6max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y6min = usgh.o3;
			y6max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		} else {
			x6min = usgh.o2;
			x6max = usgh.o2 + (usgh.n2 - 1)*usgh.d2;
			y6min = usgh.o3;
			y6max = usgh.o3 + (usgh.n3 - 1)*usgh.d3;
		}
	}

	bzero(cmd,1024);
	bzero(cname,80);
	sprintf(cname,"coordinate.%d\0",getpid());

	if(paper==0) {
        	sprintf(cmd, "( <%s xgraph n=%d,%d,%d,%d,%d,%d,0,0,0,%d,%d,%d nplot=12 marksize=12,12,12,12,12,12,12,12,12,12,12,12 mark=0,1,2,3,4,5,0,1,2,3,4,5 label1=\"x (inline)\" label2=\"y (crossline)\" width=400 height=400 grid1=dash grid2=dash title=\"+=vs3d1 *=vs3d2 x=vs3d3 tri=vgrid1 sq=vgrid2 o=vgrid3\" ; /bin/rm -f %s ) & ", cname, i1, i2, i3, 5*iv1, 5*iv2, 5*iv3, iv1, iv2, iv3, cname);
	} else {
        	sprintf(cmd, "( <%s psgraph n=%d,%d,%d,%d,%d,%d,%d,%d,%d nplot=9 mark=0,1,2,3,4,5,3,4,5 marksize=12,12,12,12,12,12,12,12,12 label1=\"x (inline)\" label2=\"y (crossline)\" wbox=5 hbox=5 grid1=dash grid2=dash title=\"+=vs3d1 *=vs3d2 x=vs3d3 tri=vgrid1 sq=vgrid2 o=vgrid3\" labelsize=14 titlesize=16 | lpr ; /bin/rm -f %s ) & ", cname, i1, i2, i3, 5*iv1, 5*iv2, 5*iv3, iv1, iv2, iv3, cname);
	}


	datafp = efopen(cname,"w");

	if(i1==1) {
		efwrite(&x1,sizeof(float),1,datafp);
		efwrite(&y1,sizeof(float),1,datafp);
	}
	if(i2==1) {
		efwrite(&x2,sizeof(float),1,datafp);
		efwrite(&y2,sizeof(float),1,datafp);
	}
	if(i3==1) {
		efwrite(&x3,sizeof(float),1,datafp);
		efwrite(&y3,sizeof(float),1,datafp);
	}
	if(iv1==1) {
		efwrite(&x4min,sizeof(float),1,datafp);
		efwrite(&y4min,sizeof(float),1,datafp);
		efwrite(&x4max,sizeof(float),1,datafp);
		efwrite(&y4min,sizeof(float),1,datafp);
		efwrite(&x4max,sizeof(float),1,datafp);
		efwrite(&y4max,sizeof(float),1,datafp);
		efwrite(&x4min,sizeof(float),1,datafp);
		efwrite(&y4max,sizeof(float),1,datafp);
		efwrite(&x4min,sizeof(float),1,datafp);
		efwrite(&y4min,sizeof(float),1,datafp);
	}
	if(iv2==1) {
		efwrite(&x5min,sizeof(float),1,datafp);
		efwrite(&y5min,sizeof(float),1,datafp);
		efwrite(&x5max,sizeof(float),1,datafp);
		efwrite(&y5min,sizeof(float),1,datafp);
		efwrite(&x5max,sizeof(float),1,datafp);
		efwrite(&y5max,sizeof(float),1,datafp);
		efwrite(&x5min,sizeof(float),1,datafp);
		efwrite(&y5max,sizeof(float),1,datafp);
		efwrite(&x5min,sizeof(float),1,datafp);
		efwrite(&y5min,sizeof(float),1,datafp);
	}
	if(iv3==1) {
		efwrite(&x6min,sizeof(float),1,datafp);
		efwrite(&y6min,sizeof(float),1,datafp);
		efwrite(&x6max,sizeof(float),1,datafp);
		efwrite(&y6min,sizeof(float),1,datafp);
		efwrite(&x6max,sizeof(float),1,datafp);
		efwrite(&y6max,sizeof(float),1,datafp);
		efwrite(&x6min,sizeof(float),1,datafp);
		efwrite(&y6max,sizeof(float),1,datafp);
		efwrite(&x6min,sizeof(float),1,datafp);
		efwrite(&y6min,sizeof(float),1,datafp);
	}
	if(iv1==1) {
		efwrite(&x4,sizeof(float),1,datafp);
		efwrite(&y4,sizeof(float),1,datafp);
	if(iv2==1) {
		efwrite(&x5,sizeof(float),1,datafp);
		efwrite(&y5,sizeof(float),1,datafp);
	}
	if(iv3==1) {
		efwrite(&x6,sizeof(float),1,datafp);
		efwrite(&y6,sizeof(float),1,datafp);
	}
	}
	efclose(datafp);

	fprintf(stderr,"\n");
	fprintf(stderr,"Displaying Command: \n");
	fprintf(stderr,"%s \n",cmd);

	system(cmd);
	bzero(cname,80);
	sprintf(cname,"data.to.plot.%d\0",getpid());
	datafp = efopen(cname,"w");
	bzero(cmd,1024);
	
	if(paper==0) {
        	sprintf(cmd, "( <%s xgraph n=%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d nplot=12 marksize=6,6,6,6,6,6,6,6,6,6,6,6 mark=0,1,2,3,4,5,0,1,2,3,4,5 linecolor=2,3,7,5,6,4,2,3,7,5,6,4 style=seismic label1=time label2=velocity width=900 height=1000 grid1=dash grid2=dash title=\"%s at x=%g y=%g\" ; /bin/rm -f %s ) & ", cname,ii*2*n1s,ii*2*n2s,ii*2*n3s,ii*2*n4s,ii*2*n5s,ii*2*n6s,ia*n1s,ia*n2s,ia*n3s,ia*n4s,ia*n5s,ia*n6s,title,xqc,yqc,cname);
	} else {
        	sprintf(cmd, "( <%s psgraph n=%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d nplot=12 marksize=6,6,6,6,6,6,6,6,6,6,6,6 style=seismic label1=time label2=velocity wbox=6 hbox=8 grid1=dash grid2=dash title=\"%s at x=%g y=%g\" titlesize=16 labelsize=14 | lpr ; /bin/rm -f %s ) & ", cname,ii*2*n1s,ii*2*n2s,ii*2*n3s,ii*2*n4s,ii*2*n5s,ii*2*n6s,ia*n1s,ia*n2s,ia*n3s,ia*n4s,ia*n5s,ia*n6s,title,xqc,yqc,cname);
	}

	fprintf(stderr,"%s \n",cmd);

	free(work);
	work = (float*) malloc(n1s*4*sizeof(float));
	
	for(it=0;it<n1s;it++) {
		work[it*4] = t1s[it];
		work[it*4+1] = v1s[it];
		work[it*4+2] = t1s[it];
		if(it==n1s-1) {
			work[it*4+3] = v1s[it];
		} else {
			work[it*4+3] = v1s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n1s,datafp);

	free(work);
	work = (float*) malloc(n2s*4*sizeof(float));
	for(it=0;it<n2s;it++) {
		work[it*4] = t2s[it];
		work[it*4+1] = v2s[it];
		work[it*4+2] = t2s[it];
		if(it==n2s-1) {
			work[it*4+3] = v2s[it];
		} else {
			work[it*4+3] = v2s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n2s,datafp);


	free(work);
	work = (float*) malloc(n3s*4*sizeof(float));
	for(it=0;it<n3s;it++) {
		work[it*4] = t3s[it];
		work[it*4+1] = v3s[it];
		work[it*4+2] = t3s[it];
		if(it==n3s-1) {
			work[it*4+3] = v3s[it];
		} else {
			work[it*4+3] = v3s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n3s,datafp);

	free(work);
	work = (float*) malloc(n4s*4*sizeof(float));
	for(it=0;it<n4s;it++) {
		work[it*4] = t4s[it];
		work[it*4+1] = v4s[it];
		work[it*4+2] = t4s[it];
		if(it==n4s-1) {
			work[it*4+3] = v4s[it];
		} else {
			work[it*4+3] = v4s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n4s,datafp);

	free(work);
	work = (float*) malloc(n5s*4*sizeof(float));
	for(it=0;it<n5s;it++) {
		work[it*4] = t5s[it];
		work[it*4+1] = v5s[it];
		work[it*4+2] = t5s[it];
		if(it==n5s-1) {
			work[it*4+3] = v5s[it];
		} else {
			work[it*4+3] = v5s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n5s,datafp);

	free(work);
	work = (float*) malloc(n6s*4*sizeof(float));
	for(it=0;it<n6s;it++) {
		work[it*4] = t6s[it];
		work[it*4+1] = v6s[it];
		work[it*4+2] = t6s[it];
		if(it==n6s-1) {
			work[it*4+3] = v6s[it];
		} else {
			work[it*4+3] = v6s[it+1];
		}
	}
	if(viplot==1) efwrite(work,sizeof(float),4*n6s,datafp);

	if(vauxplot==1) {
		free(work);
		work = (float*) malloc(2*ntmax*sizeof(float));
		if(n1s>0) vrmsout(datafp,t1s,v1s,work,n1s);
		if(n2s>0) vrmsout(datafp,t2s,v2s,work,n2s);
		if(n3s>0) vrmsout(datafp,t3s,v3s,work,n3s);
		if(n4s>0) vrmsout(datafp,t4s,v4s,work,n4s);
		if(n5s>0) vrmsout(datafp,t5s,v5s,work,n5s);
		if(n6s>0) vrmsout(datafp,t6s,v6s,work,n6s);
	
	} else if(vauxplot==2) {
		free(work);
		work = (float*) malloc(2*ntmax*sizeof(float));
		if(n1s>0) vavgout(datafp,t1s,v1s,work,n1s);
		if(n2s>0) vavgout(datafp,t2s,v2s,work,n2s);
		if(n3s>0) vavgout(datafp,t3s,v3s,work,n3s);
		if(n4s>0) vavgout(datafp,t4s,v4s,work,n4s);
		if(n5s>0) vavgout(datafp,t5s,v5s,work,n5s);
		if(n6s>0) vavgout(datafp,t6s,v6s,work,n6s);
	}

	efclose(datafp);

	system(cmd);

	exit(0);

}

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

	int it;

	vint2vrms(ts,vs,work,ns);
	bcopy(work,vs,ns*sizeof(float));
	for(it=0;it<ns;it++) {
		work[it*2] = ts[it];
		work[it*2+1] = vs[it];
	} 
	efwrite(work,sizeof(float),2*ns,datafp);
}

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

	int it;

	vint2vavg(ts,vs,work,ns);
	bcopy(work,vs,ns*sizeof(float));
	for(it=0;it<ns;it++) {
		work[it*2] = ts[it];
		work[it*2+1] = vs[it];
	} 
	efwrite(work,sizeof(float),2*ns,datafp);
}
	


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

	float tmp, ttmax;
	int i1, i;

	/* find out the maximum time */
	ttmax = t[ntv-1];
	if(ttmax>tmax) ttmax = tmax;

	if(dt>0.) {
		tmp = ttmax/dt + 1.5;
		*nt = tmp;
		for(i1=0;i1<*nt;i1++) time[i1] = i1*dt;
	} else {
		i = 0;
		for(i1=0;i1<ntv;i1++) {
			if(t[i1]<=tmax) {
				time[i] = t[i1];
				i = i + 1;
			}
		}
		*nt = i;
	}
}

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

	int i1;
	
	vrms[0] = vi[0]*vi[0]*t[0];
	for(i1=1;i1<nt;i1++) {
		vrms[i1] = vrms[i1-1] + vi[i1]*vi[i1]*(t[i1]-t[i1-1]);
	}
	for(i1=1;i1<nt;i1++) {
		vrms[i1] = sqrt(vrms[i1]/t[i1]);
	}
	vrms[0] = vi[0];
}

void vint2vavg(float *t, float *vi, float *vavg, int nt) {

	int i1;
	
	vavg[0] = vi[0]*t[0];
	for(i1=1;i1<nt;i1++) {
		vavg[i1] = vavg[i1-1] + vi[i1]*(t[i1]-t[i1-1]);
	}
	for(i1=1;i1<nt;i1++) {
		vavg[i1] = vavg[i1]/t[i1];
	}
	vavg[0] = vi[0];
}


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

	float *work, tmax;
	int *indx;
	int i1, j;


	work = (float*) emalloc(nt*sizeof(float));
	indx = (int*) emalloc(nt*sizeof(int));


	/* linearly interpolate v to nt times */
	lin1d_(t,v,&ntv,time,vi,&nt,indx);

	/* compute interval velocities via dix formula */
	work[0] = vi[0]*vi[0];
	tmax = t[ntv-1];

	for(j=1;j<nt;j++) {
		if(time[j]<=tmax) { 
			work[j] = (vi[j]*vi[j]*time[j] -  	
		   		vi[j-1]*vi[j-1]*time[j-1]) / 
		  		(time[j]-time[j-1]);
		} else {
			work[j] = work[j-1];
		}

	}
	for(j=0;j<nt;j++) {
		if(work[j]<0.) err(" Negative Interval Velocity "); 
		vi[j] = sqrt(work[j]);
	}

	free(indx);
	free(work);

}

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

	float *work, tmax;
	int *indx;
	int i1, j;


	work = (float*) emalloc(nt*sizeof(float));
	indx = (int*) emalloc(nt*sizeof(int));


	/* linearly interpolate v to nt times */
	lin1d_(t,v,&ntv,time,vi,&nt,indx);

	/* compute interval velocities via dix formula */
	work[0] = vi[0];
	tmax = t[ntv-1];

	for(j=1;j<nt;j++) {
		if(time[j]<=tmax) { 
			work[j] = (vi[j]*time[j] -  
		   		vi[j-1]*time[j-1]) / 
		  		(time[j]-time[j-1]);
		} else {
			work[j] = work[j-1];
		}

	}
	for(j=0;j<nt;j++) {
		if(work[j]<0.) err(" Negative Interval Velocity "); 
		vi[j] = work[j];
	}

	free(indx);
	free(work);

}

/* find closest VS3D card to location (xqc,yqc) */ 
/* t and v are pre-allocated arrays of length longer than 256 */
/*  
    input:
	in= 	file name of vs3d cards
	xqc=    x coordinate of qc position
	yqc=    y coordinate of qc position
	ntvmax= maximum number of t-v pairs per velocity function
	nvfmax= maximum number of velocity functions in the input
    output:
	t=      times of closest vs3d card
	v=      velocities of closest vs3d card
	np=     number of t-v pairs at this locations
	xx=     x coorinate of the closest vs3d card
	yy=     y coorinate of the closest vs3d card

    author:		z. li		      		9/14/93 
*/
	

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


	int n1, n2, nxy;
    	float *xs, *ys, *tpicks, *vpicks;
    	int *nps;
	FILE *infp;
	float dis, x, y, tmp;
	int ip, i1, idis;

	infp = efopen(in,"r");

    	/* at most 4096 input (x,y) VS3D cards with at most 256 time-vel 
		pairs each */

    	n1 = ntvmax;
    	n2 = nvfmax;

    	/* arrays used to store all VELF card's cdp, time and velocity */
    	xs = (float*)malloc(n2*sizeof(int));
    	ys = (float*)malloc(n2*sizeof(int));
    	tpicks = (float*)malloc(n1*n2*sizeof(float));
    	vpicks = (float*)malloc(n1*n2*sizeof(float));
    	nps = (int*)malloc(n2*sizeof(int));
	bzero(nps,n2*sizeof(int));

    	/* read in VS3D cards */
    	nxy = 0;
	vs3dread(infp,xs,ys,tpicks,vpicks,&nxy,nps,n1,n2);
   	if (nxy==0) err("No VS3D card input for in=%s ! Job aborted",in);
	x = xs[0];
	y = ys[0];
	dis = sqrt((xqc-x)*(xqc-x)+(yqc-y)*(yqc-y));
	idis = 0;
	*np = nps[0];
	*xx = x;
	*yy = y;
	for(i1=0;i1<*np;i1++) {
		t[i1] = tpicks[i1];
		v[i1] = vpicks[i1];
	}

	for(ip=1;ip<nxy;ip++) {
		x = xs[ip];
		y = ys[ip];
		tmp = sqrt((xqc-x)*(xqc-x)+(yqc-y)*(yqc-y));
		if(tmp<dis) {
			dis = tmp;
			idis = ip;
		}
    }
	*np = nps[idis];
	for(i1=0;i1<*np;i1++) {
		t[i1] = tpicks[idis*n1+i1];
		v[i1] = vpicks[idis*n1+i1];
	}
	*xx = xs[idis];
	*yy = ys[idis];
			
	fprintf(stderr,
	"VS3D card found for %s at x=%g y=%g ip=%d location for %d t-v pairs \n",
			in,*xx,*yy,idis+1,*np); 

	free(nps);
	free(tpicks);
	free(vpicks);	
	free(xs);
	free(ys);
	fclose(infp);

}