声波有限差分法正演模拟c语言程序

#include

#include

#define fm 30

#define dt 0.001

#define PI 3.1415926

#define Nt 401

#define Nx 200

#define Nz 200

//---------------加载震源，雷克子波-----------------------------

void fun(float source[])

{

FILE *fp;

int it,i;

float t1,t2,t0;

for(i=0;i

source[i]=0.0;

t0 = 1.0/fm;

printf("t0 = %f\n", t0);

fp=fopen("ricker.txt","w");

for(it=0;it

{

t1=it*dt;

t2=t1-t0;

source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));

fprintf(fp,"%.8f %.8f\n",t1,source[it]);

}

fclose(fp);

}

//--------------------主程序，差分算子----------------------------------------

void main()

{ FILE *fp, *fpwf;

int it,ix,iz, order,N;

int nsx,nsz;

float a0,a1,a2,a3,a4,a5,a6;

int flag;

float dx,dz;

float Un[Nx][Nz], Unm[Nx][Nz], Unp[Nx][Nz], v[Nx][Nz],record[Nx][Nt];//Un表示n时刻，Unm表示n-1时刻，Unp表示n+1时刻波场

float source[Nt];

nsx=Nx/2;

nsz=0.0;//震源位置

N=6;

order=2*N;//精度

if(order == 2)

{

a0=-2.00000;a1=1.00000;a2=0;a3=0;a4=0;a5=0;a6=0;

}

if(order == 4)

{

a0=-2.50000;a1=1.33333;a2=-8.33333e-2;a3=0;a4=0;a5=0;a6=0;

}

if(order == 6)

{

a0=-2.72222;a1=1.50000;a2=-1.50000e-1;a3=1.11111e-2;a4=0;a5=0;a6=0;

}

if(order == 8)

{

a0=-2.84722;a1=1.60000;a2=-2.00000e-1;a3=2.53968e-2;a4=-1.78571e-3;a5=0;a6=0;

}

if(order == 10)

{

a0=-2.92722;a1=1.66667;a2=-2.38095e-1;a3=3.96825e-2;a4=-4.96032e-3;a5=3.17460e-4;a6=0;

}

if(order == 12)

{

a0=-2.98278;a1=1.71429;a2=-2.67857e-1;a3=5.29101e-2;a4=-8.92857e-3;a5=1.03896e-3;a6=-6.01251e-5;

}

/*for (ix=0;ix

{ for (iz=0;iz<(Nz*2/5);iz++)

v[ix][iz]=2000;

for (iz=(Nz*2/5);iz<(Nz*3/5);iz++)

v[ix][iz]=3000;

for (iz=(Nz*3/5);iz

v[ix][iz]=4000;//给速度赋值

}

*/

for (ix=0;ix

{ for (iz=0;iz

v[ix][iz]=4000;

for (iz=Nz/2;iz

v[ix][iz]=4000;

}

dx=10.0;

dz=10.0;//赋值

for(ix=0;ix

{

for(iz=0;iz

{

Un[ix][iz]=0.0;

Unm[ix][iz]=0.0;

Unp[ix][iz]=0.0;

}

}

fun(source);

fpwf = fopen("Wavefieldall.dat", "wb");

for(it=0;it

{

for(ix=6;ix<(Nx-6);ix++)

for(iz=6;iz<(Nz-6);iz++)

{

Unp[ix][iz]=2*Un[ix][iz]-

Unm[ix][iz]+0.5*v[ix][iz]*v[ix][iz]*dt*dt/(dx*dx)*(a0*Un[ix][iz]+a1*(Un[ix+1][iz]+Un[ix-1][iz])+a2*(Un[ix+2][iz]+Un[ix-2][iz])+a3*(Un[ix+3][iz]+Un[ix-3][iz])+a4*(Un[ix+4][iz]+Un[ix-4][iz])+a5*(Un[ix+5][iz]+Un[ix-5][iz])+a6*(Un[ix+6][iz]+Un[ix-6][iz]))

+0.5*v[ix][iz]*v[ix][iz]*dt*dt/(dz*dz)*(a0*Un[ix][iz]+a1*(Un[ix][iz+1]+Un[ix][iz-

1])+a2*(Un[ix][iz+2]+Un[ix][iz-2])+a3*(Un[ix][iz+3]+Un[ix][iz-3])+a4*(Un[ix][iz+4]+Un[ix][iz-4])+a5*(Un[ix][iz+5]+Un[ix][iz-5])+a6*(Un[ix][iz+6]+Un[ix][iz-6]));

}

Unp[nsx][nsz]+=source[it];

for(ix=0;ix