/*照書P164的例題先試驗，已知G【i】【j】，B【i】【j】,PE[i],QE[i],初始e【i】,f[i],平衡節(jié)點的電壓e【1】，f【1】*/

#include<math.h>
#include<stdio.h>
#define M 50/*設(shè)定最大節(jié)點數(shù)——M-1*/
#define N 50/*設(shè)定最大支路數(shù)——N-1*/
#define M1 2*M/*設(shè)定雅克比矩陣的最大維數(shù)*/

float  G[M][M]={0},B[M][M]={0},PE[M],QUE[M],P[M],QU[M],
       DP[M],DQU[M],ai[M],bi[M],e[M],f[M]={0},de[M],df[M],
	   J[M1][M1],b1[M1],SP,SQ,DDP[N],DDQ[N];
int    n,z,npq,npv;

struct zhilu
 {int b;
  int e;
  float RG;
  float XB;
}  zl[N];


/*該子程序用于根據(jù)支路數(shù)組來求取節(jié)點導納矩陣*/
void dn(FILE *fp1)
{int i,b1,e1;
 float r,x,g,b;
 for(i=1;i<=z;i++)
    {b1=zl[i].b;
     e1=zl[i].e;
     r=zl[i].RG;
     x=zl[i].XB;
     fprintf(fp1,"z%d%d=%f+j(%f)  ==>  ",b1,e1,r,x);
     g=r/(r*r+x*x);
     b=-x/(r*r+x*x);
	 fprintf(fp1,"y%d%d=1/Z%d%d=1/(%f+j(%f))=%f+j(%f)\n",b1,e1,b1,e1,r,x,g,b);
     zl[i].RG=g;
	 zl[i].XB=b;
     G[b1][b1]+=g;
     G[e1][e1]+=g;
     G[b1][e1]+=-g;
     G[e1][b1]+=-g;
     B[b1][b1]+=b;
     B[e1][e1]+=b;
     B[e1][b1]+=-b;
     B[b1][e1]+=-b;
     }
}  


/*該子函數(shù)用于根據(jù)節(jié)點導納矩陣以及PQ點的電壓來求取各點的P和Q*/
void  PQ(FILE *fp1)
{int i,j;
 for(i=1;i<=npq;i++)
    {P[i]=0;
     QU[i]=0;
     fprintf(fp1,"P%d=",i);
     for(j=1;j<=n;j++)
        { P[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]);
          fprintf(fp1,"+%f*(%f*%f-%f*%f)+%f*(%f*%f+%f*%f)\n   ",e[i],G[i][j],e[j],B[i][j],f[j],f[i],G[i][j],f[j],B[i][j],e[j]);
        }
     fprintf(fp1,"=%f\n\n",P[i]);
     fprintf(fp1,"Q%d=",i);
     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\n\n",QU[i]);
     }
} 

/*該子函數(shù)用于根據(jù)節(jié)點導納矩陣以及PV點的電壓來求取各點的P和V^2*/
void  PV(FILE *fp1)
{int i,j;
 for(i=npq+1;i<=n-1;i++)
    {P[i]=0;
     fprintf(fp1,"P%d=",i);
     for(j=1;j<=n;j++)
	     {P[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]);
	      fprintf(fp1,"+%f*(%f*%f-%f*%f)+%f*(%f*%f+%f*%f)\n   ",e[i],G[i][j],e[j],B[i][j],f[j],f[i],G[i][j],f[j],B[i][j],e[j]);
	     }
     fprintf(fp1,"=%f\n\n",P[i]);
	 fprintf(fp1,"U%d^2=",i);
     QU[i]=e[i]*e[i]+f[i]*f[i];
     fprintf(fp1,"(%f)^2+(%f)^2=%f\n\n\n",e[i],f[i],QU[i]);
    }
} 

