#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#define T_Number 400
#define D_Number 601
#define L_Number 601
#define PI 3.141593
void zero(float **buff)
{
int i,j;
for(i=0;i<L_Number;i++)
for(j=0;j<D_Number;j++)
buff[i][j]=0.0;
}
/*float **atten(float **value,int ixs,int izs,float dz,float dx)
{
int i,j;
for(i=0;i<L_Number;i++)
for(j=0;j<D_Number;j++)
value[i][j]=exp(-0.0001*(pow((i-ixs)*dx,2)+pow((j-izs)*dz,2)));
//printf("%lf",value[290][290]);
return value;
} */
void spatial_derivative_xx(float **Dx,float **U)
{
float a0,a1,a2,a3,a4,a5;
int i,j;
a0=-2.92722;
a1=1.66667;
a2=-0.238095;
a3=0.0396825;
a4=-0.00496032;
a5=0.000317460;
for(i=5;i<L_Number-5;i++)
{
for(j=5;j<D_Number-5;j++)
{
Dx[i][j]=a0*U[i][j]+a1*(U[i+1][j]+U[i-1][j])+
a2*(U[i+2][j]+U[i-2][j])+a3*(U[i+3][j]
+U[i-3][j])+a4*(U[i+4][j]+U[i-4][j])
+a5*(U[i+5][j]+U[i-5][j]);
}
}
//printf("%f\n",Dx[299][300]);
}
void spatial_derivative_zz(float **Dz,float **U)
{
float a0,a1,a2,a3,a4,a5;
int i,j;
a0=-2.92722;
a1=1.66667;
a2=-0.238095;
a3=0.0396825;
a4=-0.00496032;
a5=0.000317460;
for(i=5;i<L_Number-5;i++)
for(j=5;j<D_Number-5;j++)
Dz[i][j]=a0*U[i][j]+a1*(U[i][j+1]+U[i][j-1])
+a2*(U[i][j+2]+U[i][j-2])+a3*(U[i][j+3]
+U[i][j-3])+a4*(U[i][j+4]+U[i][j-4])
+a5*(U[i][j+5]+U[i][j-5]);
// printf("%f",Dz[299][300]);
}
float **alloc2(int a,int b)
{
float **v;
int i;
v=(float **)malloc(a*sizeof(float *));
for(i=0;i<a;i++)
v[i]=(float *)malloc(b*sizeof(float *));
return v;
}
void main()
{
float dt,dx,dz,depth,length,time,f0,vp;
int i,j,k;
int izs,ixs;
float f[T_Number];
float **U,**U0,**U1,**uxx,**uzz,**alphy;
FILE *fp,*fp_temp,*fp_test;
dt=0.001;
dx=5.0;
dz=5.0;
depth=3000.0;
length=3000.0;
time=0.4;
f0=20.0;
vp=3000.0;
izs=(D_Number-1)/2;
ixs=(L_Number-1)/2;
U0 = (float**)alloc2(L_Number,D_Number);
U1 = (float**)alloc2(L_Number,D_Number);
U = (float**)alloc2(L_Number,D_Number);
alphy =(float**)alloc2(L_Number,D_Number);
uxx = (float**)alloc2(L_Number,D_Number);
uzz = (float**)alloc2(L_Number,D_Number);
zero(uxx);
zero(uzz);
zero(U0);
zero(U1);
zero(U);
zero(alphy);
alphy[ixs][izs]=1.0;
//alphy=atten(alphy,ixs,izs,dz,dx);
fp=fopen("wavelet.xls","w");
for (k=0;k<T_Number;k++)
{
f[k]=(1-2*PI*f0*(dt*k-1.0/f0)*PI*f0*(dt*k-1.0/f0))*exp(-PI*f0*(dt*k-1.0/f0)*PI*f0*(dt*k-1.0/f0));
fprintf(fp,"%f\t%f\n",k*dt,f[k]);
}
fclose(fp);
for (k=0;k<T_Number;k++)
{
if(k%100==0)
{printf("k=%d\n",k);}
spatial_derivative_xx( uxx, U);
spatial_derivative_zz( uzz, U);
for(i=0;i<L_Number;i++)
{ for(j=0;j<D_Number;j++)
{
U[i][j]=2.0*U1[i][j]-U0[i][j]+(1.0/2)*pow((vp*dt)/dx,2)*uxx[i][j]
+(1.0/2)*pow((vp*dt)/dz,2)*uzz[i][j]+f[k]*alphy[i][j];
}
}
printf("%f\n",U[300][300]);
for(i=0;i<L_Number;i++)
{
for(j=0;j<D_Number;j++)
{
float temp;
temp=U1[i][j];
U1[i][j]=U[i][j];
U0[i][j]=temp;
}
}
}
/* fp_temp=fopen("uxx.xls","w");
for(i=0;i<L_Number;i++)
{
for(j=0;j<D_Number;j++)
{
fprintf(fp_temp,"%f\n",uxx[i][j]);
}
}
fclose(fp_temp);
printf("%f",U[301][301]);*/
fp_temp=fopen("result1.bin","wb");
for(i=0;i<L_Number;i++)
{
for(j=0;j<D_Number;j++){
fwrite(&U[i][j],sizeof(float),1,fp_temp);}
}
fclose(fp_temp);
}
评论0