//---------------------------------------------------------------//
//                 AVO 技術(shù)正演模型實(shí)驗(yàn)                          //
//                 3 層模型（一）                                //
//                 AVO 技術(shù)正演模型實(shí)驗(yàn)                          //
//                 40炮24道                                       //
//                 地層傾角為0度                                 //
//---------------------------------------------------------------//
//---------------------------------------------------------------//
#include <stdio.h>
#include<math.h>
#define  PI   3.1415926
#define nc 3             //地層層數(shù)
#define np 4             //地層參數(shù)個(gè)數(shù)
#define  dt   0.002      //采樣間隔 
#define  zbn   51        //子波長度
#define  zf   30         //子波主頻
#define  zA   1.0        //子波最大振幅

//——————————觀測系統(tǒng)參數(shù)————————————//
#define  Nd   24         //迭加次數(shù)
#define  dx   10         //道間距
#define  ds   10         //炮間距//
#define  x_0  20        //最小偏移距
#define  L    48         //接收道數(shù)
#define  Ns   40         //炮數(shù)
#define  tl   300        //采樣點(diǎn)數(shù)
#define  ncdp  30        //cdp道數(shù)

void main()
{
	int   i,j,k;
	float par[nc][np];
	float x[Ns][L],zb[zbn];
	float par0[np]={0.0};
	FILE *fp1;
	FILE *fp2;
	FILE *fp3;
	FILE *fp4;
	FILE *fp5;
    FILE *fp6;
	FILE *fp7;
    //====  子程序說明  ====//
	void fzb(float zb[zbn]);
	void flljl(float x[Ns][L],float zb[zbn],float par[nc][np]);
	printf("程序運(yùn)行請稍后......\n\n");
    //=======================  程序主體  =====================//
	if((fp1=fopen("地層模型參數(shù)表2-1.txt","r"))==NULL)
		printf("地層參數(shù)文件沒有找到！\n");
	else
		printf("參數(shù)文件打開成功！\n");
	fp3=fopen("含油氣層參數(shù).txt","r");
//	for(i=0;i<np;i++)
	//{
		//fscanf(fp7,"%f\n",&par_0[i]);
		fscanf(fp3,"%f\t%f\t%f\t%f\n",&par0[0],&par0[1],&par0[2],&par0[3]);
     	printf("%f\t%f\t%f\t%f\n",par0[0],par0[1],par0[2],par0[3]);
	//}
	for(i=0;i<nc;i++)
	{
		fscanf(fp1,"%f\t%f\t%f\t%f\n",&par[i][0],&par[i][1],&par[i][2],&par[i][3]);
		printf("%f m/s\t%f m/s\t%f g/cm3\t%f m\n",par[i][0],par[i][1],par[i][2],par[i][3]);
		if(i>=0)
			printf("-R-----------------------------------------------------------R-\n");
	}
    //========  子波  =====//
	fzb(zb);
	fp2=fopen("子波.txt","w+");
	for(i=0;i<zbn;i++)
		fprintf(fp2,"%f\t%f\n",dt*i,zb[i]);
    //====  計(jì)算各道反射系數(shù)，形成理論記錄，抽道集  ====//
	flljl(x,zb,par);
	fclose(fp1);
	fclose(fp2);
	fclose(fp3);
}
//=======================  子波子程序  =====================//
void fzb(float zb[zbn])
{
	int i,j;
	float t;
	for(i=0;i<=zbn/2;i++)
	{
		t=dt*i;
		zb[zbn/2+i]=zA*cos(2*PI*zf*t)*exp(-0.5*pow(PI*zf*t,2));
		zb[zbn/2-i]=zA*cos(2*PI*zf*t)*exp(-0.5*pow(PI*zf*t,2));
	}
}

