// Coordinate.cpp: implementation of the CCoordinate class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "BiJin.h"
#include "Coordinate.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CCoordinate::CCoordinate()
{
	m_pPtrList = new CPtrList;
}

CCoordinate::~CCoordinate()
{
	for(int i=0;i<m_pPtrList->GetCount();i++)
	{
		POSITION pos;
		pos = m_pPtrList->FindIndex(i);
		CFoldPoint *pFoldPoint = (CFoldPoint *)m_pPtrList->GetAt(pos);
		m_pPtrList->RemoveAt(pos);
	}
	m_pPtrList->RemoveAll();
	delete m_pPtrList;
}
void CCoordinate::Getcoordinate(CStdioFile *pFile)
{
	int m_Length = pFile->GetLength();
	CFoldPoint *m_pFoldPoint = new CFoldPoint;
	CString pBuf;
	pFile->ReadString(pBuf);
	m_Length -= pBuf.GetLength();
	split(pBuf,m_dest,',');
	m_pFoldPoint->x = atof(m_dest.GetAt(0));
	m_pFoldPoint->y = atof(m_dest.GetAt(1));
	m_pPtrList->AddTail((CObject *)m_pFoldPoint);
	while(m_Length)
	{
		m_pFoldPoint = new CFoldPoint;
		pFile->ReadString(pBuf);
		m_Length -= (pBuf.GetLength()+2);
		split(pBuf,m_dest,',');
		m_pFoldPoint->x = atof(m_dest.GetAt(0));
		m_pFoldPoint->y = atof(m_dest.GetAt(1));
		m_pPtrList->AddTail((CObject *)m_pFoldPoint);
	}
	/*int m_nLength = pFile->GetLength();
	char *m_pBuf = new char[m_nLength];
	pFile->Read(m_pBuf,m_nLength);
	CFoldPoint *m_pFoldPoint;
	int m_nCount = 1;
	for(int i=0;i<m_nLength;i++)
	{
		if(m_pBuf[i] == 13)
			m_nCount++;
	}
	char *p1 = m_pBuf;
	char *p2 = m_pBuf;
	for(int j=0;j<m_nCount-1;j++)
	{
		m_pFoldPoint = new CFoldPoint;
		while((*p2) != ',')
			p2++;
		m_pFoldPoint->x = 0;
		while(p1 < p2)
		{
			m_pFoldPoint->x += (int)((*p1 - 48)*pow(10,p2-p1));
			p1++;
		}
		p1++;
		p2++;
		while((*p2) != 13)
			p2++;
		m_pFoldPoint->y  = 0;
		while(p1 < p2)
		{ 
			m_pFoldPoint->y  += (int)((*p1 - 48)*pow(10,p2-p1-1));
			p1++;
		}
		p1 +=2;
		p2 +=2;
		m_pPtrList->AddTail((CObject *)m_pFoldPoint);
	}
	m_pFoldPoint = new CFoldPoint;
	while((*p2) != ',')
		p2++;
	m_pFoldPoint->x = 0;
	while(p1 < p2)
	{
		m_pFoldPoint->x += (int)((*p1 - 48)*pow(10,p2-p1));
		p1++;
	}
	p1++;
	p2++;
	while((*p2) != '.')
		p2++;
	m_pFoldPoint->y = 0;
	while(p1 < p2)
	{
		m_pFoldPoint->y += (int)((*p1 - 48)*pow(10,p2-p1-1));
		p1++;
	}
	m_pPtrList->AddTail((CObject *)m_pFoldPoint);
	delete []m_pBuf;
	p1 = NULL;
	p2 = NULL;*/
}
void CCoordinate::split(CString source,CStringArray &dest ,char a)
{
	dest.RemoveAll();
	for(int i=0;i<source.GetLength();i++)
	{
		if(source.GetAt(i) == a)
		{
			dest.Add(source.Left(i));//去掉右邊
			for(int j=0;j<dest.GetSize()-1;j++)
			{
				dest[dest.GetSize()-1] = dest[dest.GetSize()-1].Right(dest[dest.GetSize()-1].GetLength() 
					- dest[j].GetLength() - 1);//去掉左邊
			}
		}
	}
	dest.Add(source);
	for(int j=0;j<dest.GetSize()-1;j++)
	{
		dest[dest.GetSize()-1] = dest[dest.GetSize()-1].Right(dest[dest.GetSize()-1].GetLength() 
					- dest[j].GetLength() - 1);
	}
}
bool CCoordinate::GetLowest()
{
	if(!m_pPtrList)return false;
	int size = m_pPtrList->GetCount();
	POSITION lowPos,tmpPos,lowPos1;
	double lowY;
	lowPos1 = lowPos = m_pPtrList->FindIndex(0);
	((CFoldPoint *)m_pPtrList->GetAt(lowPos))->m_bLow = true;
	for(int i=1;i<size-5;i++)
	{
		lowPos = m_pPtrList->FindIndex(i);
		CFoldPoint *pFoldPoint = (CFoldPoint *)m_pPtrList->GetAt(lowPos);
		lowY = pFoldPoint->y;
		for(int j=i;j<i+5;j++)
		{
			tmpPos = m_pPtrList->FindIndex(j);
			pFoldPoint = (CFoldPoint *)m_pPtrList->GetAt(tmpPos);
			if(pFoldPoint->y < lowY)
			{
				lowY = pFoldPoint->y;
				lowPos = tmpPos;
			}			
		}
		if(lowPos != lowPos1)
		{
			((CFoldPoint *)m_pPtrList->GetAt(lowPos))->m_bLow = true;
		}
		lowPos1 = lowPos;
	}
	lowPos = m_pPtrList->FindIndex(size-1);
	((CFoldPoint *)m_pPtrList->GetAt(lowPos))->m_bLow = true;
	return true;
}
CPtrList* CCoordinate::SPLine(CPtrList *pList,int SM)//計算pDestList的每個點的橫，縱坐標
//三次樣條插值算法
{
	CPtrList *pDestList = new CPtrList;
	CFoldPoint *pFoldHead,*pFoldTail;
	POSITION pos;
	CDoubleArray X,Y;
	double XI,YI,XX,YY;
	register long i;
	long RealSM;
	long Bei,Yu;
	CFoldPoint *pFold;
//file://賦初值
	long n=pList->GetCount();
	pos=pList->GetHeadPosition();
	for(i=0;i<n;i++)
	{
		pFold=(CFoldPoint *)pList->GetNext(pos);
		if(pFold->m_bLow == true)
		{
			X.Add(pFold->x);
			Y.Add(pFold->y);
		}
	}
	pFoldHead=(CFoldPoint *)pList->GetHead();
	pFoldTail=(CFoldPoint *)pList->GetTail();
//file://x連續
	YI=0;
	N = X.GetSize();
	RealSM = (N-1)*SM+N;
	for(i=0;i<RealSM;i++)
	{
		Bei=i/(SM+1);
		Yu=i%(SM+1);
		if(Yu!=0)
		{
			XI=X.GetAt(Bei)+(X.GetAt(Bei+1)-X.GetAt(Bei))/(SM+1)*Yu;
			SPLine4(&X,&Y,XI,YI);
			//SPLine1(&X,&Y,XI,YI);
			XX=XI;
			YY=YI;
		}
		else
		{
			XX=X.GetAt(Bei);
			YY=Y.GetAt(Bei);
		}
		pFold=new CFoldPoint;
		pFold->x=XX;
		pFold->y=YY;
		pDestList->AddTail(pFold);
	}
	return pDestList;
}
void CCoordinate::SPLine1(CDoubleArray *X,CDoubleArray *Y,double &XI,double &YI)//用直線擬合曲線，
//不保證光滑
{
	for(int i=1;i<N;i++)
	{
		if(XI>=X->GetAt(i-1) && XI<=X->GetAt(i))
		{
			double A = (X->GetAt(i) - XI)/(X->GetAt(i) - X->GetAt(i-1));
			double B = 1 - A;
			YI = A*Y->GetAt(i-1)+B*Y->GetAt(i); 
		}
	}
}
void CCoordinate::SPLine4(CDoubleArray *X,CDoubleArray *Y,double &XI,double &YI)//算法為三次樣條插值
{
	//給定數組X和Y其中y(i)=x(i),0<x1<x2....<x(n)給指定的數據Xi
	double *A,*B,*C,*G;
	//MX = G中，A為M上副對角線元素，B為M主對角線元素，C為M下副對角線元素，G為方程的常量元素
	register long I;
	double W1,W2,H; 
	A = new double[N];
	B = new double[N];
	C = new double[N];
	G = new double[N];
	for(I=1;I<N;I++)
	{
		B[I]=X->GetAt(I) - X->GetAt(I-1);//B(I)為X(I) - X(I-1)
		C[I]=(Y->GetAt(I) - Y->GetAt(I-1))/B[I];//C(I)為(Y(I)-Y(I-1))/(X(I)-X(I-1))
	}
	for(I=0;I<N;I++)
	{
		A[I]=B[I]+B[I+1];
		G[I]=6*(C[I+1] - C[I])/A[I];
		A[I]=B[I]/A[I];
	}
	for(I=1;I<N;I++)
	{
		C[I]=1-A[I];
		B[I]=2;
	}
	//上述兩個算法為三對角的分解循環，MX = G中，
	//A為M上副對角線元素，B為M主對角線元素，C為M下副對角線元素，G為方程的常量元素
	B[N-1]=X->GetAt(N-1);
	C[0]=1;
	A[N-1]=-1;
	A[0]=0;
	C[N-1]=0;
	ZG(A,B,C,&G);
	//所得到的矩陣方程(x(j)-x(j-1))*y''(j-1)/6+(x(j+1)-x(j-1))*y''(j)/3+(x(j+1)-x(j))*y''(j+1)/6
	//=(y(j+1)-y(j))/(x(j+1)-x(j)) - (y(j)-y(j-1))/(x(j)-x(j-1))
	for(I=1;I<N;I++)
	{
		if(XI>=X->GetAt(I-1) && XI<=X->GetAt(I))//GE LE
		{
			H=X->GetAt(I) - X->GetAt(I-1);
			W1=X->GetAt(I) - XI;
			W2=XI - X->GetAt(I-1);
			YI=W1*W1*W1*G[I-1]/6/H;
			YI=YI+W2*W2*W2*G[I]/6/H;
			YI=YI+W1*(Y->GetAt(I-1)-G[I-1]*H*H/6)/H;
			YI=YI+W2*(Y->GetAt(I)-G[I]*H*H/6)/H;
		}
	}
	delete []A;
	delete []B;
	delete []C;
	delete []G;
}
void CCoordinate::ZG(double *A,double *B,double *C,double **G) //MX = G,其中M
//含有上副對角元素A,主對角元素B,以及下副對角元素C,最終所得的解坐標保存到G中
{
	//追趕法
	register long I;
	C[0]=C[0]/B[0];
	for(I=1;I<N;I++)
	{
		B[I]=B[I]-A[I]*C[I-1];
		C[I]=C[I]/B[I];//分解
	}
	A[0]=0;
	C[N-1]=0;
	(*G)[0]=(*G)[0]/B[0];
	for(I=1;I<N;I++)
	{
		(*G)[I]=((*G)[I]-A[I]*(*G)[I-1])/B[I];//前代
	}
	for(I=N-2;I>-1;I--)//DO 30 I=N-1,1,-1
	{
		(*G)[I]=(*G)[I]-C[I]*(*G)[I+1];//回代
	}
}
double CCoordinate::NewY(CPtrList *pList,double XI)
{
	POSITION pos;
	CDoubleArray X,Y;
	double YI;
	register long i;
	CFoldPoint *pFold;
//file://賦初值
	long n=pList->GetCount();
	pos=pList->GetHeadPosition();
	for(i=0;i<n;i++)
	{
		pFold=(CFoldPoint *)pList->GetNext(pos);
		if(pFold->m_bLow == true)
		{
			X.Add(pFold->x);
			Y.Add(pFold->y);
		}
	}
//file://x連續
	YI=0;
	N = X.GetSize();
	SPLine4(&X,&Y,XI,YI);
	return YI;
}