#include<iostream.h>
#include<fstream.h>
#include<stdio.h>
#include<math.h>
#include<malloc.h>
#define Length 22
#define F_ROW_NUM 30
#define F_COLUMN_NUM 15

double getFalpha(int f1,int f2);
int analyse(double *x);
void getPerWave(int T,double *PerWave,double *x);
double Fvalue[F_ROW_NUM][F_ROW_NUM];
int Frow[F_ROW_NUM],Fcolumn[F_COLUMN_NUM];
double newX[Length];
double PerWave1[Length],PerWave2[Length];
double year[Length],x[Length];
int main()
{
	
	//讀入資料
	ifstream fin;
	fin.open("DATA.txt");
	int i,j;
	for(i=0;i<Length;i++)
	{
		fin>>year[i]>>x[i];
	}
	fin.close();
	//讀入F分布表
	fin.open("F分布值.txt");
	for(i=0;i<F_COLUMN_NUM+1;i++)
	{
		if(i!=0)
		{
			fin>>Fcolumn[i-1];
		}
		else
		{
			int t;
			fin>>t;
		}
	}
	for(i=0;i<F_ROW_NUM;i++)
	{
		for(j=0;j<F_COLUMN_NUM+1;j++)
		{
			if(j==0)
			{
				fin>>Frow[i];
			}
			else
			{
				fin>>Fvalue[i][j-1];
			}
		}
	}
	fin.close();
	int T1 = analyse(x);
	getPerWave(T1,PerWave1,x);
	for(i=0;i<Length;i++)
	{
		newX[i]=x[i]-PerWave1[i];
	}
	int T2=analyse(newX);
	getPerWave(T2,PerWave2,newX);
	ofstream fout;
	fout.open("output.xls",ios::ate|ios::out);
	fout<<"年份\t年最高水位\t"<<T1<<"年周期波\t新序列\t"<<T2<<"年周期波\t擬合水位\t"<<endl;
	for(i=0;i<Length;i++)
	{
		fout<<year[i]<<"\t"<<x[i]<<"\t"<<PerWave1[i]<<"\t"<<newX[i]<<"\t"<<PerWave2[i]<<"\t"<<PerWave1[i]+PerWave2[i]<<endl;
	}
	fout.close();
	cout<<"結果輸出在output.xls"<<endl;
	return 1;
}
double getFalpha(int f1,int f2)//未完善但可以用
{

	if(f1<0||f2<0)
	{
		cout<<"getFalpha Error"<<endl;
		return 0;
	}
	int i,f1Num,f2Num;
	int f1NumIsGiven=0,f2NumIsGiven=0;
	for(i=0;i<F_COLUMN_NUM;i++)
	{
		if(Fcolumn[i]==f1)
		{
			f1Num=i;
			f1NumIsGiven=1;
		}
	}
	for(i=0;i<F_ROW_NUM;i++)
	{
		if(Frow[i]==f2)
		{
			f2Num=i;
			f2NumIsGiven=1;
		}
	}
	if(f1NumIsGiven==1&&f2NumIsGiven==1)
	{
		return Fvalue[f2Num][f1Num];
	}
}

int analyse(double *x)
{
	double *Falpha;
	double *F,*f1,*f2;
	int i,j,k,maxPeriod;
	double totalValue=0,averOfAll;
	for(i=0;i<Length;i++)
	{
		totalValue+=x[i];
	}
	averOfAll=totalValue/Length;
	maxPeriod=Length/2;
	double *S1,*S2;
	S1=new double[maxPeriod-1];
	S2=new double[maxPeriod-1];
	f1=new double[maxPeriod-1];
	f2=new double[maxPeriod-1];
	F=new double[maxPeriod-1];
	Falpha=new double[maxPeriod-1]; 
	//遂周期計算
	for(i=2;i<=maxPeriod;i++)
	{
		S1[i-2]=0;S2[i-2]=0;
		int b=i;
		f1[i-2]=b-1;
		f2[i-2]=Length-b;
		int *valueNum=new int[b];
		if(Length%b==0)
		{
			for(j=0;j<b;j++)
			{
				valueNum[j]=Length/b;
			}
		}
		else
		{
			int n=Length%b;
			for(j=0;j<b;j++)
			{
				if(j<n)
				{
					valueNum[j]=Length/b+1;
				}
				else
				{
					valueNum[j]=Length/b;
				}
			}
		}
		//計算方差
		
		double *average=new double [b];
		//組平均值
		for(j=0;j<b;j++)
		{
			double sum=0;
			for(k=0;k<valueNum[j];k++)
			{
				sum+=x[j+k*b];
			}
			average[j]=sum/valueNum[j];
		}
		//S2
		for(j=0;j<b;j++)
		{
			double sum=0;
			for(k=0;k<valueNum[j];k++)
			{
				sum+=pow(x[j+k*b]-average[j],2);
			}
			S2[i-2]+=sum;
		}
		//S1
		for(j=0;j<b;j++)
		{
			double sum=0;
			for(k=0;k<valueNum[j];k++)
			{
				sum+=pow(average[j]-averOfAll,2);
			}
			S1[i-2]+=sum;
		}
		//F
		F[i-2]=(S1[i-2]/f1[i-2])/(S2[i-2]/f2[i-2]);
		Falpha[i-2]=getFalpha(f1[i-2],f2[i-2]);
		delete []valueNum;
		delete []average;
	}
	//輸出個表
	ofstream fout;
	fout.open("output.xls",ios::ate|ios::out);
	fout<<"試驗周期\tS1\tf1\tS2\tf2\tF\tF0.05"<<endl;
	for(i=2;i<=maxPeriod;i++)
	{
		fout<<i<<"\t"<<S1[i-2]<<"\t"<<f1[i-2]<<"\t"<<S2[i-2]<<"\t"<<f2[i-2]<<"\t"<<F[i-2]<<"\t"<<Falpha[i-2]<<endl;
	}
	//選周期
	double *judge=new double[maxPeriod-1];
	int T=2;//所選周期
	for(i=2;i<=maxPeriod;i++)
	{
		judge[i-2]=F[i-2]-Falpha[i-2];
		if(i>2)
		{
			if(judge[i-2]>judge[T-2])
			{
				T=i;
			}
		}
	}
	if(judge[T-2]<0)
	{
		cout<<"周期不存在"<<endl;
	}
	else
	{
		fout<<"T="<<T<<endl;
	}
	fout.close();
	delete []F;
	delete []S1;
	delete []S2;
	delete []f1;
	delete []f2;
	delete []judge;
	delete []Falpha;
	return T;
}

void getPerWave(int T,double *PerWave,double *x)
{
	int b=T;
	int i,j,k,maxPeriod;
	int *valueNum=new int[b];
	if(Length%b==0)
	{
		for(j=0;j<b;j++)
		{
			valueNum[j]=Length/b;
		}
	}
	else
	{
		int n=Length%b;
		for(j=0;j<b;j++)
		{
			if(j<n)
			{
				valueNum[j]=Length/b+1;
			}
			else
			{
				valueNum[j]=Length/b;
			}
		}
	}
		//組平均值
	for(j=0;j<b;j++)
	{
		double sum=0;
		for(k=0;k<valueNum[j];k++)
		{
			sum+=x[j+k*b];
		}
		PerWave[j]=sum/valueNum[j];
	}
	int t=0;
	for(i=0;i<Length;i++)
	{	
		if(i>0&&i%T==0)
		{
			t++;
		}
		PerWave[i]=PerWave[i-T*t];
	}
	delete []valueNum;
}