velo2ifile.c

seismic software,very useful

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


	/* output picks */

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

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

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

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

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

   			fprintf(ifilefp, 
"%9.1f %9.1f %8.0f %8.0f   %1d   %4.2f  %4.0f  %4.2f %8.0f %8.3f \n",
				xpicks[ip+h0],zpicks[ip+h0],
				veltop[ip+h0],velbot[ip+h0],dif[ip+h0],
				dens[ip+h0],qval[ip+h0],pois[ip+h0],
				velavg[ip+h0],dvdz[ip+h0]);
		}
	}
}


/* compute velocities from VELO and DVDZ cards */
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 *veltops,float *velbots,float *velavgs,float *dvdzs,
     int *vtopupdate, int *vbotupdate, int *hnums) {

	float x, z, va, g, vat, gt, vab, gb;
	int ih, ip, jh;
	float zmax, xmin, xmax, zt, zb, tmp;
	float *vth, *vbh, *vah, *dvdzh, *xvh, *zvs, *zvh;
	float *zgrid, *vagrid, *dvdzgrid, dz, *g1, *g2, *z2, z21, z2n;
	int i1,i2,n1,n2,one=1,nzgrid;
	int *ivs, *nvh;
	int iz, n, j2, j1, nvs;
	float *xnew, *znew, *vanew, *vtnew, *vbnew, *gnew, xpre, xnow;
	int *indx, inow;

	float *zsort;
	int *isort;
	char vtopyes, vbotyes;

	int it1, it2, ib1, ib2, ii;
	
	
	if(ncdp_velo==0) {
		/* if no velo card input, simply return with warning message */
		warn("no VELO card input ! \n");
	} else {
		/* find maximum depth */
		zmax = 0.;
		for(ih=0;ih<nhs;ih++) {
			for(iz=0;iz<npicks[ih];iz++) {
				if(zmax<zpicks[ih*maxpks+iz])
				  	zmax=zpicks[ih*maxpks+iz];
			}
		}


		for(i2=0;i2<ncdp_velo;i2++) {
			for(iz=0;iz<nps_velo[i2];iz++) {
				if(zmax<z_velo[i2*maxp+iz])
				  	zmax=z_velo[i2*maxp+iz];
			}
		}




		nzgrid = 1024;
		zgrid = (float *) malloc(nzgrid*sizeof(float));
		vagrid = (float *) malloc(nzgrid*ncdp_velo*sizeof(float));
		dvdzgrid = (float *) malloc(nzgrid*ncdp_velo*sizeof(float));
		g1 = (float *) malloc(nzgrid*sizeof(float));
		g2 = (float *) malloc(nzgrid*sizeof(float));

		vth = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		vbh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		vah = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		dvdzh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		xvh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		zvh = (float *) malloc(nhs*ncdp_velo*sizeof(float));
		nvh = (int *) malloc(nhs*sizeof(int));

		zvs = (float *) malloc(nhs*sizeof(float));
		ivs = (int *) malloc(nhs*sizeof(int));
		zsort = (float *) malloc(nhs*sizeof(float));
		isort = (int *) malloc(nhs*sizeof(int));

		vtnew = (float *) malloc(maxpks*sizeof(float));
		vbnew = (float *) malloc(maxpks*sizeof(float));
		indx = (int *) malloc(maxpks*sizeof(int));
		vanew = (float *) malloc(maxpks*sizeof(float));
		gnew = (float *) malloc(maxpks*sizeof(float));
		xnew = (float *) malloc(maxpks*sizeof(float));
		znew = (float *) malloc(maxpks*sizeof(float));

		z2 = (float *) malloc(nhs*sizeof(float));
		/* spline interpolation of velocities --- not a good idea 
		vt2 = (float *) malloc(maxpks*sizeof(float));
		vb2 = (float *) malloc(maxpks*sizeof(float));
		va2 = (float *) malloc(maxpks*sizeof(float));
		dvdz2 = (float *) malloc(maxpks*sizeof(float));
		*/

	

		dz = zmax/(nzgrid-2);
		for(iz=0;iz<nzgrid;iz++) zgrid[iz] = iz * dz;
		/* compute average velocity and gradient at velo location and
		   at depth zgrid */


		for(i2=0;i2<ncdp_velo;i2++) {
			n2 = i2*nzgrid;
			/* find gradient */
			x = x_velo[i2];

			if(ncdp_dvdz==0) {
				for(iz=0;iz<nzgrid;iz++) 
					dvdzgrid[iz+n2] = 0.;
			} else if (ncdp_dvdz==1) {
				n = nps_dvdz[0];
				dvdzint(zt_dvdz,zb_dvdz,dvdz_dvdz,n,
					zgrid,dvdzgrid+n2,nzgrid);
			} else {
				bisear_(&ncdp_dvdz,&one,x_dvdz,&x,&j2);
				if(x<=x_dvdz[0]) {
					n = nps_dvdz[0];
					dvdzint(zt_dvdz,zb_dvdz,dvdz_dvdz,n,
					     zgrid,dvdzgrid+n2,nzgrid);
				} else if (x>=x_dvdz[ncdp_dvdz-1]) {
					n = nps_dvdz[ncdp_dvdz-1];
					i1 = (ncdp_dvdz-1)*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,zgrid,
						dvdzgrid+n2,nzgrid);
				} else {
					j1 = j2 - 1;
					n = nps_dvdz[j1];
					i1 = j1*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,
						zgrid,g1,nzgrid);
					n = nps_dvdz[j2];
					i1 = j2*maxp;
					dvdzint(zt_dvdz+i1,zb_dvdz+i1,
						dvdz_dvdz+i1,n,
						zgrid,g2,nzgrid);
					tmp = (x-x_dvdz[j1]) /
					      (x_dvdz[j2]-x_dvdz[j1]);
				
					for(iz=0;iz<nzgrid;iz++) {
						dvdzgrid[iz+n2] = g1[iz]+
						tmp * (g2[iz]-g1[iz]);
					}
				}
			}
			/* find average velocity at zgrid */
			n = nps_velo[i2];
			n1 = i2*maxp;
			j1 = 0;
			j2 = 0;
			vconint(z_velo+n1,v_velo+n1,n,zgrid,vagrid+n2,nzgrid,
				j1,j2,dvdzgrid+n2);

		}
			


		/* compute vt and vb along horizons at velan locations */

		for(ih=0;ih<nhs;ih++) nvh[ih] = 0;
		for(i2=0;i2<ncdp_velo;i2++) {
			x = x_velo[i2];	
			jh = 0;
			/* find cross points of horizons at velan locations */

			for(ih=0;ih<nhs;ih++) {
				xmin = xpicks[ih*maxpks]-xtol;
				xmax = xpicks[ih*maxpks+npicks[ih]-1]+xtol;
				n1 = npicks[ih]; 
				n2 = ih * maxpks;
				if(x>=xmin && x<=xmax) {
					ivs[jh] = ih;
					/* find z */
					bisear_(&n1,&one,xpicks+n2,&x,&i1);
					if(x<=xpicks[n2]) {
						z = zpicks[n2];
					} else if(x>=xpicks[n1+n2-1]) {
						z = zpicks[n1+n2-1];
					} else {
						/*
						if(n1==2) {
						*/
						/* linear interpolation */
						   z = zpicks[n2+i1-1] +
						    (x-xpicks[n2+i1-1])*
						    (zpicks[n2+i1]-
						     zpicks[n2+i1-1])/
						    (xpicks[n2+i1]-
						     xpicks[n2+i1-1]);
						/*
						} else {
						*/
						/* spline interpolation */
						/* not a good idea */
						/*
						   z21 = 0.;
                                 		   z2n = 0.;
                                                   spline_(xpicks+n2,zpicks+n2,
                                         		   &n1,&z21,&z2n,z2);
                                 		   splint_(xpicks+n2,zpicks+n2,
							   z2,&n1,&x,&z);
						}
						*/
					}
					zvs[jh] = z; 
					jh = jh + 1;
				}
			}

			/* sort zvs into ascending order */
			for (iz=0;iz<jh;iz++) {
				isort[iz] = iz;
				zsort[iz] = zvs[iz];
			}

			qkisort(jh,zsort,isort);
			for(iz=0;iz<jh;iz++) zvs[iz] = zsort[isort[iz]];

			/* find vt and vb */
			for(iz=0;iz<jh;iz++) {

				z = zvs[iz];
				tmp = z/dz;
				i1 = tmp;
				it1 = i1 - 2;
				it2 = i1 - 1; 
				ib1 = i1;
				ib2 = i1 + 1;

				ilimit(&it1,0,nzgrid-1);
				ilimit(&it2,0,nzgrid-1);
				ilimit(&ib1,0,nzgrid-1);
				ilimit(&ib1,0,nzgrid-1);
				ilimit(&i1,0,nzgrid-1);
				
				ii = ivs[isort[iz]]*ncdp_velo
					+nvh[ivs[isort[iz]]];

				
				z = zgrid[it2];
				zt = zgrid[it1];
				va = vagrid[i2*nzgrid+it2];
				vat = vagrid[i2*nzgrid+it1];
				gt = dvdzgrid[i2*nzgrid+it2];

				vth[ii] = (z*va-zt*vat)/(z-zt) + 0.5*gt*(z-zt);

				z = zgrid[ib1];
				zb = zgrid[ib2];
				va = vagrid[i2*nzgrid+ib1];
				vab = vagrid[i2*nzgrid+ib2];
				gb = dvdzgrid[i2*nzgrid+ib2];

				vbh[ii] = (zb*vab-z*va)/(zb-z) - 0.5*gb*(zb-z);


				va = vagrid[i2*nzgrid+i1];
				g = dvdzgrid[i2*nzgrid+i1];


				vah[ii] = va;
				dvdzh[ii] = g;
				xvh[ii] = x;
				zvh[ii] = z;
				nvh[ivs[isort[iz]]] += 1;
				if(iz==0) {
					vth[ii] = va;
				} else if (iz==jh-1) {
					vbh[ii] = vth[ii];
				}

				/*
				fprintf(stderr,
			"vt=%6.0f vb=%6.0f va=%6.0f g=%6.2f x=%6.0f z=%6.0f \n",
					vth[ii],vbh[ii],vah[ii],dvdzh[ii],x,z);
				*/
			}
		}


		/* fing vt, bt, dvdz and va at pick positions */
		for(ih=0;ih<nhs;ih++) {
                       	vtopyes = 'n';
			for(jh=0;jh<maxhs;jh++) {
				if(hnums[ih]==vtopupdate[jh]) {
                        		vtopyes = 'y';
					break;
				}
			}
                       	vbotyes = 'n';
			for(jh=0;jh<maxhs;jh++) {
				if(hnums[ih]==vbotupdate[jh]) {
                        		vbotyes = 'y';
					break;
				}
			}

                        if(vtopyes=='n' && vbotyes=='n') continue;

			n2 = ih*maxpks;
			n1 = ih*ncdp_velo;
			nvs = nvh[ih];
			if(nvs==0) {
				fprintf(stderr, 
				"warning: no velan at horizon %d !\n",
				hnums[ih]); 
			} else if (nvs==1) {
				for(ip=0;ip<npicks[ih];ip++) {
					if(vtopyes=='y') veltops[ip+n2]=vth[n1];
                                	if(vbotyes=='y') velbots[ip+n2]=vbh[n1];
                                	velavgs[ip+n2] = vah[n1];
                                	dvdzs[ip+n2] = dvdzh[n1];
				}
			} else {
				for(ip=0;ip<npicks[ih];ip++) {
					x = xpicks[ip+ih*maxpks];
					bisear_(&nvs,&one,xvh+n1,&x,&i2);
					if(x<=xvh[n1]) {
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1];
						if(vbotyes=='y')
                                		velbots[ip+n2] = vbh[n1];
                                		velavgs[ip+n2] = vah[n1];
                                		dvdzs[ip+n2] = dvdzh[n1];
					} else if(x>=xvh[n1+nvs-1]) {
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1+nvs-1];
						if(vbotyes=='y')
                                		velbots[ip+n2] = vbh[n1+nvs-1];
                                		velavgs[ip+n2] = vah[n1+nvs-1];
                                		dvdzs[ip+n2] = dvdzh[n1+nvs-1];
					} else { 
						tmp = (x-xvh[n1+i2-1])/
						      (xvh[n1+i2]-xvh[n1+i2-1]);
						if(vtopyes=='y') 
						veltops[ip+n2] = vth[n1+i2-1] +
						  tmp*(vth[n1+i2]-vth[n1+i2-1]);
						if(vbotyes=='y')
						velbots[ip+n2] = vbh[n1+i2-1] +
						  tmp*(vbh[n1+i2]-vbh[n1+i2-1]);
                                		velavgs[ip+n2] = vah[n1+i2-1] +
						  tmp*(vah[n1+i2]-vah[n1+i2-1]);
                                		dvdzs[ip+n2] = dvdzh[n1+i2-1] +
						  tmp*
						  (dvdzh[n1+i2]-dvdzh[n1+i2-1]);
					}
				}
			}
		}

		free(indx);
		free(xnew);
		free(znew);
		free(vtnew);
		free(vbnew);
		free(vanew);
		free(gnew);

		free(zvs);
		free(ivs);
		free(isort);
		free(zsort);
		free(nvh);
		free(xvh);
		free(zvh);
		free(vbh);
		free(vth);
		free(vah);
		free(dvdzh);
		free(zgrid);
		free(vagrid);
		free(dvdzgrid);
		free(g1);
		free(g2);

		free(z2);
	}
}