#include "su.h"
#include "segy.h"
#include <signal.h>
#include <time.h>
/*********************** self documentation ******************************/
char *sdoc[] = {
" ",
" this is a program to model multiple ",
" this was made by lsz in Daqing in 2012-03-18 ",
" ",
NULL};
/**************** end self doc *******************************************/
main(int argc, char **argv)
{
float dt=0.004,tt=0.4,dx=20.0;
//float dz1=400,dz2=800,v1=1500,p1=2500,v2=2000,p2=2800,v3=3000,p3=3200;
float dz[3]={400,800,600},v[4]={1500,2000,3000,3200},p[4]={2500,2800,3200,3300};
float kz,a,dw,dk,r1,r2;
int i,j,iw,ix,it,nt0,nx=500,nt=1024,ntfft,nxfft,nw,resultx=200,resultt=1024,nj=4,l=3;
float *ricker,**wavef,**result,*r;
float **wavefbegin,**wavefend;
complex **datawx,**dataxw,**datawx2;
complex x,**rr,**rr1,**rr2;
FILE *fp1;
initargs(argc, argv);
/********cite a function*******/
nt0=(int)(tt/dt);
ricker=alloc1float(nt0);
for(i=0;i<nt0;i++) ricker[i]=0.0;
rickersrc(ricker,nt0);
/********set up wave phenomena whose value is zero****************/
wavef=alloc2float(nt,nx);
for(ix=0;ix<nx;ix++)
for(it=0;it<nt;it++)
wavef[ix][it]=0.0;
/**********cite the ricker*************/
for(ix=0;ix<nx;ix++)
if(ix==nx/2)
{for(it=0;it<100;it++)
wavef[ix][it]=ricker[it];}
/********* canshu zhunbei ******************/
warn("000000000000000000000000");
ntfft=npfar(nt);
nw=ntfft/2+1;
nxfft=npfar(nx);
warn("ntfft=%d, nw=%d, nxfft=%d\n",ntfft,nw,nxfft);
datawx=alloc2complex(nw,nxfft);
dataxw=alloc2complex(nxfft,nw);
datawx2=alloc2complex(nw,nxfft);
wavefend=alloc2float(ntfft,nxfft);
wavefbegin=alloc2float(ntfft,nxfft);
for(ix=0;ix<nxfft;ix++)
for(iw=0;iw<nw;iw++)
datawx[ix][iw]=cmplx(0.0,0.0);
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
dataxw[iw][ix]=cmplx(0.0,0.0);
for(ix=0;ix<nxfft;ix++)
for(iw=0;iw<nw;iw++)
datawx2[ix][iw]=cmplx(0.0,0.0);
for(ix=0;ix<nxfft;ix++)
for(it=0;it<ntfft;it++)
wavefend[ix][it]=0.0;
for(ix=0;ix<nxfft;ix++)
for(it=0;it<ntfft;it++)
wavefbegin[ix][it]=0.0;
for(ix=0;ix<nx;ix++)
for(it=0;it<nt;it++)
wavefbegin[ix][it]=wavef[ix][it];
/************* 2Dfft ***********************/
pfa2rc(-1,1,ntfft,nxfft,wavefbegin[0],datawx[0]);//fft1
for(ix=0;ix<nxfft;ix++)
for(iw=0;iw<nw;iw++)
{x=datawx[ix][iw];
dataxw[iw][ix]=x;}
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
dataxw[iw][ix]=ix%2? cneg(dataxw[iw][ix]):dataxw[iw][ix];
pfa2cc(-1,1,nxfft,nw,dataxw[0]);//fft2
/************* model mutiple ***************/
rr=alloc2complex(nxfft,nw);
rr1=alloc2complex(nxfft,nw);
rr2=alloc2complex(nxfft,nw);
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
{
rr[iw][ix]=cmplx(0.0,0.0);
rr1[iw][ix]=cmplx(0.0,0.0);
rr2[iw][ix]=cmplx(0.0,0.0);
}
r=alloc1float(l);
for(i=0;i<l;i++) r[i]=fabs(p[i+1]*v[i+1]-p[i]*v[i])/(p[i+1]*v[i+1]+p[i]*v[i]);
dw = 2.0*PI/(ntfft*dt);
dk = PI/(nxfft*dx);
warn("dw=%f dk=%f\n",dw,dk);
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
rr1[iw][ix]=dataxw[iw][ix];
for(i=0;i<=nj;i++) //nj multiple
{
for(j=0;j<l;j++) //layer
{for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
{
a=(iw*dw*iw*dw)/(v[j]*v[j])-((ix-nxfft/2)*dk*(ix-nxfft/2)*dk);
if(a>=0)
{
kz=sqrt(a);
x=cexp(cmplx(0.0,-kz*dz[j]));
rr1[iw][ix]=cmul(x,rr1[iw][ix]);
}
else
{
a=-a;
kz=sqrt(a);
x=cmplx(exp(-kz*dz[j]),0.0);
rr1[iw][ix]=cmul(x,rr1[iw][ix]);
}
}
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
rr2[iw][ix]=cmul(cmplx(r[j],0.0),rr1[iw][ix]);
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
rr[iw][ix]=cadd(rr[iw][ix],rr2[iw][ix]);
}
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
rr1[iw][ix]=cmul(cmplx(-1.0,0.0),rr2[iw][ix]);
}
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
dataxw[iw][ix]=rr[iw][ix];
free2complex(rr);
free2complex(rr1);
free2complex(rr2);
/************* 2Difft *************/
pfa2cc(1,1,nxfft,nw,dataxw[0]);//-fft1
for(iw=0;iw<nw;iw++)
for(ix=0;ix<nxfft;ix++)
{x=cdiv(dataxw[iw][ix],cmplx((float)nxfft,0.0));
datawx2[ix][iw]=x;}
for(ix=0;ix<nxfft;ix++)
for(iw=0;iw<nw;iw++)
datawx2[ix][iw]=ix%2? cneg(datawx2[ix][iw]):datawx2[ix][iw];
pfa2cr(1,1,ntfft,nxfft,datawx2[0],wavefend[0]);//-fft2
for(ix=0;ix<nxfft;ix++)
for(it=0;it<ntfft;it++)
wavefend[ix][it]=wavefend[ix][it]/(float)ntfft;
result=alloc2float(resultt,resultx);
for(ix=0;ix<resultx;ix++)
for(it=0;it<resultt;it++)
result[ix][it]=0.0;
for(ix=0;ix<resultx;ix++)
for(it=0;it<resultt;it++)
result[ix][it]=wavefend[ix+150][it];
fp1=fopen("result.bin","wb");
for(ix=0;ix<resultx;ix++)
for(it=0;it<resultt;it++)
fwrite(&result[ix][it],sizeof(float),1,fp1);
fclose(fp1);
free2float(wavefbegin);
free1float(ricker);
free2float(wavef);
free2float(wavefend);
free2float(result);
free2complex(datawx);
free2complex(dataxw);
free2complex(datawx2);
warn("end");
}
void rickersrc(float *ricker,int nt0)
{
float a,b,t=0.0,tt=0.4,dt=0.004;
int i;
float fp=25;
for(i=0;i<nt0;i++)
ricker[i]=0.0;
for(i=0;i<nt0;i++)
{
t=(i-nt0/2)*dt;
a=PI*t*fp;
b=a*a;
ricker[i]=(1-2*b)*exp(-b);
}
}
zhengyan.zip_地震数据_地震数据正演_地震正演_雷克子波_雷克子波程序
版权申诉
5星 · 超过95%的资源 79 浏览量
2022-09-23
12:48:45
上传
评论
收藏 2KB ZIP 举报
alvarocfc
- 粉丝: 109
- 资源: 1万+
最新资源
- 129335283047061xiazaigongjuxiang(去重软件).apk
- Android环境检测工具,检测ksu,lsp,magisk等
- WordPress后台美化插件QuarterAdmin分享
- PCB_Project单片机绘制 (2024-5-11 22-22-13).zip
- nccl-local-repo-ubuntu2204-2.21.5-cuda12.4-1.0-1-amd64
- STM32移植LVGL源码工程 stm32移植GUI-Guider源码 LVGL移植源码
- Screenshot_20240522_084328_com.tencent.mm.jpg
- 附件计算机专业课选课说明-1.xlsx
- 基于TypeScript的ahousepet-admin-web管理系统设计源码
- 《广东开放大学学习指引》期末考核要求0522.zip
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