#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define mm 200  //原數據的節點數
#define pi 3.1415926
#define M 2.0  //M為最大波峰與最大波谷絕對值的比值，取2.0或1.5
#define L 31   // 雷克子波長度
#define t0 0.002 // t0為采樣間隔
#define f0 35 // f0為最小相位的雷克子波的主頻
#define N 1000

void convolution(double *y, double *x, int nx, double *h, int ny);

void main()
{
	int i,shicha[mm],j;
	double shendu[mm];
	double v[mm],d[mm],t[mm-1],b[L],y[N],a;
	double r[N]; //存儲層反射系數值
	int n[mm-1]; //存儲第i個層的反射系數序號
    FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7,*fp8;
	fp1=fopen("shicha.txt","r");
	for(i=0;i<mm;i++)
		fscanf(fp1,"%d",&shicha[i]);
	fclose(fp1);
    fp2=fopen("shendu.txt","r");
	for(i=0;i<mm;i++)
	{	fscanf(fp2,"%lf",&shendu[i]);
	}
	fclose(fp2);
	fp3=fopen("v.dat","w");
    for(i=0;i<mm;i++)
	{
		v[i]=(1e+6)/shicha[i]; 
        fprintf(fp3,"%lf\n",v[i]);
	}
    fp4=fopen("d.txt","w");
	for(i=0;i<mm;i++)
	{	d[i]=1.62+0.00021*v[i]; 
        fprintf(fp4,"%lf\n",d[i]);
	}
	fp5=fopen("t.dat","w");
	t[0]=2*shendu[0]/v[0];
	for(i=1;i<mm;i++)
	{	t[i]=t[i-1]+2*(shendu[i]-shendu[i-1])/v[i];}
	for(i=0;i<mm;i++)
	fprintf(fp5,"%lf\n",t[i]);
    for(i=0;i<N;i++)
		r[i]=0.0;
	for(i=0;i<mm-1;i++)
	{
		n[i]=(int)(t[i]/t0); //計算第i個層的反射系數序號
		 // printf("%d\n",n[i]);
		 j=n[i];
		r[j]=(d[i+1]*v[i+1]-d[i]*v[i])/(d[i+1]*v[i+1]+d[i]*v[i]);//計算層反射系數值
	}
	fp6=fopen("r.dat","w");
    for(i=0;i<N;i++)	
		fprintf(fp6,"%lf\n",r[i]);
	fp7=fopen("子波.dat","w");
    a=(2.0/3.0)*(f0)*(f0)*log10((double)M);// 得到離散的雷克子波
	for(i=0;i<L-1;i++)
	{
		b[i]=sin(2*pi*f0*i*t0)*exp((-1)*a*(i*t0)*(i*t0));
		fprintf(fp7,"%lf\n",b[i]);
	}
	fclose(fp7);
   convolution(y,b,L,r,N);
    fp8=fopen("卷積.dat","w");
    for(i=0;i<N+L;i++) 
    fprintf(fp8, "%lf\n",y[i]);
 }

// 卷積函數
void convolution(double *y, double *x, int nx, double *h, int ny)
{
	int i,j,s;
	for(i=0;i<nx+ny-1;i++)
	{
		y[i]=0;
		for(j=0;j<ny;j++)
		{
			s=i-j;
			if(s>=0 && s<nx)
				y[i]=y[i]+x[s]*h[j];
		}
	}
}