#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "math.h"
#include "time.h"

#define N 151
#define PI 3.1415926

void Forward(double a[],double sigma0,double u[],double g[],double h);
void wavelet(double g[],double h);
void  main()
{
	int i;
	double a[N-1],sigma[N],u[2*N-1],g[2*N];
	double T=0.3,h=0.002;
	FILE * fp1,* fp2;

//	for(i=0;i<N;i++)
//		scanf("%lf",&sigma[i]);
	for(i=0;i<51;i++)
		sigma[i] = 4.0;
	for(i=51;i<101;i++)
		sigma[i] = 3.0;
	for(i=101;i<151;i++)
		sigma[i] = 5.0;
	for(i=0;i<N-1;i++)
		a[i] = 2*sigma[i+1]/(sigma[i]+sigma[i+1]);
//	h=T/N;
	wavelet(g,h);  //獲取地震子波
	Forward(a,sigma[0],u,g,h);  //u 存放正演計(jì)算結(jié)果
	if((fp1=fopen("F_result.dat","w"))==NULL)
	{
		printf("cannot open this file!\n");
		exit(0);
	}
	for(i=0;i<2*N-1;i++)
		fprintf(fp1,"%f\n",u[i]);
	fclose(fp1);
	if((fp2=fopen("wavelet.dat","w"))==NULL)
	{
		printf("cannot open this file!\n");
		exit(0);
	}
	for(i=0;i<2*N;i++)
		fprintf(fp2,"%f\n",g[i]);
	fclose(fp2);
}

/**********************the beginning of the forward program**************************/

void Forward(double a[],double sigma0,double u[],double g[],double h)  // 2T：為地震記錄長度
{
	int i,j,k,m;	
	double A[N-1];
	double B[N-1][3];        //帶型矩陣的存儲
	double U[(N-1)*(2*N-1)];
	double H[(N-1)*(2*N-1)];

// 計(jì)算矩陣A （對角存儲）
	for(i=0;i<N-1;i++)
		A[i] = -1/(2-a[i]);

// 計(jì)算矩陣B （帶型存儲）
	B[0][0] = 1.0;
	B[0][1] = a[0]/(2-a[0]);
	B[0][2] = 0.0;
	B[N-2][0] = 1.0;
	B[N-2][1] = 0.0;
	B[N-2][2] = 0.0;
	for(i=1;i<N-2;i++)
	{
		B[i][0] = 1.0;
		B[i][1] = 0.0;
		B[i][2] = a[i]/(2-a[i]);
	}

// 計(jì)算矩陣H
	for(j=0;j<(N-1)*(2*N-1);j++)
	{
		if((j%(N-1))==0)
			H[j]=-h*g[j/(N-1)+1]/sigma0;
		else
			H[j]=0.0;
	}

// 由矩陣A，B和矩陣H，計(jì)算矩陣U
	for(i=0;i<N-1;i++)
		U[i]=H[i]/A[i];

	for(i=N-1;i<2*(N-1);i++)
	{
		k=i-(N-1);
		if(k==0)
			U[i]=(H[i]-B[0][0]*U[0]-B[0][1]*U[1])/A[k];
		else if(k==N-2)
			U[i]=(H[i]-B[k][0]*U[k-1])/A[k];
		else
			U[i]=(H[i]-B[k][0]*U[k-1]-B[k][2]*U[k+1])/A[k];
	}
	for(m=0;m<=2*N-4;m++)
	{
		for (i=2*(N-1);i<3*(N-1);i++)
		{		
			k=i-2*(N-1);
			if (k==0)
				U[m*(N-1)+i]=(H[m*(N-1)+i]-A[k]*U[m*(N-1)+i-(2*N-2)]
				-B[k][0]*U[m*(N-1)+i-(N-1)]-B[k][1]*U[m*(N-1)+i-(N-2)])/A[k];
			else if(k==N-2)
				U[m*(N-1)+i]=(H[m*(N-1)+i]-A[k]*U[m*(N-1)+i-(2*N-2)]
				-B[k][0]*U[m*(N-1)+i-N])/A[k];
			else
				U[m*(N-1)+i]=(H[m*(N-1)+i]-A[k]*U[m*(N-1)+i-(2*N-2)]
				-B[k][0]*U[m*(N-1)+i-N]-B[k][2]*U[m*(N-1)+i-(N-2)])/A[k];
		}
	}

	for(i=0,j=0;i<(N-1)*(2*N-1);i++)
	{	
		if(!(i%(N-1)))
		{
			u[j] = U[i];			
			j++;
		}
	}
	printf("%d\n",clock());
  }
/**********************the end of forward*******************************/

/***************************Program of seismic wavelet*****************************/
/****************  取Ricker子波g(t)=(1-2*(pi*f*t)^2)*exp(-(pi*f*t)^2)**************/
/****************  其中 f 為Ricker子波的峰值頻率，取 f = 30 Hz ********************/

void wavelet(double g[],double h)
{
	int j,k;
	float f=40.0,T;
	
	T=1/f;
	k=(int)(T/h);
	for(j=0;j<2*N;j++)
	{
		g[j]=20*(1-2*(PI*PI*f*f*(j-k)*h*(j-k)*h))*exp(-(PI*PI*f*f*(j-k)*h*(j-k)*h));
	}
}