//有限差分時(shí)間域二階、空間域四階二維聲波方程加吸收邊界常速模擬研究 (F=30  C=2500  時(shí)好)   震源不可放在邊界
  //
  //      u(i,k,n+1)=(p*p/12)*{16[u(i+1,k,n)+u(i-1,k,n)+u(i,k+1,n)+u(i,k-1,n)]-[u(i-2,k,n)+u(i+2,k,n)+u(i,k+2,n)+u(i,k-2,n)]}
  //                 +(2-5*p*p)*u(i,k,n)-u(i,k,n-1)+C*C*DT*DT*S(m,n,k)
  //
  //       p=C*DT/h
  //
  // 左右邊界條件
  //             for(k=0;k<=NZ-1;k++)
  //                {  u2[1][k]=u1[1][k]+u1[2][k]-u0[2][k]+C*DT/DX*((u1[2][k]-u1[1][k])-(u0[3][k]-u0[2][k]));
  //                   u2[0][k]=u1[0][k]+u1[1][k]-u0[1][k]+C*DT/DX*((u1[1][k]-u1[0][k])-(u0[2][k]-u0[1][k]));
  //                   u2[NZ-2][k]=u1[NZ-2][k]+u1[NZ-3][k]-u0[NZ-3][k]+C*DT/DX*((u1[NZ-3][k]-u1[NZ-2][k])-(u0[NZ-4][k]-u0[NZ-3][k]));
  //                   u2[NZ-1][k]=u1[NZ-1][k]+u1[NZ-2][k]-u0[NZ-2][k]+C*DT/DX*((u1[NZ-2][k]-u1[NZ-1][k])-(u0[NZ-3][k]-u0[NZ-2][k]));
  //                }
  //              上下邊界條件
  //             for(i=0;i<=NX-1;i++)
  //                {  u2[i][1]=u1[i][1]+u1[i][2]-u0[i][2]+C*DT/DZ*((u1[i][2]-u1[i][1])-(u0[i][3]-u0[i][2]));
  //                   u2[i][0]=u1[i][0]+u1[i][1]-u0[i][1]+C*DT/DZ*((u1[i][1]-u1[i][0])-(u0[i][2]-u0[i][1]));
 //                   u2[i][NZ-2]=u1[i][NZ-2]+u1[i][NZ-3]-u0[i][NZ-3]+C*DT/DZ*((u1[i][NZ-3]-u1[i][NZ-2])-(u0[i][NZ-4]-u0[i][NZ-3]));
  //                   u2[i][NZ-1]=u1[i][NZ-1]+u1[i][NZ-2]-u0[i][NZ-2]+C*DT/DZ*((u1[i][NZ-2]-u1[i][NZ-1])-(u0[i][NZ-3]-u0[i][NZ-2]));
  //                }
  //                         

         
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>       //輸入輸出格式控制
#include <fstream.h>
#define F 30.0            //雷克子波主頻，單位:hz 
#define N0 2            //震源起爆時(shí)間 
#define NZ 51             //縱向最大采樣點(diǎn)號(hào) 
#define NX 51             //橫向最大采樣點(diǎn)號(hào) 
#define DT 0.002          //時(shí)間采樣間隔，單位：s 
#define NT 100             //時(shí)間最大樣點(diǎn)號(hào) 
#define DZ 9             //縱向采樣間隔 
#define DX 9             //橫向采樣間隔 
#define XS 25             //震源位置橫坐標(biāo)樣點(diǎn)號(hào) 
#define ZS 25             //震源位置縱坐標(biāo)樣點(diǎn)號(hào) 
#define C 2500            //波傳播速度，單位：m/s 
#define T1 80          //輸出t=T1*DT時(shí)的切片 


double u0[NX][NZ],u1[NX][NZ],u2[NX][NZ];
double wavelet(double f,int t0,int t,double dt)
{double wav,a,b;
 b=dt*(t-t0);
 a=pow(3.1415,2)*pow(f,2)*pow(b,2)*(-1);
 wav=(1+2*a)*exp(a);
 return wav;
 }
