#define L 5					//網絡的支路總數
#define N 5					//節點數
#define M 3					//網絡的PQ節點數
#define N0 2*(N-1)			//雅可比矩陣的行數
#define N1 N0+1
int _tmain(int argc, _TCHAR* argv[])
{
	int K=0;					//打印開關.K=1,則打印;否則,不打印
	int K1=0;					//子程序PLSC中判斷輸入電壓的形式;K1=1,則為極座標形式.否則為直角坐標形式. 
	float D=1e-4;				//有功及無功功率誤差的最大值（絕對值）.
	float *DD;					//誤差量最大者
	float p=1.0;;
	DD=&p;
	float G[f2(N,N,N)+1];				//Ybus的電導元素(實部).  
	float B[f2(N,N,N)+1];				//Ybus的電納元素(虛部).   
	float G1[L]={0.624024,0.754717,0.829876,0,0};	//第I支路的串聯電導.      
	float B1[L]={-3.900156,-2.641509,-3.112033,-63.492063,-31.746032};	//第I支路的串聯電納

	//float C1[L]={0.25,0,0.25,3.023432,1.511716};			//C1(I)第I支路的pie型對稱接地電納.?????????
	float C1[L]={0.25,0,0.25,0,0};			//C1(I)第I支路的pie型對稱接地電納

	//float C[f2(N,L,L)+1]={0,0,0,0,0, 0,0,0,3.023432,0, 0,0,0,0,0,
	//				0,0,0,-6.198035, 0,0,0,0,0,-3.099017};	//C(I,J)第I節點J支路不對稱接地電納.  ???????
	float C[f2(N,L,L)+1]={0,0,0,0,0, 0,0,0,3.023432,0, 0,0,0,0,1.511716,
					0,0,0,-3.174603,0, 0,0,0,0,-1.587302};	//C(I,J)第I節點J支路不對稱接地電納

	float CO[N]={0,0,0,0,0};				//CO(I) :第I節點的接地電納
	int S1[L]={1,1,2,2,3};					//第I支路的起始節點號
	int E1[L]={2,3,3,4,5};					//第I支路的終止節點號
	float P[N]={-1.6,-2,-3.7,5.0,0};				//第I節點的注入有功功率
	float Q[N]={-0.8,-1,-1.3,0,0};				//第I節點的注入無功功率
	float P0[N]={0};				//第I節點有功功率誤差
	
	float V0[N]={0};				//第I節點(PV節點)的電壓誤差(平方誤差)
	float V[N]={1,1,1,1.05,1.05};				//第I節點的電壓幅值
	//float E[N]={1,1,1,1.05,1.05};				//第I節點的電壓的實部
	float E[N]={1,1,1,1,1.05};				//第I節點的電壓的實部
	float F[N]={0,0,0,0,0};				//第I節點的電壓的虛部
	float JM[f2(N0,N0,N0)+1];			//Jacoby矩陣的第I行J列元素
	float A[f2(2*N-2,2*N-1,2*N-1)+1]={0};		//A(I,J)修正方程的增廣矩陣,三角化矩陣的第I行J列元素,
										//運算結束后A矩陣的最后一列存放修正的解.  
	float P1[L];				//第I支路由S1(I)節點注入的有功功率
	float Q1[L];				//第I支路由S1(I)節點注入的無功功率
	float P2[L];				//第I支路由E1(I)節點注入的有功功率
	float Q2[L];				//第I支路由E1(I)節點注入的無功功率
	float P3[L];				//第I支路的有功功率損耗
	float Q3[L];				//第I支路的無功功率損耗
	float ANGLE[N];			//第I節點電壓的角度
	float Q0[N]={0};			//第I節點有功功率誤差
	int ii,jj;
	int cnt=1;				//迭代次數
//	file4=fopen("cs_flow_1.txt","w");
	ybus(N,L,M,G,B,G1,B1,C1,C,CO,1,S1,E1);//計算電導和導納參數
	dpqc(P,Q,P0,Q0,V,V0,M,N,E,F,1,G,B,DD);
	//while (arrayMax(N,Q0)>D || arrayMax(N,P0)>D){
	while (*(DD)>D){
		cnt++;
		jmcc(M,N,N0,E,F,G,B,JM,1);//算jacob矩陣
		//算增廣的修正方程矩陣，最后一列放誤差向量，與書上的方程順序稍有不同
		for (ii=1;ii<=N0;ii++){
			for(jj=1;jj<=N0;jj++){
				A[f2(ii,jj,N1)]=JM[f2(ii,jj,N0)];
			}
		}
		for(ii=1;ii<=N0;ii++){
			if(ii<=2*M){
				if(ii%2==1)		A[f2(ii,2*N-1,2*N-1)]=Q0[f1((ii+1)/2)];
				else	A[f2(ii,2*N-1,2*N-1)]=P0[f1(ii/2)];
			}
			else{
				if(ii%2==1)		A[f2(ii,2*N-1,2*N-1)]=V0[f1((ii+1)/2)];
				else	A[f2(ii,2*N-1,2*N-1)]=P0[f1(ii/2)];
			}
		}
		sevc(A,N0,1,N1);//解線性方程組
		for(ii=1;ii<=N0;ii++){
			if(ii%2==1){
				E[f1((ii+1)/2)]-=A[f2(ii,2*N-1,2*N-1)];
			}
			else F[f1(ii/2)]-=A[f2(ii,2*N-1,2*N-1)];
		}//修正解
		dpqc(P,Q,P0,Q0,V,V0,M,N,E,F,1,G,B,DD);//重新計算誤差
	}

	/*float Atest[f2(4,5,5)+1]={1,2,0,0,1, 0,1,0,0,2, 0,0,1,0,1, 2,0,0,1,1};
	sevc ( Atest, 4, 1, 5);*/
	plsc(  N,  L,  M,  G ,  B ,  E ,  F ,\
             E1 ,  S1 ,  G1 ,  B1 ,  C1 ,  C ,\
             CO ,  P1 ,  Q1 ,  P2 ,  Q2 ,  P3 ,\
             Q3 ,  P ,  Q ,  V ,  ANGLE , K1);
	//file4=fclose("cs_flow_1");
	return 0;
}