qq_38860311
qq_38860311
2020-09-14 15:47

正演程序 希望帮忙注释下 跨专业研1新生

  • c语言

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

  • 点赞
  • 回答
  • 收藏
  • 复制链接分享

1条回答