void Alloc()
{
wavelet=alloc1float(nt); memset(wavelet,0,nt*sizeof(float));
velocity=alloc1float(nx*nz); memset(velocity,0,nx*nz*sizeof(float));
wf1=alloc1float(nx*nz); memset(wf1,0,nx*nz*sizeof(float));
wf2=alloc1float(nx*nz); memset(wf2,0,nx*nz*sizeof(float));
wf3=alloc1float(nx*nz); memset(wf3,0,nx*nz*sizeof(float));
coe=alloc1float(N/2+1); memset(coe,0,(N/2+1)*sizeof(float));
}
void output(char *filename,float *data,int nx,int nz)
{
FILE *fp;
if((fp=fopen(filename,"wr+"))==NULL)
{ printf("can not find '*filename'\n");
exit(0);
}
else
{
fwrite(data,nx*nz*sizeof(float),1,fp);
}
fclose(fp);
}
void output_wavelet(char *filename,float *data,int nt)
{
FILE *fp;
if((fp=fopen(filename,"wr+"))==NULL)
{ printf("can not find '*filename'\n");
exit(0);
}
else
{
fwrite(data,nt*sizeof(float),1,fp);
}
fclose(fp);
}
void generate_velocity(float *data,int nx,int nz)
{
int ix,iz;
for(ix=0;ix<nx;ix++)
{
for(iz=0;iz<nz;iz++)
{
if(iz<150)
data[ix*nz+iz]=2500.0;
else
data[ix*nz+iz]=3000.0;
}
}
}
void generate_wavelet(float *data,int nt,float dt,float fm)
{
float temp;
int it,itshift;
for(it=0;it<nt;it++)
{
itshift=it-(int)(1.0/fm/dt);
temp=PI*PI*fm*fm*(itshift*1.0)*dt*(itshift*1.0)*dt;
data[it]=(1.0-2.0*temp)*exp(-1.0*temp);
}
}
void add_wavelet(float *wf,float *walelet,int nx,int nz,int it,int sx,int sz)
{
wf[sx*nz+sz]=wf[sx*nz+sz]+wavelet[it];
}
void generate_coe(float *coe)
{
coe[0]=-2.0;
coe[1]=1.0;
}
void exteapolation(float *wf1,float *wf2,float *wf3,float *velocity,int nx,int nz,float dx,float dz,float dt,float *coe,int N)
{
int ix,iz;
for(ix=N/2;ix<nx-N/2;ix++)
for(iz=N/2;iz<nz-N/2;iz++)
{
/*float temp1=velocity[ix*nz+iz]*velocity[ix*nz+iz]*dt*dt/dz/dz;
float temp2=velocity[ix*nz+iz]*velocity[ix*nz+iz]*dt*dt/dx/dx;
wf3[ix*nz+iz]=2.0*wf2[ix*nz+iz]-wf1[iz*nz+iz]+temp1*(coe[1]*wf2[ix*nz+iz-1]+coe[0]*wf2[ix*nz+iz]+coe[1]*wf2[ix*nz+iz+1])+temp2*(coe[1]*wf2[(ix-1)*nz+iz]+coe[0]*wf2[ix*nz+iz]+coe[1]*wf2[(ix+1)*nz+iz]);*/
wf3[ix*nz+iz]=2.0*wf2[ix*nz+iz]-wf1[ix*nz+iz]+velocity[ix*nz+iz]*velocity[ix*nz+iz]*dt*dt/dz/dz*(coe[1]*(wf2[ix*nz+iz-1]+wf2[ix*nz+iz+1])+coe[0]*wf2[ix*nz+iz])+velocity[ix*nz+iz]*velocity[ix*nz+iz]*dt*dt/dx/dx*(coe[1]*(wf2[(ix-1)*nz+iz]+wf2[(ix+1)*nz+iz])+coe[0]*wf2[ix*nz+iz]);
}
}
void replace_wf(float *wf1,float *wf2,float *wf3,int nx,int nz)
{
int ix,iz;
for(ix=0;ix<nx;ix++)
{
for(iz=0;iz<nz;iz++)
{
wf1[ix*nz+iz]=wf2[ix*nz+iz];
wf2[ix*nz+iz]=wf3[ix*nz+iz];
}
}
}
void zhengyan_display(float *wf1,float *wf2,float *wf3,float *velocity,int nx,int nz,int nt,float dx,float dz,float dt,float *coe,int N)
{
int it;
for(it=0;it<nt;it++)
{
add_wavelet(wf2,wavelet,nx,nz,it,sx,sz);
exteapolation(wf1,wf2,wf3,velocity,nx,nz,dx,dz,dt,coe,N);
replace_wf(wf1,wf2,wf3,nx,nz);
if(it==200)
output("zhengyan.bin",wf2,nx,nz);
}
}
void Free()
{
free1float(wavelet);
free1float(velocity);
free1float(wf1);
free1float(wf2);
free1float(wf3);
free1float(coe);
}
