velpertan.c

su 的源代码库

		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw); 
		td[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp2[icdp],zdcdp2[icdp],off,VP0,eps,del+0.05,
		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);

	 	tkx[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp1[icdp],zdcdp1[icdp],off,VP0,eps,del,kz,
		kx+0.05,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);
		tkx[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp2[icdp],zdcdp2[icdp],off,VP0,eps,del,kz,
		kx+0.05,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);

	 	tkz[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp1[icdp],zdcdp1[icdp],off,VP0,eps,del,
		kz+0.05,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,
		vp,e,d,r1,r2,win,rw);
		tkz[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp2[icdp],zdcdp2[icdp],off,VP0,eps,del,
		kz+0.05,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,
		vp,e,d,r1,r2,win,rw);

	 	tde[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp1[icdp],zdcdp1[icdp],off,VP0,eps+0.05,del,
		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);
		tde[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp2[icdp],zdcdp2[icdp],off,VP0,eps+0.05,del,
		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);
/*
	 	tdV[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp1[icdp],zdcdp1[icdp],off,VP0+10,eps,del,
		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);
		tdV[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp2[icdp],zdcdp2[icdp],off,VP0+10,eps,del,
		kz,kx,nx,nz,fx,dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,
		r1,r2,win,rw);
*/

			}
	}

	/*compute average of depths and their derivatives*/
	for(icdp=0;icdp<ncdp;icdp++){
	sum1=sum2=sum3=sum4=sum5=sum6=sum7=sum8=sum9=sum10=sum11=sum12=0.0;
		for(ioff=0;ioff<noff;ioff++){
/*
		dt1[0][icdp][ioff]=-pp[0][icdp][ioff]*
			((tdV[0][icdp][ioff]-t[0][icdp][ioff])/10);
		sum1 += dt1[0][icdp][ioff]/noff;
		dt2[0][icdp][ioff]=-pp[1][icdp][ioff]*
			((tdV[1][icdp][ioff]-t[1][icdp][ioff])/10);
		sum2 += dt2[0][icdp][ioff]/noff;
*/

		dt1[1][icdp][ioff]=-pp[0][icdp][ioff]*((tkx[0][icdp][ioff]-
			t[0][icdp][ioff])/0.05);
		sum3 += dt1[1][icdp][ioff]/noff;
		dt2[1][icdp][ioff]=-pp[1][icdp][ioff]*((tkx[1][icdp][ioff]-
			t[1][icdp][ioff])/0.05);
		sum4 += dt2[1][icdp][ioff]/noff;

		dt1[2][icdp][ioff]=-pp[0][icdp][ioff]*((tde[0][icdp][ioff]-
			t[0][icdp][ioff])/0.05);
		sum5 += dt1[2][icdp][ioff]/noff;
		dt2[2][icdp][ioff]=-pp[1][icdp][ioff]*((tde[1][icdp][ioff]-
			t[1][icdp][ioff])/0.05);
		sum6 += dt2[2][icdp][ioff]/noff;

		dt1[3][icdp][ioff]=-pp[0][icdp][ioff]*((td[0][icdp][ioff]-
			t[0][icdp][ioff])/0.05);
		sum7 += dt1[3][icdp][ioff]/noff;
		dt2[3][icdp][ioff]=-pp[1][icdp][ioff]*((td[1][icdp][ioff]-
			t[1][icdp][ioff])/0.05);
		sum8 += dt2[3][icdp][ioff]/noff;

		dt1[4][icdp][ioff]=-pp[0][icdp][ioff]*((tkz[0][icdp][ioff]-
			t[0][icdp][ioff])/0.05);
		sum9 += dt1[4][icdp][ioff]/noff;
		dt2[4][icdp][ioff]=-pp[1][icdp][ioff]*((tkz[1][icdp][ioff]-
			t[1][icdp][ioff])/0.05);
		sum10 += dt2[4][icdp][ioff]/noff;

		off=off0+ioff*doff;
		zz1[icdp][ioff]=sqrt(pow(zcdp1[icdp],2)
				+r11[icdp]*pow((off/2),2)+r12[icdp]*
				pow((off/2),4)/(pow((off/2),2)+pow(zcdp1[icdp],2)));

		sum11 += zz1[icdp][ioff]/noff;
		zz2[icdp][ioff]=sqrt(pow(zcdp2[icdp],2)
				+r21[icdp]*pow((off/2),2)+r22[icdp]*
				pow((off/2),4)/(pow((off/2),2)+pow(zcdp2[icdp],2)));
		sum12 += zz2[icdp][ioff]/noff;
		}
		dt1avg[0][icdp]=sum1;dt2avg[0][icdp]=sum2;
		dt1avg[1][icdp]=sum3;dt2avg[1][icdp]=sum4;
		dt1avg[2][icdp]=sum5;dt2avg[2][icdp]=sum6;
		dt1avg[3][icdp]=sum7;dt2avg[3][icdp]=sum8;
		dt1avg[4][icdp]=sum9;dt2avg[4][icdp]=sum10;
		zz1avg[icdp]=sum11;zz2avg[icdp]=sum12;	
	}

	for(ipar1=0;ipar1<5;ipar1++)
		for(ipar2=0;ipar2<5;ipar2++)
			for(icdp=0;icdp<ncdp;icdp++)
			for(ioff=0;ioff<noff;ioff++){
				A1[ipar1][ipar2] +=
				(dt1[ipar1][icdp][ioff]-dt1avg[ipar1][icdp])
				*(dt1[ipar2][icdp][ioff]-dt1avg[ipar2][icdp]);
			}
	for(ipar1=0;ipar1<5;ipar1++)
		for(ipar2=0;ipar2<5;ipar2++)
			for(icdp=0;icdp<ncdp;icdp++)
			for(ioff=0;ioff<noff;ioff++){
				A2[ipar1][ipar2] +=
				(dt2[ipar1][icdp][ioff]-dt2avg[ipar1][icdp])
				*(dt2[ipar2][icdp][ioff]-dt2avg[ipar2][icdp]);
			}
	
	for(ipar1=0;ipar1<5;ipar1++)
		for(icdp=0;icdp<ncdp;icdp++)
		for(ioff=0;ioff<noff;ioff++){
			b1[ipar1] +=
			(dt1[ipar1][icdp][ioff]-dt1avg[ipar1][icdp])
			*(zz1[icdp][ioff]-zz1avg[icdp]);
		}
		
	for(ipar1=0;ipar1<5;ipar1++)
		for(icdp=0;icdp<ncdp;icdp++)
		for(ioff=0;ioff<noff;ioff++){
			b2[ipar1] +=
			(dt2[ipar1][icdp][ioff]-dt2avg[ipar1][icdp])
			*(zz2[icdp][ioff]-zz2avg[icdp]);
		}
	
	/*compute A and b*/
	printf("Matrix A\n");
	for(ipar1=0;ipar1<5;ipar1++){
		for(ipar2=0;ipar2<5;ipar2++){
			A[ipar1][ipar2]=A1[ipar1][ipar2]+A2[ipar1][ipar2];
			printf("%f \t",A[ipar1][ipar2]);
		}
	printf("\n \n");
	b[ipar1]=-b1[ipar1]-b2[ipar1];
	}

	printf("Vector b\n");
	printf("%f %f %f %f %f\n",b[0],b[1],b[2],b[3],b[4]);
	conjugate_gradient(A,x,b,5,0.0001);

	printf("\n \n");
		printf("d_VP0=%e d_kz=%e d_kx=%e d_epsilon=%e d_delta=%e\n",x[0],x[4],x[1],x[2],x[3]);
	printf("\n \n");
	/* deallocate space */
        free1float(xcip);
        free1float(r11cip);
        free1float(r12cip);
        free1float(r21cip);
        free1float(r22cip);
        free1float(xcdp);
        free1float(zcdp1);
        free1float(zcdp2);
        free1float(zdcdp1);
        free1float(zdcdp2);
        free1float(r12);
        free1float(r22);
        free1float(r21);
        free1float(r11);
        free1float(zint1);
        free1float(xint1);
        free1float(zint2);
        free1float(xint2);
	free3float(t);
	free3float(tkz);
	free3float(tkx);
	free3float(td);

	return EXIT_SUCCESS;
}

