#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轴步长间隔
#define dz 10 //z轴步长间隔
#define xmax 500 //x轴总步长
#define zmax 300 //z轴总步长
#define T 1.5 // 传播时间长度
#define M 5
/************************速度模型******************************/
float v_m(int i,int j)
{
if(j<1000)
return 2000;
else if(j<2000)
return 2500;
else
return 3000;
/*if(j<-0.2*i+1500)
return 2000;
else if(j<0.2*i+1500)
return 2500;
else
return 3000;*/
}
void v_w()//写入速度模型文件
{
ofstream file;
int i,j;
float v;
file.open("velocity_model1.dat",ios::binary);
for(i=0;i<xmax;i++)
for(j=0;j<zmax;j++)
{
v=(float)v_m(i*dx,j*dz);
file.write((char*)&v,sizeof(float));
}
file.close();
}
/********************子波函数****************************/
float Ricker(double t) //ricker子波
{
float k;
double x=pow(PI*f0*(t-t0),2);
return k=A*exp(-x)*(1-2*x);
}
double gauss(double t)//gauss子波
{
double k;
double x=pow(PI*f0*(t-t0),2);
return k=-A*x*exp(-x)*(3-2*x)/(t-t0);
}
/*************************差分系数的计算*********************************/
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 C(float c[],int m)//系数函数
{
float aa[20][20],b[20]={0};
for(int i=0;i<m;i++)
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},c[M+1]={0};
float tx[xmax+2*M][1550];
float tz[zmax+2*M][1551];
float v[xmax+2*M][zmax+2*M];
C(c,M);
v_w();
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);}
}
/********************声波的模拟计算************************/
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[xmax/2+M][M]+=Ricker(k*dt);//加入震源子波
for(i=M;i<xmax+M;i++)
{
for(j=M;j<zmax+M;j++)
{
tz[j][k]=p3[250][j];//对vsp记录值的记录
tx[i][k]=p3[i][50];//suface 50
p1[i][j]=p2[i][j]; //数组的替换
p2[i][j]=p3[i][j];
}
}
/*****************************海绵边界************************************/
/*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<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)*/
/**************记录值的二进制文件写入******************/
if(k==400||k==700||k==1100)//对在k=400、700、1100时,波的写入
{
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++)//对vsp的写入
for(j=M;j<zmax+M;j++)
{
out1.write((char*) &p2[i][j],sizeof(float));
}
out1.close();
}
}
ofstream out2("tx.dat",ios::binary);//对surface的写入
for(i=M;i<xmax+M;i++)
for(k=0;k*dt<T;k++)
{
out2.write((char*) &tx[i][k],sizeof(float));
}
out2.close();
ofstream out3("tz.dat",ios::binary);
for(j=M;j<zmax+M;j++)
for(k=0;k*dt<T;k++)
{
out3.write((char*) &tz[j][k],sizeof(float));
}
out3.close();
return 0;
}
评论1