/*讀取PQ節點的注入功率*/
 for(i=1;i<=npq;i++)
    {fscanf(fp,"%f",PE+i);
     fscanf(fp,"%f",QUE+i);
    }
 
 /*讀取PV節點的P && V并得出U^2*/
 for(i=npq+1;i<=n-1;i++)
    {fscanf(fp,"%f",PE+i);
     fscanf(fp,"%f",e+i);
     QUE[i]=e[i]*e[i];
    }
 
 /*讀取平衡節點的e && f*/
 fscanf(fp,"%f",e+n);
 fscanf(fp,"%f",f+n);
 

/*讀取各支路的起點號b，終點號e，支路阻抗R && X*/
 for(i=1;i<=z;i++)
    {fscanf(fp,"%d",&zl[i].b);
     fscanf(fp,"%d",&zl[i].e);
     fscanf(fp,"%f",&zl[i].RG);
     fscanf(fp,"%f",&zl[i].XB);
    }

 /*讀取計算精度E*/
 fscanf(fp,"%f",&E);
 
 /*讀取最大迭代次數km*/
 fscanf(fp,"%d",&km);

  /*給PQ節點的e賦初值*/
 for(i=1;i<=npq;i++)
    e[i]=1;
 
 printf("                          歡迎使用潮流精靈2007!\n"); 
 fprintf(fp1,"                          歡迎使用潮流精靈2007!\n"); 


 /*求取節點導納矩陣并輸出*/
 fprintf(fp1,"求節點導納矩陣：\n"); 
 dn(fp1);
 fprintf(fp1,"\n導納矩陣中的各元素:Y11=......Y12=......Ynn=......\n\n"); 
 
 fprintf(fp1,"形成導納矩陣YB:\nYB=\n");
 for(i=1;i<=n;i++)
   {for(j=1;j<=n;j++)
      fprintf(fp1,"   %f+j(%f)",G[i][j],B[i][j]);
    fprintf(fp1,"\n"); 
   }
 fprintf(fp1,"\n\n\n"); 

 
 
 
 /*該循環為迭代循環，直到達到規定精度或最大迭代次數為止*/
 do
{ /*輸出分隔線*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  

 fprintf(fp1,"第%d次迭代：K=%d\n\n",K,K);

 fprintf(fp1,"1、計算各PQ、PV節點功率的不平衡量，及PV節點電壓的不平衡量:\n\n");
 fprintf(fp1,"取：\n");
 for(i=1;i<=npq;i++)
    fprintf(fp1,"   U%d(%d)=e%d(%d)+f%d(%d)=%f+j%f\n",i,K,i,K,i,K,e[i],f[i]);
 for(i=npq+1;i<=n-1;i++)
    fprintf(fp1,"   U%d(%d)=e%d(%d)+f%d(%d)=%f+j%f\n",i,K,i,K,i,K,e[i],f[i]);
 fprintf(fp1,"   節點%d是平衡節點,保持U%d=e%d+f%d=%f+j%f為定值。\n\n",n,n,n,n,e[n],f[n]);
 
 
 
 /*計算PQ節點的P和Q*/
 fprintf(fp1,"a/計算各PQ、PV節點功率：\n");
 
 fprintf(fp1,"Pi(%d)=求和... Qi(%d)=求和... Ui(%d)=ei(%d)^2+fi(%d)^2\n\n",K,K,K,K,K);
 fprintf(fp1,"下面的符號別忘了帶角標!\n\n");
 
 PQ(fp1);
 /*計算PV節點的P和V*/
 PV(fp1);


 /*計算PQ節點的DP,DQ*/
 fprintf(fp1,"于是：DPi=Pi-Pi(%d);DQi=Qi-Qi(%d);DUi^2=Ui^2-(ei(%d)^2+fi(%d)^2)\n\n",K,K,K,K);
 for(i=1;i<=npq;i++)
    {DP[i]=PE[i]-P[i];
     DQU[i]=QUE[i]-QU[i];
     fprintf(fp1,"DP%d=P%d-P%d(%d)=%f-(%f)=%f\n",i,i,i,K,PE[i],P[i],DP[i]);
     fprintf(fp1,"DQ%d=Q%d-Q%d(%d)=%f-(%f)=%f\n",i,i,i,K,QUE[i],QU[i],DQU[i]);
    }
 /*計算PV節點的DP,D(U^2)*/
 for(i=npq+1;i<=n-1;i++)
    {DP[i]=PE[i]-P[i];
     DQU[i]=QUE[i]-QU[i];
     fprintf(fp1,"DP%d=P%d-P%d(%d)=%f-(%f)=%f\n",i,i,i,K,PE[i],P[i],DP[i]);
     fprintf(fp1,"D(U%d)^2=(U%d)^2-(e%d(%d)^2+f%d(%d)^2)=%f-((%f)^2+(%f)^2)=%f\n",i,i,i,K,i,K,QUE[i],e[i],f[i]);
    }
 fprintf(fp1,"\n\n");

 fprintf(fp1,"2、計算雅可比矩陣中各元素:\n\n");


 /*計算PQ節點的aii和bii*/
 for(i=1;i<=npq;i++)
   {ai[i]=(P[i]*e[i]+QU[i]*f[i])/(e[i]*e[i]+f[i]*f[i]);
    bi[i]=(P[i]*f[i]-QU[i]*e[i])/(e[i]*e[i]+f[i]*f[i]);  
    fprintf(fp1,"I%d=(P%d-jQ%d)/U%d*=(%f-j%f)/(%f-j%f)=%f+j%f=a%d%d+jb%d%d\n",i,i,i,i,P[i],QU[i],e[i],f[i],ai[i],bi[i],i,i,i,i); 
   }
 fprintf(fp1,"\n\n");

/*計算PV節點的aii和bii*/
 for(i=npq+1;i<=n-1;i++)
   {fprintf(fp1,"Q%d=",i);
    QU[i]=0;
    for(j=1;j<=n;j++)
       {QU[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
        fprintf(fp1,"+%f*(%f*%f-%f*%f)-%f*(%f*%f+%f*%f)\n   ",f[i],G[i][j],e[j],B[i][j],f[j],e[i],G[i][j],f[j],B[i][j],e[j]);
       }
     fprintf(fp1,"=%f\n",QU[i]);
    ai[i]=(P[i]*e[i]+QU[i]*f[i])/(e[i]*e[i]+f[i]*f[i]);
    bi[i]=(P[i]*f[i]-QU[i]*e[i])/(e[i]*e[i]+f[i]*f[i]);  
    fprintf(fp1,"I%d=(P%d-jQ%d)/U%d*=(%f-j%f)/(%f-j%f)=%f+j%f=a%d%d+jb%d%d\n",i,i,i,i,P[i],QU[i],e[i],f[i],ai[i],bi[i],i,i,i,i); 
   }
 fprintf(fp1,"\n\n");


fprintf(fp1,"雅可比矩陣的各個元素分別為：\n\n");
 
/*求PQ節點形成的雅可比矩陣*/
 PQykb(fp1);


/*求PV節點形成的雅可比矩陣*/
 PVykb(fp1);

fprintf(fp1,"3、所以可以得到K=%d時的雅可比矩陣:\n",K);
 /*輸出雅克比矩陣*/
fprintf(fp1,"J(%d)=\n",K);
for(i=1;i<=2*(n-1);i++)
   {for(j=1;j<=2*(n-1);j++)
      fprintf(fp1,"  %f",J[i][j]);
   fprintf(fp1,"\n");
   }
fprintf(fp1,"\n\n");

fprintf(fp1,"至此，可以建立修正方程組如下：....\n\n");


/*用列主元高斯消去法來接上面行成的線性方程組*/
for(i=1;i<=n-1;i++)
   {b1[2*i-1]=DP[i];
    b1[2*i]=DQU[i];
   }
lzy(2*(n-1),J,b1,fp1);
 
/*把解方程的結果賦給df與de并求最大不平衡量的絕對值dm*/
dm=0;
for(i=1;i<=n-1;i++)
   {df[i]=b1[2*i-1];
    de[i]=b1[2*i];
    if(fabs(df[i])>dm)
       dm=fabs(df[i]);
    if(fabs(de[i])>dm)
       dm=fabs(de[i]);
   }

fprintf(fp1,"解得：\n");
for(i=1;i<=n-1;i++)
   {fprintf(fp1,"   df%d=%f\n",i,df[i]);
    fprintf(fp1,"   de%d=%f\n",i,de[i]);
   }
fprintf(fp1,"\n\n");

fprintf(fp1,"因為：ei(%d)=ei(%d)+dei(%d); fi(%d)=fi(%d)+dfi(%d)\n所以：\n",K+1,K,K,K+1,K,K);
/*根據不平衡量求f與e的修正值并輸出*/
for(i=1;i<=n-1;i++)
   {f[i]+=df[i];
    e[i]+=de[i];
   }

for(i=1;i<=n-1;i++)
   {fprintf(fp1,"   e%d(%d)=e%d(%d)+de%d(%d)=%f+%f=%f\n",i,K+1,i,K,i,K,e[i]-de[i],de[i],e[i]);
    fprintf(fp1,"   f%d(%d)=f%d(%d)+df%d(%d)=%f+%f=%f\n",i,K+1,i,K,i,K,f[i]-df[i],df[i],f[i]);
   }
fprintf(fp1,"\n\n");

K++;

}
while(K<km && dm>E);

 
 /*輸出分隔線*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 
 /*輸出各節點的電壓*/
 fprintf(fp1,"各節點電壓為:\n");
 for(i=1;i<=n;i++)
    fprintf(fp1,"   U%d=e%d+jf%d=%f+j(%f)\n",i,i,i,e[i],f[i]);


/*輸出分隔線*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 /*計算平衡功率S~n*/
 fprintf(fp1,"平衡節點功率為:\n");
 i=n;
 PE[i]=0;
 for(j=1;j<=n;j++)
    PE[i]+=e[i]*(G[i][j]*e[j]-B[i][j]*f[j])+f[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
 QUE[i]=0;
 for(j=1;j<=n;j++)
    QUE[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
 fprintf(fp1,"S~%d=......\n",i);
 fprintf(fp1,"   =(%f+j%f)[",e[i],f[i]);
 for(j=1;j<=n;j++)
    fprintf(fp1,"(%f-j(%f))(%f-j(%f))",G[i][j],B[i][j],e[j],f[j]);
    
 fprintf(fp1,"]\n   =%f+j(%f)\n",PE[i],QUE[i]);




/*輸出分隔線*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n");  
 
 /*計算線路功率*/
 xlgl(fp1);

 
 /*求網絡總損耗*/
 fprintf(fp1,"\n法二:\n");
 for(i=npq+1;i<=n-1;i++)
    {fprintf(fp1,"Q%d=",i);
	 QUE[i]=0;
     for(j=1;j<=n;j++)
     {QUE[i]+=f[i]*(G[i][j]*e[j]-B[i][j]*f[j])-e[i]*(G[i][j]*f[j]+B[i][j]*e[j]);
      fprintf(fp1,"+%f*(%f*%f-%f*%f)-%f*(%f*%f+%f*%f)\n   ",f[i],G[i][j],e[j],B[i][j],f[j],e[i],G[i][j],f[j],B[i][j],e[j]);
     }
	 fprintf(fp1,"=%f\n\n\n",QUE[i]);
    }
 SP=0;SQ=0;
 fprintf(fp1,"DDS~=");
 for(i=1;i<=n;i++)
    {SP+=PE[i];
     SQ+=QUE[i];
	 fprintf(fp1,"+(%f+j%f)",PE[i],QUE[i]);
    }
 fprintf(fp1,"\n    =%f+%f",SP,SQ);

 fprintf(fp1,"\n\n******************************************THE END**********************************************");
   
 fclose(fp);
 fclose(fp1);
 
 printf("\n\n\n\n           潮流計算已完成！請打開輸出文件out.txt查看計算結果！\n\n\n\n\n\n                       ");
 
}