/*該子函數(shù)用于求PQ節(jié)點形成雅可比矩陣*/
void  PQykb(FILE *fp1)
{int i,j;
 for(i=1;i<=npq;i++)
   for(j=1;j<=n-1;j++) 
      {J[2*i-1][2*j-1]=-B[i][j]*e[i]+G[i][j]*f[i];
         J[2*i-1][2*j]=G[i][j]*e[i]+B[i][j]*f[i];
         J[2*i][2*j-1]=-J[2*i-1][2*j];
           J[2*i][2*j]=J[2*i-1][2*j-1];
        if(i==j)
           {J[2*i-1][2*j-1]+=bi[i];
              J[2*i-1][2*j]+=ai[i]; 
              J[2*i][2*j-1]+=ai[i];
                J[2*i][2*j]+=-bi[i];
             fprintf(fp1,"H%d%d=偏P%d/偏f%d=-B%d%d*e%d+G%d%d*f%d+b%d%d=-(%f)*%f+%f*%f+%f=%f\n",i,j,i,j,i,j,i,i,j,i,i,i,B[i][j],e[i],G[i][j],f[i],bi[i],J[2*i-1][2*j-1]);
             fprintf(fp1,"N%d%d=偏P%d/偏e%d=G%d%d*e%d+B%d%d*f%d+a%d%d=%f*%f+%f*%f+%f=%f\n",i,j,i,j,i,j,i,i,j,i,i,i,G[i][j],e[i],B[i][j],f[i],ai[i],J[2*i-1][2*j]);
             fprintf(fp1,"J%d%d=偏Q%d/偏f%d=-G%d%d*e%d-B%d%d*f%d+a%d%d=-(%f)*%f-(%f)*%f+%f=%f\n",i,j,i,j,i,j,i,i,j,i,i,i,G[i][j],e[i],B[i][j],f[i],ai[i],J[2*i][2*j-1]);
             fprintf(fp1,"L%d%d=偏Q%d/偏e%d=-B%d%d*e%d+G%d%d*f%d-b%d%d=-(%f)*%f+%f*%f-(%f)=%f\n\n",i,j,i,j,i,j,i,i,j,i,i,i,B[i][j],e[i],G[i][j],f[i],bi[i],J[2*i][2*j]);
           }  
         else 
             {fprintf(fp1,"H%d%d=偏P%d/偏f%d=-B%d%d*e%d+G%d%d*f%d=-(%f)*%f+%f*%f=%f\n",i,j,i,j,i,j,i,i,j,i,B[i][j],e[i],G[i][j],f[i],J[2*i-1][2*j-1]);
              fprintf(fp1,"N%d%d=偏P%d/偏e%d=G%d%d*e%d+B%d%d*f%d=%f*%f+%f*%f=%f\n",i,j,i,j,i,j,i,i,j,i,G[i][j],e[i],B[i][j],f[i],J[2*i-1][2*j]);
              fprintf(fp1,"J%d%d=偏Q%d/偏f%d=-G%d%d*e%d-B%d%d*f%d=-(%f)*%f-(%f)*%f=%f\n",i,j,i,j,i,j,i,i,j,i,G[i][j],e[i],B[i][j],f[i],J[2*i][2*j-1]);
              fprintf(fp1,"L%d%d=偏Q%d/偏e%d=-B%d%d*e%d+G%d%d*f%d=-(%f)*%f+%f*%f=%f\n\n",i,j,i,j,i,j,i,i,j,i,B[i][j],e[i],G[i][j],f[i],J[2*i][2*j]);
             }
      }

}

