//有限差分二階二維聲波方程加吸收邊界變速剖面研究
#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();
   
}