sss.rar_弹性方程_弹性波

#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; }