/*該子函數(shù)用于求PV節(jié)點形成雅可比矩陣*/
void  PVykb(FILE *fp1)
{int i,j;
 for(i=npq+1;i<=n-1;i++)
   for(j=1;j<=n-1;j++) 
      {J[2*i-1][2*j-1]=-B[i][j]*e[i]+G[i][j]*f[i];
         J[2*i-1][2*j]= G[i][j]*e[i]+B[i][j]*f[i];
         J[2*i][2*j-1]=0;
           J[2*i][2*j]=0;
        if(i==j)
           {J[2*i-1][2*j-1]+=bi[i];
              J[2*i-1][2*j]+=ai[i]; 
              J[2*i][2*j-1]+=2*f[i];
                J[2*i][2*j]+=2*e[i];
            fprintf(fp1,"H%d%d=偏P%d/偏f%d=-B%d%d*e%d+G%d%d*f%d+b%d%d=-(%f)*%f+%f*%f+%f=%f\n",i,j,i,j,i,j,i,i,j,i,i,i,B[i][j],e[i],G[i][j],f[i],bi[i],J[2*i-1][2*j-1]);
            fprintf(fp1,"N%d%d=偏P%d/偏e%d=G%d%d*e%d+B%d%d*f%d+a%d%d=%f*%f+%f*%f+%f=%f\n",i,j,i,j,i,j,i,i,j,i,i,i,G[i][j],e[i],B[i][j],f[i],ai[i],J[2*i-1][2*j]);
            fprintf(fp1,"R%d%d=偏U%d^2/偏f%d=*f%d=2*%f=%f\n",i,j,i,j,i,f[i],J[2*i][2*j-1]);
            fprintf(fp1,"S%d%d=偏U%d^2/偏e%d=*e%d=2*%f=%f\n\n",i,j,i,j,i,e[i],J[2*i][2*j]);
           }  
		else
		   {fprintf(fp1,"H%d%d=偏P%d/偏f%d=B%d%d*e%d+G%d%d*f%d=-(%f)*%f+%f*%f=%f\n",i,j,i,j,i,j,i,i,j,i,B[i][j],e[i],G[i][j],f[i],J[2*i-1][2*j-1]);
            fprintf(fp1,"N%d%d=偏P%d/偏e%d=%d%d*e%d+B%d%d*f%d=%f*%f+%f*%f=%f\n",i,j,i,j,i,j,i,i,j,i,G[i][j],e[i],B[i][j],f[i],J[2*i-1][2*j]);
            fprintf(fp1,"R%d%d=偏U%d^2/偏f%d=0\n",i,j,i,j);
            fprintf(fp1,"S%d%d=偏U%d^2/偏e%d=0\n\n",i,j,i,j);
		   }
       }
}


/*該子程序為列主元高斯消去法解線性方程組*/
void lzy(int n,float a[M1][M1],float b[M1],FILE *fp1)
{float d,t,s;
 int i,j,k,l;
 
 for(k=1;k<=n-1;k++)
     {d=a[k][k];l=k;
      for(i=k+1;i<=n;i++)
         if(fabs(a[i][k])>fabs(d))
           {d=a[i][k];
            l=i;
           }
      if(d==0)
         {printf("\nJI YI!!!");
          break;
         }
      if(l!=k)
        for(j=k;j<=n;j++) 
           {t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;
            t=b[l];b[l]=b[k];b[k]=t;
           }

      for(i=k+1;i<=n;i++)
         a[i][k]=a[i][k]/a[k][k];
      for(i=k+1;i<=n;i++)
         for(j=k+1;j<=n;j++)
            a[i][j]=a[i][j]-a[i][k]*a[k][j];
      for(i=k+1;i<=n;i++)  
         b[i]=b[i]-a[i][k]*b[k];
     
      }


   b[n]=b[n]/a[n][n];
   for(i=n-1;i>=1;i--)
       {s=0;
         for(j=i+1;j<=n;j++)
            s+=a[i][j]*b[j];
        b[i]=(b[i]-s)/a[i][i];
       }
   
}

