yxcf.cpp

有限差分二阶二维声波方程加吸收边界变速剖面研究.c++实现！

//有限差分二阶二维声波方程加吸收边界变速剖面研究
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>       //输入输出格式控制
#include <fstream.h>
#define F 30.0            //雷克子波主频，单位:hz
#define N0 2            //震源起爆时间
#define NZ 50            //纵向最大采样点号
#define NX 50            //横向最大采样点号
#define DT 0.002          //时间采样间隔，单位：s
#define NT 150          //时间最大样点号
#define DZ 12             //纵向采样间隔
#define DX 12             //横向采样间隔
#define XS 25             //震源位置横坐标样点号
#define ZS 6             //震源位置纵坐标样点号
#define C 3000            //波传播速度，单位：m/s
#define N 4           //地震记录深度样点
double u0[NX][NZ],u1[NX][NZ],u2[NX][NZ];
double out[NX][NT];
double wavelet(double f,int t0,int t,double dt)
{  double wav,a,b;
   b=dt*(t-t0);
   a=pow(3.14,2)*pow(f,2)*pow(b,2)*(-1);
   wav=(1+2*a)*exp(a);
   return wav;
 }
double shot(int e)
{  double m,n,mn,s[NT];

   m=pow(C,2)*pow(DT,2);
   n=DX*DZ;
   s[e]=wavelet(F,N0,e,DT);
   mn=m*s[e]/n;
   return mn;
}

void printout()
{ int i,k;
   ofstream outfile("p.sgy",ios::binary);
  
   outfile<<setiosflags(ios::left);
   outfile<<"p.sgy";    //从第1道到第51道;  采样点数：151;  采样间隔：2毫秒。
   outfile<<endl;
  for(i=0;i<=NX-1;i++)
   {  for(k=0;k<=NT-1;k++)
        outfile<<setw(20)<<out[i][k];
      outfile<<endl;
   }
    cout<<"write the file success!"<<"\n";
    outfile.close();
}
void main()
{ int i,k,n,data;
 
  double rx,rz;
  double a[NX][NZ],b[NX][NZ];
  int c[NX][NZ];
 
  //第一步：速度文件的载入及相关整理（替换） 
  
  rx=DT/DX;
  rz=DT/DZ;
  //读速度文件
   ifstream infile;
   infile.open("v.dat");
   for(i=0;i<NX;i++)
   for(k=0;k<NZ;k++)
   {  infile>>data;
      c[i][k]=data;
   }
   
   infile.close();
   for(i=0;i<=NX-1;i++)
  for(k=0;k<=NZ-1;k++)
  {    a[i][k]=pow(c[i][k],2)*pow(rx,2);
       b[i][k]=pow(c[i][k],2)*pow(rz,2);
  }
    

  for(n=0;n<=NT-1;n++)
  { //第二步：赋初值,初始时刻的全波场值均为零，P(i, j, 0)=0）
    //      时刻dt时，在炮点S (XS, ZS)给出一个脉冲震源S(t)，其它点波场P =0； 
    //初始条件t=0时
    if(n==0)
     {for(i=0;i<=NX-1;i++)
      for(k=0;k<=NZ-1;k++)
         u0[i][k]=0.0;

      }
    else if(n==1)
    { for(i=0;i<=NX-1;i++)
      for(k=0;k<=NZ-1;k++)
         u1[i][k]=0.0;

    }
    //第三步：边界条件处理及7点差分计算波场延拓  
    //当n>=2时
    else if(n>=2)
    {
       for(i=1;i<=NX-2;i++)
       {  //上边界吸收
          u2[i][0]=(2-2*rz*c[i][0]-b[i][0])*u1[i][0]+2*rz*c[i][0]*(1+rz*c[i][0])*u1[i][1]-b[i][0]*u1[i][2]+(2*rz*c[i][0]-1)*u0[i][0]-2*rz*c[i][0]*u0[i][1];
          //下边界吸收 
          u2[i][NZ-1]=(2-2*rz*c[i][NZ-1]-b[i][NZ-1])*u1[i][NZ-1]+2*rz*c[i][NZ-1]*(1+rz*c[i][NZ-1])*u1[i][NZ-2]-b[i][NZ-1]*u1[i][NZ-3]+(2*rz*c[i][NZ-1]-1)*u0[i][NZ-1]-2*rz*c[i][NZ-1]*u0[i][NZ-2];

        }
       for(k=1;k<=NZ-2;k++)
       {  //左边界吸收 
          u2[0][k]=(2-2*c[0][k]*rx-a[0][k])*u1[0][k]+2*c[0][k]*rx*(1+rx*c[0][k])*u1[1][k]-a[0][k]*u1[2][k]+(2*rx*c[0][k]-1)*u0[0][k]-2*rx*c[0][k]*u0[1][k];
          //右边界吸收 
          u2[NX-1][k]=(2-2*rx*c[NX-1][k]-a[NX-1][k])*u1[NX-1][k]+2*rx*c[NX-1][k]*(1+rx*c[NX-1][k])*u1[NX-2][k]-a[NX-1][k]*u1[NX-3][k]+(2*rx*c[NX-1][k]-1)*u0[NX-1][k]-2*rx*c[NX-1][k]*u0[NX-2][k];

        }
       for(i=1;i<=NX-2;i++)
       for(k=1;k<=NZ-2;k++)
       { if((k==ZS&&i==XS)==1)       //炮点S[ZS][XS]条件  
            u2[i][k]=shot(n-1)+2*(1-a[i][k]-b[i][k])*u1[i][k]+a[i][k]*(u1[i+1][k]+u1[i-1][k])+b[i][k]*(u1[i][k+1]+u1[i][k-1])-u0[i][k];
         else
            u2[i][k]=2*(1-a[i][k]-b[i][k])*u1[i][k]+a[i][k]*(u1[i+1][k]+u1[i-1][k])+b[i][k]*(u1[i][k+1]+u1[i][k-1])-u0[i][k];
        }
        //第四步：四个角点的处理 
      
       u2[0][0]=0.5*u2[1][0]+0.5*u2[0][1];
       u2[NX-1][0]=0.5*u2[NX-2][0]+0.5*u2[NX-1][1];
       u2[0][NZ-1]=0.5*u2[0][NZ-2]+0.5*u2[1][NZ-1];
       u2[NX-1][NZ-1]=0.5*u2[NX-2][NZ-1]+0.5*u2[NX-1][NZ-2];
       for(i=0;i<=NX-1;i++)
           out[i][n]=u2[i][N];
         // 数组循环覆盖 
       for(i=0;i<=NX-1;i++)
       for(k=0;k<=NZ-1;k++)
        u0[i][k]=u1[i][k];
       for(i=0;i<=NX-1;i++)
       for(k=0;k<=NZ-1;k++)
        u1[i][k]=u2[i][k];
    }
  }
   for(i=0;i<=NX-1;i++)
      {  out[i][0]=0.0;
         out[i][1]=0.0;
      }
   printout();
   
}