/******************************************************************************/
/*****  interpolate depths and compute slope of the reflector at each cdp *****/
/******************************************************************************/
void interpolate (int ncdp, float *xcdp, float *xint, float *zint, float *zcdp, float *zdcdp, int np, char *method)
{
float (*zind)[4];

                /* else, if monotonic interpolation */
                if (method[0]=='m') {
                                zind = (float (*)[4])ealloc1float(np*4);
                                cmonot(np,xint,zint,zind);
                                intcub(0,np,xint,zind,ncdp,xcdp,zcdp);
                                intcub(1,np,xint,zind,ncdp,xcdp,zdcdp);

                /* else, if Akima interpolation */
                } else if (method[0]=='a') {
                                zind = (float (*)[4])ealloc1float(np*4);
                                cakima(np,xint,zint,zind);
                                intcub(0,np,xint,zind,ncdp,xcdp,zcdp);
                                intcub(1,np,xint,zind,ncdp,xcdp,zdcdp);

                /* else, if cubic spline interpolation */
                } else if (method[0]=='s') {
                                zind = (float (*)[4])ealloc1float(np*4);
                                csplin(np,xint,zint,zind);
                                intcub(0,np,xint,zind,ncdp,xcdp,zcdp);
                                intcub(1,np,xint,zind,ncdp,xcdp,zdcdp);
				}
}

