#include
#include
#include
#include
#define pi 3.1415926
#define fmain 15
void make_ricker(float dt,float lt,int wave_length,float *ricker);
void make_am(float *am);
void extrap(float **sv_record,float *am,float **wf1,float **wf2,float **wf3,int nx,int nz,float dx,float dt,float dz);
void Read_sv_record_boundary(char *filename,float **sv_record,int nx,int nz);
void output(float **wf_record,int nx,int nz);
void output1(float **se_record,int nx,int lt);
main()
{
int nx,nz;
int ix,iz,it;
int sz,sx;
float dx,dz,dt;
float lt,time_record;
nx=401,nz=301,dx=10,dz=10,dt=1;
lt=5000,time_record=1000;
sx=200,sz=0;//paodian//
float **v1=NULL,**v2=NULL;
v1=alloc2float(nx,nz);
v2=alloc2float(nx,nz);
float sv_record=NULL;
Read_sv_record_boundary("./diqian-401-301",sv_record,nx,nz);
float **wf1,wf2,**wf3,**wf_record,**se_record;
wf1=alloc2float(nx,nz);
wf2=alloc2float(nx,nz);
wf3=alloc2float(nx,nz);
wf_record=alloc2float(nx,nz);
se_record=alloc2float(nx,lt+100);
for(iz=0;iz<nz;iz++)
for(ix=0;ix<nx;ix++)
{
wf1[iz][ix]=0;
wf2[iz][ix]=0;
wf3[iz][ix]=0;
wf_record[iz][ix]=0;
}
for(it=0;it<lt+100;it++)
for(ix=0;ix<nx;ix++)
{
se_record[it][ix]=0;
}
int wave_length=200;
float *ricker;
ricker=alloc1float(lt);
make_ricker(dt,lt,wave_length,ricker);
float *am;
am=alloc1float(6);
make_am(am);
for(it=0;it<lt+100;it++)
{
if(it<=200)
{
wf2[sz][sx]=wf2[sz][sx]+ricker[it];
}
extrap(v2,am,wf1,wf2,wf3,nx,nz,dx,dt,dz);
if(it==time_record);
{
for(iz=0;iz<nz;iz++)
for(ix=0;ix<nx;ix++)
wf_record[iz][ix]=wf2[iz][ix];
}
for(ix=0;ix<nx;ix++)
se_record[it][ix]=wf3[100][ix];
for(iz=0;iz<nz;iz++)
for(ix=0;ix<nx;ix++)
{
wf1[iz][ix]=wf2[iz][ix];
wf2[iz][ix]=wf3[iz][ix];
}
}
output(wf_record,nx,nz);
output1(se_record,nx,lt);
free1float(ricker);
free2float(sv_record);
free1float(am);
free2float(wf1);
free2float(wf2);
free2float(wf3);
free2float(wf_record);
free2float(se_record);
}
void Read_sv_record_boundary(char *filename,float **sv_record,int nx,int nz)
{
int ix,iz;
FILE*fp;
if((fp=fopen(filename,"r"))==NULL)
{
printf("cannot open this file\n,");
exit(0);
}
for(ix=0;ix<nx;ix++)
{
for(iz=0;iz<nz;iz++)
{
fread(&sv_record[iz][ix],4,1,fp);
}
}
fclose(fp);
}
void make_ricker(float dt,float lt,int wave_length,float *ricker)
{
float sdt,time,tp1,tp2;
int it;
sdt=dt*0.001;
for(it=0;it<lt;it++)
{
time=(it+1)*sdt-(wave_length/2.0)*sdt;
tp1=pi*fmain*time;
tp2=tp1*tp1;
ricker[it]=(1.0-2.0*tp2)*exp(-tp2);
}
}
void make_am(float *am)
{
am[0]=-2.92722;
am[1]=1.66667;
am[2]=-0.238095;
am[3]=0.0396825;
am[4]=-0.00496032;
am[5]=0.000317460;
}
void extrap(float **sv_record,float am,float **wf1,float **wf2,float **wf3,int nx,int nz,float dx,float dt,float dz)
{
int ix,iz;
float dt_r,a,b;
dt_r=dt*0.001;
for(iz=5;iz<nz-5;iz++)
for(ix=5;ix<nx-5;ix++)
{
wf3[iz][ix]=2*wf2[iz][ix]-wf1[iz][ix]+
sv_record[iz][ix]*sv_record[iz][ix]*dt_r*dt_r/(dz*dz)
(am[0]*wf2[iz][ix]+
am[1]*(wf2[iz+1][ix]+wf2[iz-1][ix])+
am[2]*(wf2[iz+2][ix]+wf2[iz-2][ix])+
am[3]*(wf2[iz+3][ix]+wf2[iz-3][ix])+
am[4]*(wf2[iz+4][ix]+wf2[iz-4][ix])+
am[5]*(wf2[iz+5][ix]+wf2[iz-5][ix])
)+
sv_record[iz][ix]*sv_record[iz][ix]*dt_r*dt_r/(dx*dx)*
(am[0]*wf2[iz][ix]+
am[1]*(wf2[iz][ix+1]+wf2[iz][ix-1])+
am[2]*(wf2[iz][ix+2]+wf2[iz][ix-2])+
am[3]*(wf2[iz][ix+3]+wf2[iz][ix-3])+
am[4]*(wf2[iz][ix+4]+wf2[iz][ix-4])+
am[5]*(wf2[iz][ix+5]+wf2[iz][ix-5]));
}
}
void output(float **wf_record,int nx,int nz)
{
int ix,iz;
FILE *fp2;
fp2=fopen("wf_record.bin","wb+");
for(ix=0;ix<nx;ix++)
for(iz=0;iz<nz;iz++)
{
fwrite(&wf_record[iz][ix],4,1,fp2);
}
fclose(fp2);
}
void output1(float **se_record,int nx,int lt)
{
int ix,it;
FILE *fp1;
fp1=fopen("se_record.bin","wb+");
for(ix=0;ix<nx;ix++)
for(it=0;it<lt+100;it++)
{
fwrite(&se_record[it][ix],4,1,fp1);
}
fclose(fp1);
}
https://blog.csdn.net/binbigdata/article/details/80630524
正演代码,生成了波场快照与地震数据,但你这个代码我无法运行