velo2hfile.c

来自「seismic software,very useful」· C语言 代码 · 共 921 行 · 第 1/2 页

   			fprintf(hfilefp, 
"%9.1f %9.1f %8.0f %8.0f   %1d   %4.2f  %4.0f  %4.2f \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]);
		}
	}
}


/* compute velocities from VELO and DVDZ cards */
/* this if for hfile */
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) {

	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, j;
	float *xnew, *znew, *vanew, *vtnew, *vbnew, *gnew, xpre, xnow;
	int *indx, inow, ihbot;

	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];
			}
		}

		/* compute average v and dvdz grids */
		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));

		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 (velocity above) and vb (velocity below) 
			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 {
						/* 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]);
					}
					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;
                                j1 = tmp;
				if(iz<jh-1) {
					tmp = zvs[iz+1]/dz;
					j2 = tmp;
				} else { 
                                	j2 = j1 + 1;
				}
                                ilimit(&j1,0,nzgrid-1);
                                ilimit(&j2,0,nzgrid-1);
                                ii = ivs[isort[iz]]*ncdp_velo
                                        +nvh[ivs[isort[iz]]];
				if(ztol>0.) {
					n = nps_velo[i2];
					n1 = i2*maxp;
					tmp = fabs(z-z_velo[n1]);
					i1 = 0;	
					for(j=1;j<n;j++) {
						if(fabs(z-z_velo[n1+j])<tmp) {
						    i1 = j;
						    tmp=fabs(z-z_velo[n1+j]); 	
						}
					}
					if(tmp<ztol) {
						va = v_velo[n1+i1];
						if(i1+1<n) {
							vab = v_velo[n1+i1+1];
							zb = z_velo[n1+i1+1];
						} else {
                                			va=vagrid[i2*nzgrid+j1];
                                			z = zgrid[j1];
                                			vab = 
							  vagrid[i2*nzgrid+j2];
                                			zb = zgrid[j2];
						}
					} else {
                                		va = vagrid[i2*nzgrid+j1];
                                		z = zgrid[j1];
                                		vab = vagrid[i2*nzgrid+j2];
                                		zb = zgrid[j2];
					}
				} else {
                                	va = vagrid[i2*nzgrid+j1];
                                	z = zgrid[j1];
                                	vab = vagrid[i2*nzgrid+j2];
                                	zb = zgrid[j2];
				}

                                tmp = z/dz;
                                j1 = tmp;
                                tmp = zb/dz;
                                j2 = tmp;
                                ilimit(&j1,0,nzgrid-1);
                                ilimit(&j2,0,nzgrid-1);
				g = 0.;
				for(j=j1;j<j2+1;j++) {
					g = g + dvdzgrid[i2*nzgrid+j];
				}
				if(j2>j1) g  = g / (j2-j1+1);
			    
                                vth[ii] = (zb*vab-z*va)/(zb-z) - 0.5*g*(zb-z);
                                vbh[ii] = (zb*vab-z*va)/(zb-z) + 0.5*g*(zb-z);
                                vah[ii] = va;
                                dvdzh[ii] = g;
                                xvh[ii] = x;
                                zvh[ii] = zvs[iz];
                                nvh[ivs[isort[iz]]] += 1;

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

			}
		}


		/* fing vt, vb 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;

			/* find velocity at top of a layer */
			if(vtopyes=='y') {
				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 { 
					for(ip=0;ip<npicks[ih];ip++) {
						findv(nvs,xvh+n1,vth+n1,
							xpicks[ip+n2],
							&veltops[ip+n2]);  
					}
				}
			}
			/* now for the bottom velocity of the layer */
			if(vbotyes=='y') {
				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 { 
					for(ip=0;ip<npicks[ih];ip++) {
						findv(nvs,xvh+n1,vbh+n1,
							xpicks[ip+n2],
							&velbots[ip+n2]);  
					}
				}
			}
		}

		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);
	}
}

/* find velocity at x position of an horizon velocity function */ 
void findv(int ns, float *xs, float *vs, float x, float *v) { 

	int one=1, i2;
	float tmp;

	if (ns==1) {
		*v = vs[0];
	} else {
	       	bisear_(&ns,&one,xs,&x,&i2);
		if(x<=xs[0]) { 
			*v = vs[0];
		} else if(x>=xs[ns-1]) {
			*v = vs[ns-1];
		} else { 
			tmp = (x-xs[i2-1])/(xs[i2]-xs[i2-1]);
			*v = vs[i2-1] + tmp*(vs[i2]-vs[i2-1]);
		}
	}
}

/* find an interface index ihbot directly below an x position of an 
	interfaceat of index ihxz */  
/* *ihbot is the index found. If -1, index not found	*/

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

	int ih, n1, n2, one=1;
	float xmin, xmax, z;
	float *zxs, *zsort;
	int *ixs, *ixsave, *isort;
	int i1, i2, nxs, iz;

	zxs = (float*) malloc(nhs*sizeof(float));
	ixs = (int*) malloc(nhs*sizeof(int));
	ixsave = (int*) malloc(nhs*sizeof(int));
	isort = (int*) malloc(nhs*sizeof(int));
	zsort = (float*) malloc(nhs*sizeof(float));

	/* find cross point at horizons */
	nxs = 0;
	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) {
			ixs[nxs] = 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 {
			/* 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]);
			}
			zxs[nxs] = z; 
			nxs = nxs + 1;
		}
	}

	/* sort zxs into ascending order */
	for (iz=0;iz<nxs;iz++) {
        	isort[iz] = iz;
                zsort[iz] = zxs[iz];
        }
        if(nxs>0) qkisort(nxs,zsort,isort);
        for(iz=0;iz<nxs;iz++) zxs[iz] = zsort[isort[iz]];
        for(iz=0;iz<nxs;iz++) ixsave[iz] = ixs[iz];
        for(iz=0;iz<nxs;iz++) ixs[iz] = ixsave[isort[iz]];

	*ihbot = -1;

	for(iz=0;iz<nxs;iz++) {
		if(ihxz==ixs[iz]) {
			if(iz+1 < nxs)
				*ihbot = ixs[iz+1];
		break;
		}
	}

	free(ixs);
	free(ixsave);
	free(zxs);
	free(zsort);
	free(isort);
}