/*****************************************************************************/
/*************** Compute time for a given offset and cdp point ***************/
/*****************************************************************************/
float velpertan_time(float x00,float z00,float *xxbound,float *zzbound,float x,float z,
	float zd,float off,float VP0,float eps,float del,float kz,float kx,
	int nx,int nz,float fx,float dx,float fz,float dz,int nt,float dt,
	float amin, float amax,float tol,float *p,float *vp,float *e,float *d,
	float r1,float r2,int *win,float rw){

int iz,ix,NRS;
float a_normal,a_vert=0.0,da,X,da_2,angle,T,t1,t2,c1,c2;
float **a3333, **a1111, **a1133, **a1313, **a1113, **a3313;
float **e_temp, **d_temp;
float X1,X2,p1,p2;
Ray *ray;

a3333=alloc2float(nz,nx);
a1111=alloc2float(nz,nx);
a1133=alloc2float(nz,nx);
a1313=alloc2float(nz,nx);
a1113=alloc2float(nz,nx);
a3313=alloc2float(nz,nx);
e_temp=alloc2float(nz,nx);
d_temp=alloc2float(nz,nx);

a_normal=atan(1/zd);
if (a_normal<0){ a_normal=(a_normal*180/3.1415-90); 
	a_vert=180+a_normal;
	}
if (a_normal>0){ a_normal=(a_normal*180/3.1415+90);
	a_vert=180-a_normal;
	}
if (fabs(a_normal)<=90)
warn("slope of reflector is infeasible; please pick reflector again");

/************************ setup the velocity field ****************************/
                for (ix=0;ix<nx;ix++)
                        for (iz=0;iz<nz;iz++){
			if((iz*dz)<=zzbound[ix]) a3333[ix][iz] = vp[ix*nz+iz];
			else
                        a3333[ix][iz] = VP0+kx*(fx+dx*ix-x00)+kz*(fz+dz*iz-z00);

				a1313[ix][iz] = 0.0;
				a1113[ix][iz] = 0.0;
				a3313[ix][iz] = 0.0;
			}
		smooth2(nz,nx,r1,r2,win,rw,a3333);

                for (ix=0; ix<nx; ++ix)
                        for (iz=0;iz<nz; ++iz)
                                a3333[ix][iz]=pow(a3333[ix][iz],2);

                for (ix=0; ix<nx; ++ix)
                        for (iz=0; iz<nz; ++iz){
			if((iz*dz)<=zzbound[ix]){
			e_temp[ix][iz]=e[ix*nz+iz];
			d_temp[ix][iz]=d[ix*nz+iz];
			}
			else {
			e_temp[ix][iz]=eps;
			d_temp[ix][iz]=del;
			}
		}
		smooth2(nz,nx,r1,r2,win,rw,e_temp);
		smooth2(nz,nx,r1,r2,win,rw,d_temp);

                for (ix=0; ix<nx; ++ix)
                        for (iz=0; iz<nz; ++iz){
			a1111[ix][iz] = a3333[ix][iz]+2*e_temp[ix][iz]* 
							a3333[ix][iz];
			a1133[ix][iz] = sqrt(2*d_temp[ix][iz]*a3333[ix][iz]*
                                    a3333[ix][iz]+a3333[ix][iz]*a3333[ix][iz]);
                }

/************ find the time for a given offset and diffractor point **********/
	/* find the closest two rays that bracket the desired offset */
	X=0; angle=0;da=-10.0;
	do{
	da=da+10.0;
	ray = makeRay(x,z,a_normal+da,nt,dt,a1111,a3333,a1133,a1313,
			a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);          
	NRS=ray->nrs;	
	if(ray->rs[NRS-1].z>0) printf("error in computation\n");
	X1=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
	freeRay(ray);
	ray = makeRay(x,z,a_normal-da,nt,dt,a1111,a3333,a1133,a1313,
				a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);   
	NRS=ray->nrs;	
	X2=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
	freeRay(ray);
	X=X2-X1;
	if (fabs(X-off)<0.003) break;
	}while(X<off);
	da_2=da-10.0;

	if(fabs(X-off)>0.003) 
	angle = zbrentou(da,da_2,tol,off,x,z,a_normal,nt,dt,a1111,a3333,a1133,a1313,a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);

	ray = makeRay(x,z,a_normal+angle,nt,dt,a1111,a3333,a1133,a1313,
				a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);   
	NRS=ray->nrs;	

	X1=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
	t1=-((ray->rs[NRS-1].t-ray->rs[NRS-2].t)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].t;
	p1=1/ray->rs[0].v;
	if(a_normal>0.0) c1=cos((a_vert-angle)*3.1415/180);
		else 	 c1=cos((a_vert+angle)*3.1415/180);
	freeRay(ray);
	ray = makeRay(x,z,a_normal-angle,nt,dt,a1111,a3333,a1133,a1313,
				a1113,a3313,nx,dx,fx,nz,dz,fz,amax,amin);   
	NRS=ray->nrs;	

	X2=-((ray->rs[NRS-1].x-ray->rs[NRS-2].x)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].x;
	t2=-((ray->rs[NRS-1].t-ray->rs[NRS-2].t)/(ray->rs[NRS-1].z-
		ray->rs[NRS-2].z))*ray->rs[NRS-2].z+ray->rs[NRS-2].t;
	p2=1/ray->rs[0].v;
	if(a_normal>0.0) c2=cos((a_vert+angle)*3.1415/180);
		else 	 c2=cos((a_vert-angle)*3.1415/180);
	freeRay(ray);
	T=t1+t2;
	*p=1/(p1*c1+p2*c2);
	printf("est_offset=%f true_offset=%f time=%f cip=%f \n",X2-X1,off,T,x);

	/* free elements of the stiffness tensor */
	free2float(a3333);
	free2float(a1111);
	free2float(a1133);
	free2float(a1313);
	free2float(a1113);
	free2float(a3313);
	free2float(e_temp);
	free2float(d_temp);

	return T;