/*該子函數(shù)用于計算線路功率*/
void xlgl(FILE *fp1)
{int i,j,b1,e1;
 float g0,b0,g,b,P1,Q1,P2,Q2,SP=0,SQ=0;
 fprintf(fp1,"線路上功率及損耗為：\n");
 for(i=1;i<=z;i++)
 {g0=0;
  b0=0;
  b1=zl[i].b;
  e1=zl[i].e;
  g=zl[i].RG;
  b=zl[i].XB;
  for(j=1;j<=z;j++)
     if (zl[j].b==0 && zl[j].e==b1)
        {g0=zl[j].RG;
         b0=zl[j].XB;
        }
  P1=e[b1]*(e[b1]*g0-f[b1]*b0+(e[b1]-e[e1])*g-(f[b1]-f[e1])*b)+f[b1]*(e[b1]*b0+f[b1]*g0+(e[b1]-e[e1])*b+(f[b1]-f[e1])*g);
  Q1=f[b1]*(e[b1]*g0-f[b1]*b0+(e[b1]-e[e1])*g-(f[b1]-f[e1])*b)-e[b1]*(e[b1]*b0+f[b1]*g0+(e[b1]-e[e1])*b+(f[b1]-f[e1])*g);;
  fprintf(fp1,"s~%d%d=U%d[U*%dy*%d0+(U*%d-U*%d)y*%d%d]\n",b1,e1,b1,b1,b1,b1,e1,b1,e1);
  fprintf(fp1,"    =(%f+j(%f)){(%f-j(%f))(%f-j(%f))+[(%f-j(%f))-(%f-j(%f))](%f-j(%f))}\n",e[b1],f[b1],e[b1],f[b1],g0,b0,e[b1],f[b1],e[e1],f[b1],g,b);
  fprintf(fp1,"    =%f+j(%f)\n",P1,Q1);
  g0=0;
  b0=0; 
  for(j=1;j<=z;j++)
    if (zl[j].b==0 && zl[j].e==e1)
      {g0=zl[j].RG;
       b0=zl[j].XB;
      }
  P2=e[e1]*(e[e1]*g0-f[e1]*b0+(e[e1]-e[b1])*g-(f[e1]-f[b1])*b)+f[e1]*(e[e1]*b0+f[e1]*g0+(e[e1]-e[b1])*b+(f[e1]-f[b1])*g);
  Q2=f[e1]*(e[e1]*g0-f[e1]*b0+(e[e1]-e[b1])*g-(f[e1]-f[b1])*b)-e[e1]*(e[e1]*b0+f[e1]*g0+(e[e1]-e[b1])*b+(f[e1]-f[b1])*g);;
  fprintf(fp1,"s~%d%d=U%d[U*%dy*%d0+(U*%d-U*%d)y*%d%d]\n",e1,b1,e1,e1,e1,e1,b1,e1,b1);
  fprintf(fp1,"    =(%f+j(%f)){(%f-j(%f))(%f-j(%f))+[(%f-j(%f))-(%f-j(%f))](%f-j(%f))}\n",e[e1],f[e1],e[e1],f[e1],g0,b0,e[e1],f[e1],e[b1],f[e1],g,b);
  fprintf(fp1,"    =%f+j(%f)\n",P2,Q2);
  DDP[i]=P1+P2;
  DDQ[i]=Q1+Q2;
  fprintf(fp1,"dS%d%d=S%d%d+S%d%d=(%f+j%f)+(%f+j%f)=%f+j%f\n\n",b1,e1,b1,e1,e1,b1,P1,Q1,P2,Q2,DDP[i],DDQ[i]);
 }
 for(i=1;i<=z;i++)
 {SP+=DDP[i];
  SQ+=DDQ[i];
 }
 
 /*輸出分隔線*/
 for(i=0;i<18;i++)
    fprintf(fp1,"------");
 fprintf(fp1,"\n"); 
 
 fprintf(fp1,"有兩種方法可以求出網(wǎng)絡(luò)總損耗，任選其一:\n\n"); 
 fprintf(fp1,"法一：DDS~=求和dS＝%f+j%f\n",SP,SQ);
     
}





main()
{float dm,E,g0,b0;
 int   i,j,K=0,km;
 FILE  *fp,*fp1;
 
 if ((fp=fopen("input.txt","r"))==NULL)
   {printf("\nCannot open this file!");
    exit(0);
   }

 if ((fp1=fopen("output.txt","w"))==NULL)
   {printf("\nCannot open this file!");
    exit(0);
   }
 
 /*讀取節(jié)點數(shù)n，支路數(shù)z，PQ和PV節(jié)點數(shù)*/
 fscanf(fp,"%d",&n);
 fscanf(fp,"%d",&z);
 fscanf(fp,"%d",&npq);
 fscanf(fp,"%d",&npv);
 if (n!=(npq+npv+1))
    printf("                            節(jié)點數(shù)目輸入有誤!請仔細檢查!");