// 有限差分波動方程不加邊界條件二階變速剖面
   // u2[i][k]=2*(1-a-b)*u1[i][k]+a*(u1[i+1][k]+u1[i-1][k])+b*(u1[i][k+1]+u1[i][k-1])-shot(i,k,g)-u0[i][k];  
   //
   //  eponge邊界吸收函數
   //    xx=0.305/IABMAX;
   //    for(i=0;i<=IABMAX-1;i++)
   //      {  aaa=pow((IABMAX-i),2);
   //         eponge[i]=exp(-xx*aaa);
   //      }
   //
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>       //輸入輸出格式控制
#include <fstream.h>
#define F 35.0            //雷克子波主頻，單位:hz 
#define N0 2            //震源起爆時間 
#define NZ 241             //縱向最大采樣點號 
#define NX 91             //橫向最大采樣點號 
#define DT 0.002          //時間采樣間隔，單位：s 
#define NT 200             //時間最大樣點號 
#define DZ 12             //縱向采樣間隔 
#define DX 12            //橫向采樣間隔 
#define XS 45             //震源位置橫坐標樣點號 
#define ZS 28             //震源位置縱坐標樣點號 
#define C 2000            //波傳播速度，單位：m/s 
#define T1 49          //輸出t=T1*DT時的切片 
#define IABMAX 20           //邊界吸收寬度 
#define N 25

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 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 print()
{
   int i,k;
   ofstream outfile("profile.out");
   outfile<<NT<<endl;
   outfile<<NX-2*IABMAX<<endl;
   outfile<<DX<<endl;
   outfile<<setiosflags(ios::left)<<setiosflags(ios::scientific);
   for(i=IABMAX;i<=NX-IABMAX-1;i++)
   {  for(k=0;k<=NT-1;k++) 
        outfile<<setw(20)<<out[i][k]<<endl;
     
   }
    cout<<"write the file success!"<<"\n";
    outfile.close();     // 調用成員函數close 
}
void main()
{ int i,k,n,g;
  double rx,rz;
  double a[NX][NZ],b[NX][NZ];
  double eponge[IABMAX];
  double aaa,xx;
  rx=DT/DX;
  rz=DT/DZ;
  for(i=0;i<=NX-1;i++)
  for(k=0;k<=NZ-1;k++)
  {   if(k<=40)
         {  a[i][k]=pow(1000.0,2)*pow(rx,2);
            b[i][k]=pow(1000.0,2)*pow(rz,2);
         }
      else if(k<=50)
		  {  a[i][k]=pow(2000.0,2)*pow(rx,2);
             b[i][k]=pow(2000.0,2)*pow(rz,2);
          }
      else  
		  {  a[i][k]=pow(2500.0,2)*pow(rx,2);
             b[i][k]=pow(2500.0,2)*pow(rz,2);
          }
      
  }
  //eponge邊界吸收函數
   xx=0.305/IABMAX;
   for(i=0;i<=IABMAX-1;i++)
      {  aaa=pow((IABMAX-i),4);
         eponge[i]=exp(-xx*aaa)*1e-10;
      }
  for(n=0;n<=NT-1;n++)
  {  //初始條件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;

    }
    //當n>=2時 
    else if(n>=2)
    {  g=n-1;
       u2[0][0]=2*(1-a[0][0]-b[0][0])*u1[0][0]+a[0][0]*u1[1][0]+b[0][0]*u1[0][1]-u0[0][0];
       u2[0][NZ-1]=2*(1-a[0][NZ-1]-b[0][NZ-1])*u1[0][NZ-1]+a[0][NZ-1]*u1[1][NZ-1]+b[0][NZ-1]*u1[0][NZ-2]-u0[0][NZ-1];
       u2[NX-1][0]=2*(1-a[NX-1][0]-b[NX-1][0])*u1[NX-1][0]+a[NX-1][0]*u1[NX-2][0]+b[NX-1][0]*u1[NX-1][1]-u0[NX-1][0];
       u2[NX-1][NZ-1]=2*(1-a[NX-1][NZ-1]-b[NX-1][NZ-1])*u1[NX-1][NZ-1]+a[NX-1][NZ-1]*u1[NX-2][NZ-1]+b[NX-1][NZ-1]*u1[NX-1][NZ-2]-u0[NX-1][NZ-1];
       for(i=1;i<=NX-2;i++)
       { u2[i][0]=2*(1-a[i][0]-b[i][0])*u1[i][0]+a[i][0]*(u1[i+1][0]+u1[i-1][0])+b[i][0]*u1[i][1]-u0[i][0];
         u2[i][NZ-1]=2*(1-a[i][NZ-1]-b[i][NZ-1])*u1[i][NZ-1]+a[i][NZ-1]*(u1[i+1][NZ-1]+u1[i-1][NZ-1])+b[i][NZ-1]*u1[i][NZ-2]-u0[i][NZ-1];
        }
       for(k=1;k<=NZ-2;k++)
       { u2[0][k]=2*(1-a[0][k]-b[0][k])*u1[0][k]+a[0][k]*u1[1][k]+b[0][k]*(u1[0][k+1]+u1[0][k-1])-u0[0][k];
         u2[NX-1][k]=2*(1-a[NX-1][k]-b[NX-1][k])*u1[NX-1][k]+a[NX-1][k]*u1[NX-2][k]+b[NX-1][k]*(u1[NX-1][k+1]+u1[NX-1][k-1])-u0[NX-1][k];
        }
       for(i=1;i<=NX-2;i++)
       for(k=1;k<=NZ-2;k++)
	   {   if(i==XS&&k==ZS)
              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])-shot(i,k,g)-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];
       }

        //左右邊界吸收 
               for(k=0;k<=NZ-1;k++)
               for(i=0;i<=IABMAX-1;i++)
			   {  
				  u2[i][k]=u2[i][k]*eponge[i];
                  u2[i][k]=u2[i][k]*eponge[i];
                  u2[NX-1-i][k]=u2[NX-1-i][k]*eponge[i];
                  u2[NX-1-i][k]=u2[NX-1-i][k]*eponge[i];
			   }
               //上下邊界吸收 
               for(k=0;k<=IABMAX-1;k++)
               for(i=0;i<=NX-1;i++)
               {  
                  u2[i][k]=u2[i][k]*eponge[k];
                  u2[i][k]=u2[i][k]*eponge[k];
                  u2[i][NZ-1-k]=u2[i][NZ-1-k]*eponge[k];
                  u2[i][NZ-1-k]=u2[i][NZ-1-k]*eponge[k];
               } 
        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;
   }
   print();  
}