acoustic forward1.zip_交错_地震正演_声波交错_声波正演模拟_空间 阶

#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); }