//=======================  計(jì)算各道反射系數(shù)，形成理論記錄，抽道集子程序  =====================//
void flljl(float x[Ns][L],float zb[zbn],float par[nc][np])
{
	int i,i1,i2,j,k,kn,r1,r2,l2,flag;
	int j1;
	float y[tl+zbn-1],yy[tl];
	float t_0[nc-1]={0.0},t[nc-1],vr,sum=0.0;
	float R1[nc-1],x1,sa2[nc-1],ta2,vp[2],vs[2],ro[2],K,da,a,A,B,C;//sina2->sa2和一些過渡參數(shù)
	float R[tl]={0.0};
	float par_0[np]={1896.0,50.0,1.50,100.0};//含油氣層參數(shù)
	int X,k2;
	float sa2_j,sa_j,ti0,t00,vi;
	float cdp_1[tl][ncdp]={0.0},cdp_2[tl][ncdp]={0.0},cdp_3[tl][ncdp]={0.0},tp1[tl]={0.0};//角度疊加約束
	float cdp[tl][ncdp]={0.0},cdp0[tl][ncdp*Nd],tp[tl];
	double zx[tl],zc=2045.0,zza=-0.06 ,zzb=0.06,zN[tl];
	FILE *fp1;
	FILE *fp2;
	FILE *fp3;
	FILE *fp4;
	FILE *fp5;
	FILE *fp6;
	FILE *fp7;
	FILE *fp8;
	FILE *fp9;
    //===== 計(jì)算基本參數(shù) =====//
	fp1=fopen("反射系數(shù).txt","w+");
	fp2=fopen("地震記錄.dat","wb+");
	fp3=fopen("地震記錄(文本).dat","w+");
	fp4=fopen("抽道得的道集.dat","wb+");
	fp5=fopen("抽道疊加后30個(gè)cdp點(diǎn).dat","wb+");
	fp6=fopen("抽道得的道集（文本）.dat","w+");
	fp7=fopen("含油氣層參數(shù).txt","r");
	fp8=fopen("角度疊加(3組角度范圍＋全角度疊加共4個(gè)30CDP).dat","wb+");
	fp9=fopen("查_角度疊加(3組角度范圍＋全角度疊加共4個(gè)30CDP).dat","w+");
	//===== 讀取含油氣參數(shù) =====//
	for(i=0;i<np;i++)
	{
		fscanf(fp7,"%f\n",&par_0[i]);printf("%f\n",par_0[i]);
	}
	//fscanf(fp7,"%f\t%f\t%f\t%f ",&par_0[0],&par_0[1],&par_0[2],&par_0[3]);
	//——————————噪音———————————//
	zx[0]=8388607.0/8388608.0;
	for(i=1;i<tl;i++)
		zx[i]=zc*zx[i-1]-int(zc*zx[i-1]);
    for(i=0;i<tl;i++)
		zN[i]=float(zza+(zzb-zza)*zx[i]);
	//——————————噪音———————————//
	for(i=0;i<Ns;i++)
		for(j=0;j<L;j++)
		{
			x1=x_0+dx*j;
			X=ds*i+x1;
			if(X/dx>=30&&X/dx<40)//含油氣參數(shù)替換
				for(k2=0;k2<np;k2++)
				{
					par[1][k2]=par_0[k2];
					if(i==0)printf("****%d\t%f\n",j,par_0[k2]);
				}
			for(i1=0;i1<nc-1;i1++)
			{
				t_0[i1]=0.0;//注意清零?。?			}
			sum=0.0;
			for(i1=0;i1<nc-1;i1++)
			{
	        	t[i1]=par[i1][3]/par[i1][0];
	        	for(i2=0;i2<=i1;i2++)
					t_0[i1]+=2*(t[i2]);
	        	sum+=t[i1]*par[i1][0]*par[i1][0];
				if(i==0&&j==0)printf("i=%d\tt[i1]=%f\tt_0[i1]=%f\n",i,t[i1],t_0[i1]);
			}
        	vr=sqrt((sum*2)/t_0[nc-2]);
			//if(i==0)printf("%d\t%f\n",j,t_0[nc-2]);
			for(i1=0;i1<tl+zbn-1;i1++)//每道地震記錄前歸零！
				y[i1]=0.0;
			for(i1=0;i1<tl;i1++)
				R[i1]=0.0;
			for(k=0;k<nc-1;k++)
			{
				sa2[k]=pow(x1*par[k][0],2)/pow(vr,2)/(pow(vr*t_0[k],2)+pow(x1,2));
				//sa2[k]=pow(par[k][0],2)*(x1-vr*t_0[k]*atan(x1/(vr*t_0[k])))/(x1*pow(vr,2));
				//if(i==0&&k==1)printf("%d\t%f\n",j,vr);
				if(sa2[k]<0)printf("%d\t,%d\t,%d\tsa2[]=%f\n",i,j,k,sa2[k]);
				ta2=sa2[k]/(1-sa2[k]);
				vp[0]=(par[k+1][0]+par[k][0])/2.0;
				vp[1]=par[k+1][0]-par[k][0];
				vs[0]=(par[k+1][1]+par[k][1])/2.0;
				vs[1]=par[k+1][1]-par[k][1];
				ro[0]=(par[k+1][2]+par[k][2])/2.0;
				ro[1]=par[k+1][2]-par[k][2];
				K=(pow(par[k+1][1],2)/pow(par[k+1][0],2)+pow(par[k][1],2)/pow(par[k][0],2));//K=vs2/vp2
				//過渡參數(shù)計(jì)算↑↓……
				A=(vp[1]/vp[0]+ro[1]/ro[0])/2.0;
				B=vp[1]/(vp[0]*2.0)-4*K*vs[1]/vs[0]-2*K*ro[1]/ro[0];
				C=vp[1]/(2*vp[0]);
				R1[k]=A+B*sa2[k]+C*sa2[k]*ta2;
				//printf("%d,%d,%d\n",i,j,k);
				kn=t_0[k]/dt;
				R[kn]=R1[k];
				
			}
			if(i==0)
			{
				fprintf(fp1,"%f\t%f\t%f\t%f\n",asin(sqrt(sa2[0]))*180/PI,R1[0],asin(sqrt(sa2[1]))*180/PI,R1[1]);
			//	printf("%d\t%f\t%f\tcha=%f\n",j,asin(sqrt(sa2[0]))*180/PI,R1[0],par[0][0]);
			}
            //查??！反射系數(shù)計(jì)算
			for(i1=0;i1<tl;i1++)
				for(j1=0;j1<zbn;j1++)
				{
					y[i1+j1]+=R[i1]*zb[j1];
				}
			for(i1=0;i1<tl;i1++)
			{
				yy[i1]=y[i1+zbn/2]+zN[i1];//
				fprintf(fp3,"%f\t",yy[i1]+0.005*i1);
			}
			fwrite(yy,4,tl,fp2);
			fprintf(fp3,"\n");
			//———  抽道、迭加求和↓  ———//
			for(r1=0;r1<ncdp;r1++)//選道是flag=1,否flag=0
			{
				flag=r1/2;
				if(i>=flag)
					if(j+L/Nd*i-L+L/Nd==r1)//?
					{
						l2=r1*Nd+i-flag;//-r1%Nd;//因該是除以疊加次數(shù)的余數(shù)??！
					    if(r1+1==1)printf("CDP道：i=%d\tj=%d\tcdp=%d\tl2=%d\n",i+1,j+1,r1+1,l2+1);		
						//if(r1==0)printf("X=%d\t%d\n",X,X/dx);
						for(r2=0;r2<tl;r2++) 
						{
					    	cdp0[r2][l2]=yy[r2];
							//———  角度疊加～～確定角度范圍↓...———//
							t00=r2*dt;
							for(k=0;k<nc-1;k++)
							{
								if(t00<=t_0[k])
								{
									vi=par[k][0];
									break;
								}
								else
									vi=par[nc-2][0];
							}
							sa2_j=pow(vi,2)*(x1-vr*t00*atan(x1/(vr*t00)))/(x1*pow(vr,2));
							sa_j=asin(sqrt(sa2_j))*180/PI;
							//三組角度范圍角度疊加//
							if(sa_j>=0&&sa_j<4)
								cdp_1[r2][r1]+=yy[r2];
							else
								cdp_1[r2][r1]+=0.0;
							if(sa_j>=4&&sa_j<8)
								cdp_2[r2][r1]+=yy[r2];
							else
								cdp_2[r2][r1]+=0.0;
						    if(sa_j>=8&&sa_j<35)
								cdp_3[r2][r1]+=yy[r2];
							else
								cdp_3[r2][r1]+=0.0;
						    //———  角度疊加～～確定角度范圍↑...———//
					    	cdp[r2][r1]+=yy[r2];						
						}
					}
			}
	        //———  抽道、迭加求和↑  ———//
		}//6炮24道循環(huán)結(jié)束
		printf("X總測線長度：%d\t總炮數(shù)：%d\n",X,X/ds);
	/*for(j=0;j<ncdp;j++) 
		for(i=0;i<tl;i++) 
		{
			cdp[i][j]=cdp[i][j]/Nd;//振幅平均
			cdp_1[i][j]=cdp[i][j]/Nd;
			cdp_2[i][j]=cdp[i][j]/Nd;
			cdp_3[i][j]=cdp[i][j]/Nd;
		}*/
/*	for(j=0;j<ncdp*Nd;j++)
	{
		for(i=0;i<tl;i++)
			tp[i]=cdp0[i][j];
		fwrite(tp,4,tl,fp4);   
	}*/
	for(j=0;j<ncdp*Nd;j++)
	{
		for(i=0;i<tl;i++)
		{
			tp[i]=cdp0[i][j];
			fprintf(fp6,"%f\t",tp[i]+0.005*i);
		}
		fprintf(fp6,"\n");
		fwrite(tp,4,tl,fp4);   
	}
	for(j=0;j<ncdp;j++) 
	{
		for(i=0;i<tl;i++)
		{
			tp[i]=cdp[i][j];
		}
		fwrite(tp,4,tl,fp5);
	}
	//角度疊加寫入數(shù)據(jù)
	flag=1;
    while(flag<5)
	{
		for(j=0;j<ncdp;j++) 
		{
			for(i=0;i<tl;i++)
			{
				if(flag==1) { tp[i]=cdp_1[i][j];fprintf(fp9,"%f\t",tp[i]+0.005*i);}
				if(flag==2) tp[i]=cdp_2[i][j];
				if(flag==3) tp[i]=cdp_3[i][j];
				if(flag==4) tp[i]=cdp[i][j];
			}
			if(flag==1)fprintf(fp9,"\n");
			fwrite(tp,4,tl,fp8);
		}
		fwrite(tp1,4,tl,fp8);
		printf("%d\n",flag);
		flag++;
	}


	fclose(fp1);
	fclose(fp2);
	fclose(fp3);
	fclose(fp4);
	fclose(fp5);
	fclose(fp6);
	fclose(fp7);
	fclose(fp8);
	fclose(fp9);
}