double shot(int a,int b,int e)
{double m,n,mn,s[NT];
  if(a==XS&&b==ZS)
   { m=pow(C,2)*pow(DT,2);
     n=DX*DZ;
     s[e]=wavelet(F,N0,e,DT);
     mn=m*s[e]/n;
    }
  else mn=0.0;
 return mn;
}
void printu2()
{
   int i,k;
   ofstream outfile("sjxs.out");
   outfile<<NZ<<endl;
   outfile<<NX<<endl;
   outfile<<DX<<endl;
   outfile<<setiosflags(ios::left)<<setiosflags(ios::scientific);
   for(i=0;i<=NX-1;i++)
   {   for(k=0;k<=NZ-1;k++)
        outfile<<setw(20)<<u2[i][k]<<endl;
        
   }
    cout<<"write the file success!"<<"\n";
    outfile.close();     // 調(diào)用成員函數(shù)close 
}
void main()
{  int i,k,n,g;
   double rx,rz,a,b,w;
   //第一步：速度文件的載入及相關(guān)整理（替換）  
   rx=DT/DX;
   rz=DT/DZ;
   a=pow(C,2)*pow(rx,2);
   b=pow(C,2)*pow(rz,2);


   for(n=0;n<=NT-1;n++)
     {  //第二步：賦初值;初始時(shí)刻的全波場(chǎng)值均為零，P(i, j, 0)=0） 
         //  時(shí)刻dt時(shí)，在炮點(diǎn)S (XS, ZS)給出一個(gè)脈沖震源S(t)，其它點(diǎn)波場(chǎng)P =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;

           }
       //第三步：邊界條件處理及11點(diǎn)差分計(jì)算波場(chǎng)延拓  

       else if(n>=2)
           {  g=n-1;
               for(i=2;i<=NX-3;i++)
               for(k=2;k<=NZ-3;k++)
                 {  w=(2-2.5*a-2.5*b)*u1[i][k]-u0[i][k];
                    if(i==XS&&k==ZS)
                       u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w+shot(i,k,g);
                    else
                       u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w;
                 }
                u2[0][0]=a/12*(16*u1[1][0]-u1[2][0])+b/12*(16*u1[0][1]-u1[0][2])+(2-2.5*a-2.5*b)*u1[0][0]-u0[0][0];

               u2[1][0]=a/12*(16*(u1[0][0]+u1[2][0])-u1[3][0])+b/12*(16*u1[1][1]-u1[1][2])+(2-2.5*a-2.5*b)*u1[1][0]-u0[1][0];

               u2[0][1]=a/12*(16*u1[1][1]-u1[2][1])+b/12*(16*(u1[0][0]+u1[0][2])-u1[0][3])+(2-2.5*a-2.5*b)*u1[0][1]-u0[0][1];

               u2[1][1]=a/12*(16*(u1[0][1]+u1[2][1])-u1[3][1])+b/12*(16*(u1[1][0]+u1[1][2])-u1[1][3])+(2-2.5*a-2.5*b)*u1[1][1]-u0[1][1];

               u2[NX-1][0]=a/12*(16*u1[NX-2][0]-u1[NX-3][0])+b/12*(16*u1[NX-1][1]-u1[NX-1][2])+(2-2.5*a-2.5*b)*u1[NX-1][0]-u0[NX-1][0];

               u2[NX-2][0]=a/12*(16*(u1[NX-3][0]+u1[NX-1][0])-u1[NX-4][0])+b/12*(16*u1[NX-2][1]-u1[NX-2][2])+(2-2.5*a-2.5*b)*u1[NX-2][0]-u0[NX-2][0];

               u2[NX-1][1]=a/12*(16*u1[NX-2][1]-u1[NX-3][1])+b/12*(16*(u1[NX-1][0]+u1[NX-1][2])-u1[NX-1][3])+(2-2.5*a-2.5*b)*u1[NX-1][1]-u0[NX-1][1];

               u2[NX-2][1]=a/12*(16*(u1[NX-3][1]+u1[NX-1][1])-u1[NX-4][1])+b/12*(16*(u1[NX-2][0]+u1[NX-2][2])-u1[NX-2][3])+(2-2.5*a-2.5*b)*u1[NX-2][1]-u0[NX-2][1];

               u2[0][NZ-1]=a/12*(16*u1[1][NZ-1]-u1[2][NZ-1])+b/12*(16*u1[0][NZ-2]-u1[0][NZ-3])+(2-2.5*a-2.5*b)*u1[0][NZ-1]-u0[0][NZ-1];

               u2[0][NZ-2]=a/12*(16*u1[1][NZ-2]-u1[2][NZ-2])+b/12*(16*(u1[0][NZ-3]+u1[0][NZ-1])-u1[0][NZ-4])+(2-2.5*a-2.5*b)*u1[0][NZ-2]-u0[0][NZ-2];

               u2[1][NZ-1]=a/12*(16*(u1[0][NZ-1]+u1[2][NZ-1])-u1[3][NZ-1])+b/12*(16*u1[1][NZ-2]-u1[1][NZ-3])+(2-2.5*a-2.5*b)*u1[1][NZ-1]-u0[1][NZ-1];

               u2[1][NZ-2]=a/12*(16*(u1[0][NZ-2]+u1[2][NZ-2])-u1[3][NZ-2])+b/12*(16*(u1[1][NZ-3]+u1[1][NZ-1])-u1[1][NZ-4])+(2-2.5*a-2.5*b)*u1[1][NZ-2]-u0[1][NZ-2];

               u2[NX-1][NZ-1]=a/12*(16*u1[NX-2][NZ-1]-u1[NX-3][NZ-1])+b/12*(16*u1[NX-1][NZ-2]-u1[NX-1][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-1]-u0[NX-1][NZ-1];

               u2[NX-2][NZ-1]=a/12*(16*(u1[NX-3][NZ-1]+u1[NX-1][NZ-1])-u1[NX-4][NZ-1])+b/12*(16*u1[NX-2][NZ-2]-u1[NX-2][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-1]-u0[NX-2][NZ-1];

               u2[NX-1][NZ-2]=a/12*(16*u1[NX-2][NZ-2]-u1[NX-3][NZ-2])+b/12*(16*(u1[NX-1][NZ-3]+u1[NX-1][NZ-1])-u1[NX-1][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-2]-u0[NX-1][NZ-2];

               u2[NX-2][NZ-2]=a/12*(16*(u1[NX-3][NZ-2]+u1[NX-1][NZ-2])-u1[NX-4][NZ-2])+b/12*(16*(u1[NX-2][NZ-3]+u1[NX-2][NZ-1])-u1[NX-2][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-2]-u0[NX-2][NZ-2];

               for(i=2;i<=NX-3;i++)
                  u2[i][0]=a/12*(16*(u1[i-1][0]+u1[i+1][0])-(u1[i-2][0]+u1[i+2][0]))+b/12*(16*u1[i][1]-u1[i][2])+(2-2.5*a-2.5*b)*u1[i][0]-u0[i][0];

               for(i=2;i<=NX-3;i++)
                  u2[i][NZ-1]=a/12*(16*(u1[i-1][NZ-1]+u1[i+1][NZ-1])-(u1[i-2][NZ-1]+u1[i+2][NZ-1]))+b/12*(16*u1[i][NZ-2]-u1[i][NZ-3])+(2-2.5*a-2.5*b)*u1[i][NZ-1]-u0[i][NZ-1];

               for(k=2;k<=NZ-3;k++)
                  u2[0][k]=a/12*(16*u1[1][k]-u1[2][k])+b/12*(16*(u1[0][k-1]+u1[0][k+1])-(u1[0][k-2]+u1[0][k+2]))+(2-2.5*a-2.5*b)*u1[0][k]-u0[0][k];

               for(k=2;k<=NZ-3;k++)
                  u2[NX-1][k]=a/12*(16*u1[NX-2][k]-u1[NX-3][k])+b/12*(16*(u1[NX-1][k-1]+u1[NX-1][k+1])-(u1[NX-1][k-2]+u1[NX-1][k+2]))+(2-2.5*a-2.5*b)*u1[NX-1][k]-u0[NX-1][k];

               for(i=2;i<=NX-3;i++)
                  u2[i][1]=a/12*(16*(u1[i-1][1]+u1[i+1][1])-(u1[i-2][1]+u1[i+2][1]))+b/12*(16*(u1[i][0]+u1[i][2])-u1[i][3])+(2-2.5*a-2.5*b)*u1[i][1]-u0[i][1];

               for(i=2;i<=NX-3;i++)
                  u2[i][NZ-2]=a/12*(16*(u1[i-1][NZ-2]+u1[i+1][NZ-2])-(u1[i-2][NZ-2]+u1[i+2][NZ-2]))+b/12*(16*(u1[i][NZ-3]+u1[i][NZ-1])-u1[i][NZ-4])+(2-2.5*a-2.5*b)*u1[i][NZ-2]-u0[i][NZ-2];

               for(k=2;k<=NZ-3;k++)
                  u2[1][k]=a/12*(16*(u1[0][k]+u1[2][k])-u1[3][k])+b/12*(16*(u1[1][k-1]+u1[1][k+1])-(u1[1][k-2]+u1[1][k+2]))+(2-2.5*a-2.5*b)*u1[1][k]-u0[1][k];

               for(k=2;k<=NZ-3;k++)
                  u2[NX-2][k]=a/12*(16*(u1[NX-3][k]+u1[NX-1][k])-u1[NX-4][k])+b/12*(16*(u1[NX-2][k-1]+u1[NX-2][k+1])-(u1[NX-2][k-2]+u1[NX-2][k+2]))+(2-2.5*a-2.5*b)*u1[NX-2][k]-u0[NX-2][k];

               /* 左右邊界條件 */
               for(k=0;k<=NZ-1;k++)
                  {  u2[1][k]=u1[1][k]+u1[2][k]-u0[2][k]+C*DT/DX*((u1[2][k]-u1[1][k])-(u0[3][k]-u0[2][k]));
                     u2[0][k]=u1[0][k]+u1[1][k]-u0[1][k]+C*DT/DX*((u1[1][k]-u1[0][k])-(u0[2][k]-u0[1][k]));
                     u2[NZ-2][k]=u1[NZ-2][k]+u1[NZ-3][k]-u0[NZ-3][k]+C*DT/DX*((u1[NZ-3][k]-u1[NZ-2][k])-(u0[NZ-4][k]-u0[NZ-3][k]));
                     u2[NZ-1][k]=u1[NZ-1][k]+u1[NZ-2][k]-u0[NZ-2][k]+C*DT/DX*((u1[NZ-2][k]-u1[NZ-1][k])-(u0[NZ-3][k]-u0[NZ-2][k]));
                  }
               /* 上下邊界條件*/
               for(i=0;i<=NX-1;i++)
                  {  u2[i][1]=u1[i][1]+u1[i][2]-u0[i][2]+C*DT/DZ*((u1[i][2]-u1[i][1])-(u0[i][3]-u0[i][2]));
                     u2[i][0]=u1[i][0]+u1[i][1]-u0[i][1]+C*DT/DZ*((u1[i][1]-u1[i][0])-(u0[i][2]-u0[i][1]));
                     u2[i][NZ-2]=u1[i][NZ-2]+u1[i][NZ-3]-u0[i][NZ-3]+C*DT/DZ*((u1[i][NZ-3]-u1[i][NZ-2])-(u0[i][NZ-4]-u0[i][NZ-3]));
                     u2[i][NZ-1]=u1[i][NZ-1]+u1[i][NZ-2]-u0[i][NZ-2]+C*DT/DZ*((u1[i][NZ-2]-u1[i][NZ-1])-(u0[i][NZ-3]-u0[i][NZ-2]));
                  }
               if(n==T1)
                  printu2();
               // 數(shù)組循環(huá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];
           }

      }

    

}