#include <iostream>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <fstream>
using namespace std;
#define PI 3.1415926
#define f0 20 //子波主频
#define t0 0.08 //子波延迟,单位：毫秒
#define A 100 //振幅
#define dt 0.001 //时间间隔
#define dx 10 //x间隔 米 共2000m
#define dz 10 //z间隔 米 共2000m
#define xmax 400
#define zmax 400
#define T 1.5 //改数组大小
#define M 5
float v_m(int i,int j)//速度模型
{ return 2400;
}
void l_gauss(float a[20][20],float b[],float x[],int n)//列主元高斯消去法
{
int i,j,k;
float B[20][20];
for(i=0;i<n;i++)//生成增广矩阵
{
B[i][n]=b[i];
for(j=0;j<n;j++)
B[i][j]=a[i][j];
}
float q,t;
for(k=0;k<n;k++)
{
int imax,p;
imax=k;
for(p=k+1;p<n;p++)//找出列主元下标
if(B[p][k]>B[imax][k])
imax=p;
for(p=k;p<=n;p++)//将列主元所在的行换到首行
{
t=B[imax][p];
B[imax][p]=B[k][p];
B[k][p]=t;
}
{for(i=k+1;i<n;i++)//化成上三角阵
{
q=B[i][k]/B[k][k];
for(j=k;j<=n;j++)
B[i][j]=B[i][j]-q*B[k][j];
}
}
float s(0);
for(i=n-1;i>=0;i--)//回代得到x
{
s=0;
for(j=i+1;j<n;j++)
s=s+B[i][j]*x[j];
x[i]=(B[i][n]-s)/B[i][i];
}
}
}
float Ricker(double t) //ricker子波公式
{
float k;
double x=pow(PI*f0*(t-t0),2);
return k=A*exp(-x)*(1-2*x);
}
float C(float c[],int m)//计算2M阶差分系数
{
float aa[20][20],b[20]={0};
for(int i=0;i<m;i++)//生成矩阵AA
for(int j=0;j<M;j++)
{
aa[i][j]=pow(j+1,2*(i+1));
}
b[0]=1;
l_gauss(aa,b,c,m);
return 1;
}
int main()//计算个时间面上的值
{
int i,j,k;
float p1[xmax+2*M][zmax+2*M]={0};
float p2[xmax+2*M][zmax+2*M]={0};
float p3[xmax+2*M][zmax+2*M]={0};
float c[M+1]={0};
float tx[xmax+2*M][1550],tz[zmax+2*M][1551];
float v[xmax+2*M][zmax+2*M];
C(c,M);
for(i=0;i<xmax+2*M;i++)
{
for(j=0;j<zmax+2*M;j++)
{
v[i][j]=v_m(i*dx,j*dz);}
}
/**************************ce****************************************/
for(k=0;k*dt<T;k++)
{
for(i=M;i<350+M;i++)
{
for(j=M;j<350+M;j++)
{
float asum=0,bsum=0;
for(int m=0;m<M;m++)
{
asum+=c[m] * (p2[i+m+1][j]-2*p2[i][j]+p2[i-m-1][j]);
bsum+=c[m] * (p2[i][j+m+1]-2*p2[i][j]+p2[i][j-m-1]);
}
p3[i][j]=v[i][j]*v[i][j]*dt*dt*(asum/(dx*dx)+bsum/(dz*dz))+2*p2[i][j]-p1[i][j];
}
}
for(i=350+M;i<xmax+M;i++)
for(j=M;j<zmax+M;j++)
p3[i][j]=(float)(-(p2[i][j]-p2[i-1][j])/dz*dt*v[i][j]+p2[i][j]);
for(j=350+M;j<zmax+M;j++)
for(i=M;i<xmax+M;i++)
p3[i][j]=(float)(-(p2[i][j]-p2[i][j-1])/dx*dt*v[i][j]+p2[i][j]);
p3[200+M][200+M]+=Ricker(k*dt);
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
p1[i][j]=p2[i][j];
p2[i][j]=p3[i][j];
}
}
/**************************************************************/
/*****************************sp************************************/
/*for(k=0;k*dt<T;k++)
{
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
float asum=0,bsum=0;
for(int m=0;m<M;m++)
{
asum+=c[m] * (p2[i+m+1][j]-2*p2[i][j]+p2[i-m-1][j]);
bsum+=c[m] * (p2[i][j+m+1]-2*p2[i][j]+p2[i][j-m-1]);
}
p3[i][j]=v[i][j]*v[i][j]*dt*dt*(asum/(dx*dx)+bsum/(dz*dz))+2*p2[i][j]-p1[i][j];
}
}
p3[180+M][180+M]+=Ricker(k*dt);
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
p1[i][j]=p2[i][j]; //替换数组
p2[i][j]=p3[i][j];
if(i>=200)
{float f=0.0000001;p2[i][j]=p2[i][j]*exp(-f*(i-200)*(i-200));}
if(j>=200)
{float f=0.0000001;p2[i][j]=p2[i][j]*exp(-f*(j-200)*(j-200));}
}
}
/*********************************************************************************/
/******************************************************************************/
/*for(k=0;k*dt<T;k++)
{
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
float asum=0,bsum=0;
for(int m=0;m<M;m++)
{
asum+=c[m] * (p2[i+m+1][j]-2*p2[i][j]+p2[i-m-1][j]);
bsum+=c[m] * (p2[i][j+m+1]-2*p2[i][j]+p2[i][j-m-1]);
}
p3[i][j]=v[i][j]*v[i][j]*dt*dt*(asum/(dx*dx)+bsum/(dz*dz))+2*p2[i][j]-p1[i][j];
}
}
p3[zmax/2+M][zmax/2+M]+=Ricker(k*dt);
dt2=dt*dt;v2=v[i-m][j-m]*v[i-m][j-m];
f=3*x*x*v/(2*bd*bd*bd*d*d*d)*log(1e3);
df=3*x*v/(bd*bd*bd*d*d*d)*log(1e3);
fx=calc_f(v[i-m][j-m],((i-m)-(nx-bd)+1)*dx,dx);
dfx=calc_df(v[i-m][j-m],((i-m)-(nx-bd)+1)*dx,dx);
pa3[i][j]=2*pa2[i][j]-pa1[i][j]+dt2*(v2*x_deri-2*fx*(pa2[i][j]-pa1[i][j])/dt-fx*fx*pa2[i][j]);
a3[i][j]=2*a2[i][j]-a1[i][j]+dt2*((-v2*dfx*(p2[i+1][j]-p2[i-1][j])/2.0/dx)-2*fx*(a2[i][j]-a1[i][j])/dt-fx*fx*a2[i][j]);
pb3[i][j]=(a2[i][j]-fx*pb2[i][j])*dt+pb2[i][j];
pc3[i][j]=dt2*v2*z_deri+2*pc2[i][j]-pc1[i][j];
p3[i][j]=pa3[i][j]+pb3[i][j]+pc3[i][j];
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
p1[i][j]=p2[i][j]; //替换数组
p2[i][j]=p3[i][j];
}
}
/******************************************************************************/
if(k==600||k==800||k==1100)//分别输出500/800/1200的数据
{
ostringstream v;
v<<k;
string s1,s2,s3;
s1=v.str();
s2=".dat";
s3=s1+s2;
ofstream out1(s3.c_str(),ios::binary);
for(i=M;i<xmax+M;i++)
for(j=M;j<zmax+M;j++)
{
out1.write((char*) &p2[i][j],sizeof(float));
}
out1.close();
}
}
return 0;
}
