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

#include#include#define fm 30#define dt 0.001#define PI 3.#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;iNt;i+) sourcei=0.0;t0 = 1.0/fm;printf(t0 = %fn, t0); fp=fopen(ricker.txt,w);for(it=0;itNt;it+)t1=it*dt;t2=t1-t0;sourceit=(1.0-2.0*pow(PI*fm*t1,2.0)*exp(-pow(PI*fm*t1,2.0);fprintf(fp,%.8f %.8fn,t1,sourceit);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 UnNxNz, UnmNxNz, UnpNxNz, vNxNz,recordNxNt;/Un表示n时刻，Unm表示n-1时刻，Unp表示n+1时刻波场 float sourceNt;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;ixNx;ix+) for (iz=0;iz(Nz*2/5);iz+) vixiz=2000;for (iz=(Nz*2/5);iz(Nz*3/5);iz+) vixiz=3000;for (iz=(Nz*3/5);izNz;iz+) vixiz=4000;/给速度赋值 */ for (ix=0;ixNx;ix+) for (iz=0;izNz/2;iz+)vixiz=4000; for (iz=Nz/2;izNz;iz+)vixiz=4000; dx=10.0;dz=10.0;/赋值for(ix=0;ixNx;ix+)for(iz=0;izNz;iz+)Unixiz=0.0;Unmixiz=0.0;Unpixiz=0.0;fun(source);fpwf = fopen(Wavefieldall.dat, wb);for(it=0;itNt;it+)for(ix=6;ix(Nx-6);ix+)for(iz=6;iz(Nz-6);iz+) Unpixiz=2*Unixiz-Unmixiz+0.5*vixiz*vixiz*dt*dt/(dx*dx)*(a0*Unixiz+a1*(Unix+1iz+Unix-1iz)+a2*(Unix+2iz+Unix-2iz)+a3*(Unix+3iz+Unix-3iz)+a4*(Unix+4iz+Unix-4iz)+a5*(Unix+5iz+Unix-5iz)+a6*(Unix+6iz+Unix-6iz) +0.5*vixiz*vixiz*dt*dt/(dz*dz)*(a0*Unixiz+a1*(Unixiz+1+Unixiz-1)+a2*(Unixiz+2+Unixiz-2)+a3*(Unixiz+3+Unixiz-3)+a4*(Unixiz+4+Unixiz-4)+a5*(Unixiz+5+Unixiz-5)+a6*(Unixiz+6+Unixiz-6); Unpnsxnsz+=sourceit;for(ix=0;ixNx;ix+)for(iz=0;izNz;iz+) Unmixiz=Unixiz;Unixiz=Unpixiz;for(ix=0;ixNx;ix+) recordixit=Unixnsz; /-记录文件-/-记录Un文件-for(ix=0; ixNx; ix+)fwrite(Unix, sizeof(float), Nz, fpwf);if(it=400) fp=fopen(wavefield.dat,wb);for(ix=0; ixNx; ix+)fwrite(Unix, sizeof(float), Nz, fp);fclose(fp);if(it%20=0) printf(it=%dn, it);fclose(fpwf); /-记录record文件- fp=fopen(record.dat,wb); for(ix=0;ixNx;ix+) fwrit