suinvzco3d.c

su 的源代码库

        {
                temp1= sqrt((ix*dxo)*(ix*dxo)+(iy*dyo)*(iy*dyo))/dxo;
                ix2=temp1;
                temp1=temp1-ix2;

                t3d[iy][ix][iz]= (1.0-temp1)*t2d[ix2][iz]+
                        temp1*t2d[ix2+1][iz];

		p3d[iy][ix][iz]= (1.0-temp1)*p2d[ix2][iz]+
                        temp1*p2d[ix2+1][iz];

		sigma3d[iy][ix][iz]= (1.0-temp1)*sigma2d[ix2][iz]+
                        temp1*sigma2d[ix2+1][iz];

		a23d[iy][ix][iz]= (1.0-temp1)*a2d[ix2][iz]+
                        temp1*a2d[ix2+1][iz];

		dzda23d[iy][ix][iz]= (1.0-temp1)*dzda2d[ix2][iz]+
                        temp1*dzda2d[ix2+1][iz];

	}
	
}


void SingleRay(int nx,int nz,float a, float *v,float dx,float dz,float *x,float *p3,float *sigma,float *t,float *dzda2,int ray_end_iz)
/*****************************************************************************
Ray equation solution of one ray in v(z) medium in polar coordinate
******************************************************************************
Input:
nx              number of x samples
nz              number of z smaples
a		polar angle
v		array for v(z)               
dx              x sample interval
dz              z sample interval

Output:
x		array for ray position of x 
p3		array for slowness p3
sigma		array for sigma
t		array for traveltime
dzda2		array for the z derivative to propagation angle a2
******************************************************************************
Notes: Used internally in function: amplitude()
******************************************************************************
Author:  Meng Xu, Center for Wave Phenomena, 07/8/95
******************************************************************************/
{
	int i,ii;
	float fx,zx,tt,sig,temp,step;

	sig=tt=temp=0;
	fx=0;
	zx=dx*nx;

	for(i=0;i<nz;++i)
	{
	sigma[i]=0;
	p3[i]=0;
	x[i]=0;
	t[i]=0;
	dzda2[i]=0;
	}

	
	for(i=0;i<nz;++i)
	{	
		if(1./v[i]/v[i]-sin(a)*sin(a)/v[0]/v[0]>0)  
			{
			p3[i]=sqrt(ABS(1./v[i]/v[i]-sin(a)*sin(a)/v[0]/v[0]));

			sigma[i]=sig;
			sig+=dz/p3[i];

			x[i]=sigma[i]*sin(a)/v[0];
			t[i]=tt;
			tt+=dz/v[i]/v[i]/p3[i];
			dzda2[i]=temp;
			temp-=dz/p3[i]*sin(a)/v[0]*p3[0]/p3[i];
			ray_end_iz=nz-1;
			if(p3[i]<=TINY) 
				{
				fprintf(stderr,"\nTurning point: z=%f\n",i*dz);
				continue;
				}
			if(x[i]<fx||x[i]>zx) {step=i;}
			}
		else
			{
			fprintf(stderr,"\nTurning point: ang=%f,z=%f",a,i*dz);
			ray_end_iz=i;

			for(ii=i; ii<nz; ++ii) 
				{
				sigma[i]=0;
       		                x[i]=-1;
       		                t[i]=0;
                       		dzda2[i]=0;
				}
                        if(x[i]<fx||x[i]>zx) {step=i;}
			i=nz;
			}
	}
}
		


	
void para2d(int na,int nx,int nz,float dx,float dz,float *v,float **t2d,
		float **sigma2d,float **p2d,float **dzda2d,float **a2d)
