#include "PowerFlow_Type.h"
#include "fstream.h"
#include "math.h"
#include "iostream.h"
#include "iomanip.h"



////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////獲取數(shù)據(jù)/////////////////////////////////////////////////////////////////////

F_GetData_Result_Type PowerFlow_Type::F_GetData()
{

	int Nn,Ng,Nb;
	double V0,Theta0,epsilon;
	Node_Type *Node;
	GeneratorNode_Type *GeneratorNode;
	Branch_Type *Branch;

	int i;

	ifstream fin("source2.txt");
	fin>>Nn>>Ng>>Nb>>V0>>Theta0>>epsilon;
	Node=new Node_Type[Nn+1];
	GeneratorNode=new GeneratorNode_Type[Ng+1];
	Branch=new Branch_Type[Nb+1];

	for(i=1;i<=Nn;i++)
	{
		fin>>Node[i].flag>>Node[i].n>>Node[i].P>>Node[i].Q>>Node[i].V>>Node[i].Theta;
	}
	for(i=1;i<=Ng;i++)
	{
		fin>>GeneratorNode[i].n>>GeneratorNode[i].P>>GeneratorNode[i].Q;
	}
	for(i=1;i<=Nb;i++)
	{
		fin>>Branch[i].n1>>Branch[i].n2>>Branch[i].R>>Branch[i].X>>Branch[i].YK;
	}


	F_GetData_Result_Type F_GetData_Result;
	{
		F_GetData_Result.Nn=Nn;
		F_GetData_Result.Ng=Ng;
		F_GetData_Result.Nb=Nb;
		F_GetData_Result.V0=V0;
		F_GetData_Result.Theta0=Theta0;
		F_GetData_Result.epsilon=epsilon;
		F_GetData_Result.Node=Node;
		F_GetData_Result.GeneratorNode=GeneratorNode;
		F_GetData_Result.Branch=Branch;
	}
	return F_GetData_Result;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////統(tǒng)計(jì)PV節(jié)點(diǎn)和平衡節(jié)點(diǎn)/////////////////////////////////////////////////////////////////

F_FormPVandBalanceNode_Result_Type PowerFlow_Type::F_FormPVandBalanceNode(F_FormPVandBalanceNode_Source_Type F_FormPVandBalanceNode_Source)
{
	int Nn=F_FormPVandBalanceNode_Source.Nn;
	Node_Type *Node=F_FormPVandBalanceNode_Source.Node;

	int Npv,Np;
	PVNode_Type *PVNode;
	BalanceNode_Type *BalanceNode;


	int i,j;

	struct PVNodeArray_Type 
	{
		int flag;
		double V;
	}*PVNodeArray;
	PVNodeArray=new PVNodeArray_Type[Nn+1];

	for (i=1;i<=Nn;i++)
	{
		PVNodeArray[i].flag=0;
	}

	for (i=1;i<=Nn;i++)
	{
		if (Node[i].flag==1)
		{
			PVNodeArray[Node[i].n].flag=1;
			PVNodeArray[Node[i].n].V=Node[i].V;
		}
	}

	Npv=0;
	for (i=1;i<=Nn;i++)
	{
		if (PVNodeArray[i].flag==1)
		{
			Npv=Npv+1;
		}
	}
	PVNode=new PVNode_Type[Npv+1];
	j=1;
	for (i=1;i<=Nn;i++)
	{
		if (PVNodeArray[i].flag==1)
		{
			PVNode[j].n=i;
			PVNode[j].V=PVNodeArray[i].V;
			j=j+1;
		}
	}
	delete PVNodeArray;


	struct BalanceNodeArray_Type 
	{
		int flag;
		double V,Theta;
	}*BalanceNodeArray;
	BalanceNodeArray=new BalanceNodeArray_Type[Nn+1];

	for (i=1;i<=Nn;i++)
	{
		BalanceNodeArray[i].flag=0;
	}

	for (i=1;i<=Nn;i++)
	{
		if (Node[i].flag==2)
		{
			BalanceNodeArray[Node[i].n].flag=1;
			BalanceNodeArray[Node[i].n].V=Node[i].V;
			BalanceNodeArray[Node[i].n].Theta=Node[i].Theta;
			
		}
	}
	
	Np=0;
	for (i=1;i<=Nn;i++)
	{
		if (BalanceNodeArray[i].flag==1)
		{
			Np=Np+1;
		}
	}
	BalanceNode=new BalanceNode_Type[Np+1];
	j=1;
	for (i=1;i<=Nn;i++)
	{
		if (BalanceNodeArray[i].flag==1)
		{
			BalanceNode[j].n=i;
			BalanceNode[j].V=BalanceNodeArray[i].V;
			BalanceNode[j].Theta=BalanceNodeArray[i].Theta;
			j=j+1;
		}
	}
	delete BalanceNodeArray;


	F_FormPVandBalanceNode_Result_Type F_FormPVandBalanceNode_Result;
	{
		F_FormPVandBalanceNode_Result.Npv=Npv;
		F_FormPVandBalanceNode_Result.PVNode=PVNode;
		F_FormPVandBalanceNode_Result.Np=Np;
		F_FormPVandBalanceNode_Result.BalanceNode=BalanceNode;
	}
	return F_FormPVandBalanceNode_Result;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////支路排序/////////////////////////////////////////////////////////////////////
void PowerFlow_Type::F_BranchPaiXu(F_BranchPaiXu_Source_Type F_BranchPaiXu_Source)
{
	int Nb=F_BranchPaiXu_Source.Nb;
	Branch_Type *Branch=F_BranchPaiXu_Source.Branch;
	
	int i,j;
	
	for (i=1;i<=Nb;i++)//同一支路的兩個(gè)節(jié)點(diǎn)前小后大排序
	{
		if (abs(Branch[i].n1)>abs(Branch[i].n2))
		{
			int temp;
			temp=Branch[i].n1;
			Branch[i].n1=Branch[i].n2;
			Branch[i].n2=temp;
		}
	}
	
	for (i=1;i<=Nb-1;i++)//支路排序
	{
		for (j=i+1;j<=Nb;j++)
		{
			if (abs(Branch[i].n1)>abs(Branch[j].n1))
			{
				Branch_Type branchTemp;
				branchTemp=Branch[i];
				Branch[i]=Branch[j];
				Branch[j]=branchTemp;
			}
			else if (abs(Branch[i].n1)==abs(Branch[j].n1))
			{
				if (abs(Branch[i].n2)>abs(Branch[j].n2))
				{
					Branch_Type branchTemp;
					branchTemp=Branch[i];
					Branch[i]=Branch[j];
					Branch[j]=branchTemp;
				}
			}
			else 
				continue;
		}
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////節(jié)點(diǎn)排序/////////////////////////////////////////////////////////////////////

void PowerFlow_Type::F_NodePaixu(F_NodePaiXu_Source_Type F_NodePaiXu_Source)
{
	int Nn=F_NodePaiXu_Source.Nn,Ng=F_NodePaiXu_Source.Ng;
	Node_Type *Node=F_NodePaiXu_Source.Node;
	GeneratorNode_Type *GeneratorNode=F_NodePaiXu_Source.GeneratorNode;

	int i,j;
	for (i=1;i<=Nn-1;i++)
	{
		for (j=i+1;j<=Nn;j++)
		{
			if (Node[i].n>Node[j].n)
			{
				Node_Type Temp;
				Temp=Node[i];
				Node[i]=Node[j];
				Node[j]=Temp;
			}
		}
	}

	for (i=1;i<=Ng-1;i++)
	{
		for (j=i+1;j<=Ng;j++)
		{
			if (GeneratorNode[i].n>GeneratorNode[j].n)
			{
				GeneratorNode_Type Temp;
				Temp=GeneratorNode[i];
				GeneratorNode[i]=GeneratorNode[j];
				GeneratorNode[j]=Temp;
			}
		}
	}
}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////形成導(dǎo)納矩陣第一步,不考慮接地支路的影響///////////////////////////////////////////////////////
F_YForm1_Result_Type PowerFlow_Type::F_YForm1(F_YForm1_Source_Type F_YForm1_Source)
{

	Branch_Type *Branch=F_YForm1_Source.Branch;
	int Nb=F_YForm1_Source.Nb;
	int Nn=F_YForm1_Source.Nn;

	Yii_Type *Yii=new Yii_Type[Nn+1];
	Yii_Type *Yii1=new Yii_Type[Nn+1];
	Yij_Type *Yij=new Yij_Type[Nb+1];
	Yij_Type *Yij1=new Yij_Type[Nb+1];
	int *NYseq=new int[Nn+2];
	int *NYsum=new int[Nn+1];


	int n;
	int i,j;
	double R,X,YK,Zmag2,Gij,Bij,b_ij;
	
	for (i=0;i<=Nn;i++)
	{
		Yii[i].G=Yii[i].B=Yii1[i].G=Yii1[i].B=0;
		NYsum[i]=0;
	}
	for (n=1;n<=Nb;n++)
	{
		i=abs(Branch[n].n1);
		j=abs(Branch[n].n2);
		R=Branch[n].R;
		X=Branch[n].X;
		YK=Branch[n].YK;
		
		Zmag2=R*R+X*X;
		Gij=R/Zmag2;
		Bij=-X/Zmag2;
		
		b_ij=-1.0/X;
		
		if (Branch[n].n1<0||Branch[n].n2<0)
		{
			Yij[n].G=-Gij/YK;
			Yij[n].B=-Bij/YK;
			
			Yii[i].G+=Gij/YK;
			Yii[i].B+=Bij/YK;
			Yii[j].G+=Gij/YK;
			Yii[j].B+=Bij/YK;
			
			Yij1[n].G=0;
			Yij1[n].B=-b_ij/YK;
			
			Yii1[i].B+=b_ij/YK;
			Yii1[j].B+=b_ij/YK;
		}
		else
		{
			Yij[n].G=-Gij;
			Yij[n].B=-Bij;
			
			Yii[i].G+=Gij;
			Yii[i].B+=Bij;
			Yii[j].G+=Gij;
			Yii[j].B+=Bij;
			
			Yij1[n].G=0;
			Yij1[n].B=-b_ij;
			
			Yii1[i].B+=b_ij;
			Yii1[j].B+=b_ij;
		}
		
		Yij[n].n2=Yij1[n].n2=j;
		NYsum[i]+=1;
	}
	NYseq[1]=1;
	for (i=2;i<=Nn;i++)
	{
		NYseq[i]=NYseq[i-1]+NYsum[i-1];
	}
	NYseq[Nn+1]=NYseq[Nn];


	F_YForm1_Result_Type F_YForm1_Result;
	{
		F_YForm1_Result.Yii=Yii;
		F_YForm1_Result.Yii1=Yii1;
		F_YForm1_Result.Yij=Yij;
		F_YForm1_Result.Yij1=Yij1;
		F_YForm1_Result.NYsum=NYsum;
		F_YForm1_Result.NYseq=NYseq;
	}
	return F_YForm1_Result;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////形成導(dǎo)納矩陣第二步,追加接地支路//////////////////////////////////////////////////////////
void PowerFlow_Type::F_YForm2(F_YForm2_Source_Type F_YForm2_Source)
{

	int Nn=F_YForm2_Source.Nn;
	int Nb=F_YForm2_Source.Nb;
	Branch_Type *Branch=F_YForm2_Source.Branch;
	Yii_Type *Yii=F_YForm2_Source.Yii;
	Yii_Type *Yii1=F_YForm2_Source.Yii1;
	Yij_Type *Yij=F_YForm2_Source.Yij;
	Yij_Type *Yij1=F_YForm2_Source.Yij1;

	int n;
	int i,j;
	double YK,Gij,Bij,b_ij;
	
	for (n=1;n<=Nb;n++)
	{
		i=Branch[n].n1;
		j=Branch[n].n2;
		YK=Branch[n].YK;
		if (i<0)
		{
			i=abs(i);
			
			Gij=Yij[n].G;
			Bij=Yij[n].B;
			b_ij=Yij1[n].B;
			
			Yii[i].G+=(1.0-1.0/YK)*Gij;
			Yii[i].B+=(1.0-1.0/YK)*Bij;
			Yii1[i].B+=(1.0-1.0/YK)*b_ij;
			
			Yii[j].G+=(1.0-YK)*Gij;
			Yii[j].B+=(1.0-YK)*Bij;
			Yii1[j].B=Yii1[j].B+(1.0-YK)*b_ij;
		}
		else if(j<0)
		{
			j=abs(j);
			Gij=Yij[n].G;
			Bij=Yij[n].B;
			b_ij=Yij1[n].B;
			
			Yii[j].G+=(1.0-1.0/YK)*Gij;
			Yii[j].B+=(1.0-1.0/YK)*Bij;
			Yii1[j].B+=(1.0-1.0/YK)*b_ij;
			
			Yii[i].G+=(1.0-YK)*Gij;
			Yii[i].B+=(1.0-YK)*Bij;
			Yii1[i].B=Yii1[j].B+(1.0-YK)*b_ij;
		}
		
		else
		{
			Bij=YK/2.0;
			b_ij=YK/2.0;
			
			Yii[i].B+=Bij;
			Yii[j].B+=Bij;
			
			Yii1[i].B+=b_ij;
			Yii1[j].B+=b_ij;			
		}
	}
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////形成因子表//////////////////////////////////////////////////////////////////////////

F_YinZi_Result_Type PowerFlow_Type::F_YinZi(F_YinZi_Source_Type F_YinZi_Source)
{
	int flag=F_YinZi_Source.flag,Nn=F_YinZi_Source.Nn,Npv=F_YinZi_Source.Npv,Np=F_YinZi_Source.Np;
	PVNode_Type *PVnode=F_YinZi_Source.PVnode;
	BalanceNode_Type *BalanceNode=F_YinZi_Source.BalanceNode;
	int *NYseq=F_YinZi_Source.NYseq,*NYsum=F_YinZi_Source.NYsum;
	Yii_Type *Yii=F_YinZi_Source.Yii;
	Yij_Type *Yij=F_YinZi_Source.Yij;

	int *NUsum=new int[Nn+1];
	int *NUseq=new int[Nn+2];
	NUseq[1]=1;
	double *D=new double[Nn+1];
	U_Type *U=new U_Type[(Nn+1)*(Nn+1)/2];

	int i,j,k,n2;
	int i_pv,n_pv;
	int i_Balance,n_Balance;
	double *B=new double[Nn+1];
	double Btemp;


	
	i_pv=1;n_pv=PVnode[i_pv].n;
	i_Balance=1;n_Balance=BalanceNode[i_Balance].n;
	
	for (i=1;i<=Nn;i++)
	{
		if (i==n_Balance&&i_Balance<=Np)
		{
			i_Balance++;
			n_Balance=BalanceNode[i_Balance].n;
			D[i]=0.0;
			NUsum[i]=0;
			NUseq[i+1]=NUseq[i];
			continue;
		}
		else if (flag==2&&i==n_pv&&i_pv<=Npv)
		{
			i_pv++;
			n_pv=PVnode[i_pv].n;
			NUsum[i]=0;
			NUseq[i+1]=NUseq[i];
			D[i]=0.0;
			continue;
		}
		else
		{
			for (j=i+1;j<=Nn;j++)
			{
				B[j]=0.0;
			}
			B[i]=Yii[i].B;

			for (j=NYseq[i];j<NYseq[i+1];j++)
			{
				B[Yij[j].n2]=Yij[j].B;
			}
			for (j=1;j<=Np;j++)
			{
				B[BalanceNode[j].n]=0.0;
			}
			if (flag==2)
			{
				for (j=1;j<=Npv;j++)
				{
					B[PVnode[j].n]=0.0;
				}
			}



			for (j=1;j<i;j++)
			{
				int flag_Uji=0; 
				double Uji; 
				for (k=NUseq[j];k<NUseq[j+1];k++)
				{
					if (U[k].n2==i)
					{
						flag_Uji=1; 
						Uji=U[k].value; 
						Btemp=Uji/D[j];
						break;
					}
				}

				if (flag_Uji==1)
				{
					for (k=NUseq[j];k<NUseq[j+1];k++)
					{
						n2=U[k].n2;
						B[n2]-=Btemp*U[k].value;
					}
				}
			}

			for (j=i+1;j<=Nn;j++)
			{
				B[j]=B[j]/B[i];
			}
			B[i]=1.0/B[i];


			D[i]=B[i];
			int count=NUseq[i];
			for (j=i+1;j<=Nn;j++)
			{
				if (B[j]!=0)
				{
					U[count].value=B[j];
					U[count].n2=j;
					count++;
				}
			}
			NUseq[i+1]=count;
			NUsum[i]=NUseq[i+1]-NUseq[i];
		}
	}

	F_YinZi_Result_Type F_YinZi_Result;
	{
		F_YinZi_Result.D=D;
		F_YinZi_Result.U=U;
		F_YinZi_Result.NUseq=NUseq;
		F_YinZi_Result.NUsum=NUsum;
	}
	return F_YinZi_Result;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////傳電壓初值//////////////////////////////////////////////////////////////////////////
F_PassV0_Result_Type PowerFlow_Type::F_PassV0(F_PassV0_Source_Type F_PassV0_Source)
{
	int Nn=F_PassV0_Source.Nn;