C++实现合成地震记录

#include<iostream>
#include<math.h>
#include<fstream>
using namespace std;
#define pi 3.14
#define dt 0.002
#define xl 0.060
#define hl 0.300
#define fm 30
void main()
{
double h[200];
double t,m,n;
int i,j;
m=xl/dt+1;
n=hl/dt+1;
int b=(int)m/2;
double x[31], y[246];
cout<<"设计的层数为层"<<'\n'<<"各层密度分别为2.0 2.3 2.3 2.6 2.0"<<'\n';
cout<<"各层速度分别为2000 2500 2100 2700 3000"<<'\n'<<"各层厚度分别为100 100 100 100 100"<<'\n';
ofstream out1("wavelet.txt");
for(i=0;i<=15;i++)///////////////////生成雷克子波
{
t=i*dt;
x[15-i]=(1.0-2.0*pow(pi*fm*t,2.0))*exp(-pow(pi*fm*t,2.0));
x[15+i]=x[15-i];
// cout<<t*1000<<'\t'<<x[i+b]<<'\n';
}
for(i=0;i<31;i++)
{
cout<<x[i]<<'\n';
out1<<x[i]<<'\n';
}
out1<< flush; out1.close();
double p[5]={2.0, 2.3, 2.3, 2.6, 2.0};/////////////////////////////////生成反射系数double v[5]={2000, 2500, 2100, 2700, 3000};
double th[5]={100, 100, 100, 100, 100};
double *w1, *w2, *w3, *w4, r[4],a[4]={0};
w1=p;
w2=v;
w3=r;
for(i=0;i<4;i++)
{
r[i]=(w1[i+1]*w2[i+1]-w1[i]*w2[i])/(w1[i+1]*w2[i+1]+w1[i]*w2[i]);
// cout<<*(e+i);
// cout<<r[i]<<'\t';
}
a[0]=2*th[0]/(v[0]*dt);
for(i=1;i<4;i++)
{
a[i]=a[i-1]+2*th[i]/(v[i]*dt);
// cout<<a[i]<<" ";
}
w4=a;
ofstream out2("reflect.txt");
for(i=0;i<200;i++)
{
if(i==(int)w4[0])
h[i]=w3[0];
else if(i==(int)w4[1])
h[i]=w3[1];
else if(i==(int)w4[2])
h[i]=w3[2];
else if(i==(int)w4[3])
h[i]=w3[3];
else
h[i]=0.0;
t=i*dt;
out2<<t*1000<<'\t'<<h[i]<<'\n';
// cout<<h[i]<<'\n';
}
out2<<flush; out2.close();
// cout<<endl;
////////////////////////////////////////////实现卷积运算ofstream out3("convolution.txt");
for(i=0;i<246;i++)
{
y[i]=0.0;
}
for(i=0;i<231;i++)
{
y[i]=0.0;
for(j=0;j<m;j++)
{
if((i-j)>=0&&(i-j)<=200)
y[i]=y[i]+x[j]*h[i-j];
}
// cout<<t*1000<<'\t'<<y[i]<<'\n';
}
for(i=0;i<231;i++)
{
y[i]=y[i+15];
t=i*dt;
out3<<t*1000<<'\t'<<y[i]<<'\n';
}
out3<<flush; out3.close();
}。