/*****************************************************************************
Calculate ray parameters on a 2D grid
******************************************************************************
Input:
nx             number of output grids in x
nz             number of output grids in z
dx             output sample rate in x
dz             output sample rate in z
v	       velocity array v(z)
	
Output:
t2d             2D array for traveltime
sigma2d         2D array for sigma
p2d             2D array for p3
dzda2d          2D array for dz/da
a2d             2D array for a2

******************************************************************************
Notes: Used internally in function: amplitude()
******************************************************************************
Author:  Meng Xu, Center for Wave Phenomena, 07/8/95
******************************************************************************/
{
	int i,k,ii;
	int iflag;
	int i00=0;
	int imax;
	int *ray_end_iz;

	float sx;
	float x0;
	float da;	
	float a,**p3,**sigma,**t,**x,**dzda2,**a2;

	da=3.14159265/na/2.;

	p3=ealloc2float(nz,na);
	sigma=ealloc2float(nz,na);
	t=ealloc2float(nz,na);
	x=ealloc2float(nz,na);
	dzda2=ealloc2float(nz,na);
	a2=ealloc2float(nz,na);
	ray_end_iz=ealloc1int(na);
	

	for(i=0,a=0.;i<na;++i,a+=da)
	{
	SingleRay(nx,nz,a,v,dx,dz,x[i],p3[i],
		sigma[i],t[i],dzda2[i],ray_end_iz[i]);
	for(k=0;k<nz;++k) a2[i][k]=a;
	}

	/*interpolation from ray coordinate to (x,z) coordinate*/ 
	for(k=1;k<nz;++k)
	  for(i=0;i<nx;++i)
		{
		x0=dx*i;
		iflag=0;
		for(ii=0;ii<na;++ii)
			{
			if(ii==0) {i00=0;}
			if(x[ii][k]>x0&&x[i00][k]<=x0)
				{
				iflag=ii; 
				ii=na;
				}
			}
		if(iflag==0) 
			{
			for(ii=0;ii<na;++ii)
                        	{
                        	if(ray_end_iz[ii]==k) {imax=ii;}
				}
			iflag=na-1;
			sx=0;
			}
		else	{
			sx=(x[iflag][k]-x0)/(x[iflag][k]-x[iflag-1][k]);
			}	
			
		t2d[i][k]=(1.0-sx)*t[iflag][k] + 
					sx*t[iflag-1][k];
		sigma2d[i][k]=(1.0-sx)*sigma[iflag][k] +     
                                        sx*sigma[iflag-1][k];
		p2d[i][k]=(1.0-sx)*p3[iflag][k] +     
                                        sx*p3[iflag-1][k];
		a2d[i][k]=(1.0-sx)*a2[iflag][k] +
                                        sx*a2[iflag-1][k];
		dzda2d[i][k]=(1.0-sx)*dzda2[iflag][k] +
                                        sx*dzda2[iflag-1][k];
if(t2d[i][k]==0) fprintf(stderr,"z=%f\n",k*dz);
		}	
	free2float(p3);
	free2float(sigma);
	free2float(t);
	free2float(x);
	free2float(dzda2);
	free2float(a2);	

}

void amplitude(float fx, float fy, int nx, int ny, int nz, 
	float dx, float dy, float dz,
    	float fxo, float fyo, int nxo, int nyo, int nzo, 
	float dxo, float dyo, float dzo,
    	int na, float *v, float off, float ***t3d,float ***amp, float ***p3d)
/*****************************************************************************
Calculate Jacobian and Beylkin determinant to give the amplitude weight
******************************************************************************
Input:
fx		first x position in input data
fy		first y position in input data
nx=nxo		number of velocity grids in x
ny=nyo		number of velocity grids in y
nz		number of velocity grids in z
fxo		first x position in output
fyo		first y position in output		
nxo             number of output grids in x
nyo             number of output grids in y
nzo             number of output grids in z
dxo             output sample rate in x
dyo             output sample rate in y
dzo             output sample rate in z
off		offset
na		number of ray angles
v		1D array of velocity v(z)

Output:
t3d             3D array for traveltime
amp         3D array for amplitude weight
p3d		slowness  in z direction p_3
******************************************************************************
Author:  Meng Xu, Center for Wave Phenomena, 07/8/95
******************************************************************************/
{
        int i,j,k;
	int nxr;
	int ix,jl,jr,il,ir;

	float xi,sx;
	float r,a1,a2,a1s,a1g,a2s,a2g,temp,temp1,temp2;
	float det[3][3],ps[3],pg[3];
        float **t2d,**p2d,**sigma2d,**dzda2d,**a2d;
	float ***jacob,***asag,***h;
	float ***sigma3d,***dzda23d,***a23d;

	nxr = NINT(sqrt(nxo*dxo*nxo*dxo+nyo*dyo*nyo*dyo)/dxo)+1;

        p2d=ealloc2float(nzo,nxr);
        t2d=ealloc2float(nzo,nxr);
        sigma2d=ealloc2float(nzo,nxr);
        dzda2d=ealloc2float(nzo,nxr);
        a2d=ealloc2float(nzo,nxr);

	jacob=ealloc3float(nzo,nxo,nyo);
	asag=ealloc3float(nzo,nxo,nyo);
	h=ealloc3float(nzo,nxo,nyo);
        sigma3d=ealloc3float(nzo,nxo,nyo);
        dzda23d=ealloc3float(nzo,nxo,nyo);
        a23d=ealloc3float(nzo,nxo,nyo);


	para2d(na,nxr,nzo,dxo,dzo,v,t2d,
                sigma2d,p2d,dzda2d,a2d);

/*	for(i=0;i<nx;++i)
        for(k=0;k<nz;++k)
	{
        fprintf(stderr,"%f\n",t2d[i][k]);	
	r=sqrt(i*dx*i*dx+k*dz*k*dz);
	fprintf(stderr,"%f\n",r/v[0]);
	}
*/
	para3d(
                nxo,nyo,nzo,
                dxo,dyo,dzo,
                t3d,sigma3d,p3d,
                dzda23d,a23d,
                t2d,sigma2d,p2d,dzda2d,a2d);

	/*compute J and thus A*/
	xi=off/dxo;
	ix=xi;
	sx=xi-ix;
	
	for(i=0;i<nyo;++i)
	 for(j=0;j<nxo;++j)
	  for(k=1;k<nzo;++k) {
		if(j==0)
			a1=3.14159265/2.;
		else 
		    a1=atan(i*dyo/j*dxo);

		a2=a23d[i][j][k];
		det[0][0]=sin(a2)*cos(a1)/v[0];
		det[0][1]=sin(a1)*sin(a2)/v[0];
		det[0][2]=p3d[i][j][k];
		det[1][0]=-sigma3d[i][j][k]*sin(a1)*sin(a2)/v[0];
		det[1][1]=sigma3d[i][j][k]*cos(a1)*sin(a2)/v[0];
		det[1][2]=0.0;
		det[2][0]=sigma3d[i][j][k]*cos(a1)*cos(a2)/v[0];
		det[2][1]=sigma3d[i][j][k]*sin(a1)*cos(a2)/v[0];
		det[2][2]=dzda23d[i][j][k];
		temp=4.*PI*sqrt(v[0]*ABS(determ(det)));
		if(temp>0.0){
			jacob[i][j][k]=100.*sqrt(sin(a2))/temp;
		}	
	}
	for(i=0;i<nyo;++i)
         for(j=0;j<nxo;++j)
          for(k=0;k<nzo;++k)
		asag[i][j][k]=jacob[i][j][k]*jacob[i][ABS(j-ix)][k];

	/*compute h*/
	for(i=1;i<nyo;++i)
         for(j=0;j<nxo;++j)
          for(k=1;k<nzo;++k)
		{
		if(t3d[i][j][k]==0)
                        {
                        amp[i][j][k]=0;
                        }
                else{
		if(j==0) a1s=PI/2.;
		else a1s=atan(i*dyo/j/dxo);
		if(ABS(j*dxo-off)==0) a1g=PI/2;
		else a1g=atan(i*dyo/ABS(j*dxo-off));
		a2s=a23d[i][j][k];
		a2g=a23d[i][ABS(j-ix)][k];
		if(j==nxo-1){jl=j-1; jr=j;}
		else {jl=j-1; jr=j+1;}
		temp1=(sigma3d[i][jl][k]-sigma3d[i][jr][k])/dxo/2.;	
		if(ABS(j-ix)==0){jl=ABS(j-ix); jr=jl+1;} 
                else {jl=ABS(j-ix)-1; jr=jl+1;}    
                temp2=(sigma3d[i][jl][k]-sigma3d[i][jr][k])/dxo;
		
		ps[0]=1./v[0]*cos(a1s)*sin(a2s);
		ps[1]=1./v[0]*sin(a1s)*sin(a2s);
		ps[2]=p3d[i][j][k];
		pg[0]=1./v[0]*cos(a1g)*sin(a2g);
		pg[1]=1./v[0]*sin(a1g)*sin(a2g);
                pg[2]=p3d[i][ABS(j-ix)][k];

		det[0][0]=ps[0]+pg[0];
		det[0][1]=ps[1]+pg[1];
		det[0][2]=ps[2]+pg[2];

		det[1][0]=-(1.+ps[0]*temp1)/sigma3d[i][j][k]
			-(1.+pg[0]*temp2)/sigma3d[i][ABS(j-ix)][k];
		det[1][1]=-(ps[1]*temp1)/sigma3d[i][j][k]
			-(pg[1]*temp2)/sigma3d[i][ABS(j-ix)][k];
		det[1][2]=(ps[0]*(1.+ps[0]*temp1)/sigma3d[i][j][k]
			+ps[1]*(ps[1]*temp1)/sigma3d[i][j][k])/ps[2]
			+(pg[0]*(1.+pg[0]*temp2)/sigma3d[i][ABS(j-ix)][k]
                        +pg[1]*(pg[1]*temp2)/sigma3d[i][ABS(j-ix)][k])/pg[2];
		if(i==nyo-1){il=i-1; ir=i;}
                else {il=i-1; ir=i+1;}
                temp1=(sigma3d[il][j][k]-sigma3d[ir][j][k])/dyo/2.;
                temp2=(sigma3d[il][ABS(j-ix)][k]-sigma3d[ir][ABS(j-ix)][k])
			/2./dyo;

		det[2][0]=-(ps[0]*temp1)/sigma3d[i][j][k]
                        -(pg[0]*temp2)/sigma3d[i][ABS(j-ix)][k];
                det[2][1]=-(1.+ps[1]*temp1)/sigma3d[i][j][k]
                        -(1.+pg[1]*temp2)/sigma3d[i][ABS(j-ix)][k];
                det[2][2]=(ps[0]*(ps[0]*temp1)/sigma3d[i][j][k]
                        +ps[1]*(1.+ps[1]*temp1)/sigma3d[i][j][k])/ps[2]
                        +(pg[0]*(pg[0]*temp2)/sigma3d[i][ABS(j-ix)][k]
                        +pg[1]*(1.+pg[1]*temp2)/sigma3d[i][ABS(j-ix)][k])/pg[2];
		h[i][j][k]=determ(det);

		r=sqrt((ps[0]+pg[0])*(ps[0]+pg[0])+(ps[1]+pg[1])*(ps[1]+pg[1])
		   +(ps[2]+pg[2])*(ps[2]+pg[2]));
		amp[i][j][k]=h[i][j][k]/asag[i][j][k]/r;
		}

	}
	 for(i=1;i<nxo;++i)
          for(k=1;k<nzo;++k)
		amp[0][i][k]=amp[1][i][k]; 

	free2float(p2d);
        free2float(t2d);
        free2float(sigma2d);
        free2float(dzda2d);
        free2float(a2d);

}