GenPower=malloc((N+1)*sizeof(struct GeneratorPower)); /*為發電機所發功率數組分配內存*/  
  fprintf(fp,"\t系統潮流計算結果：\n(1)迭代過程紀錄：\n迭代次數\t    \t有功迭代\t\t無功迭代\n\t\t   ΔPmax\tP-Node\t ΔQmax\t\tQ-Node\n");
  for(k=0;1;k++)
  {
    fprintf(fp,"   %2d:\t",k);
    if(Kp==1)
    {
        NodePower(1,N,NodeVol,NodePow,Yii,Yij,NYseq);/*節點功率計算函數*/
        Iteration(1,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI1,&MaxError,&ErrNode);/*迭代*/
        fprintf(fp,"\t%10.7lf\t %d\t",MaxError,ErrNode);
        if(MaxError>=epsilon)
            FormulaSolution(1,U1,D1,NUsum1,DI1,N,NodeVol,V0);/*線性方程求解*/
        else
            Kp=0;
    }
    else
        fprintf(fp,"\t\t\t\t");
    if(Kq==1)
    {
        NodePower(2,N,NodeVol,NodePow,Yii,Yij,NYseq);/*節點功率計算函數*/
        Iteration(2,Generator,Load,PVNode,NodeVol,NodePow,GenPower,N,DI2,&MaxError,&ErrNode);/*迭代*/
        fprintf(fp,"%10.7lf\t %d\n",MaxError,ErrNode);
        if(MaxError>=epsilon)
            FormulaSolution(2,U2,D2,NUsum2,DI2,N,NodeVol,V0);/*線性方程求解*/
        else
            Kq=0;

    }
    else
        fprintf(fp,"\n");
    if(Kp==0&&Kq==0)
        break;
    if(k>2000)
    {
        fprintf(fp,"\n迭代次數超過2000次，系統不收斂!\n");
        fprintf(fp,"本系統為%d個節點、%d條支路系統，含有%d個補償電容和%d個并聯電抗\n",N,Nb,Nc,Nr);
        now=time(NULL);
        fprintf(fp,"\t\t處理時間： %s\n",ctime(&now));
        /*fclose(fp);
        exit(1);*/
        goto DATA;
    }

    
    printf("COMPUTING......\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
  }


  fprintf(fp,"\n總計迭代次數: %d 次!\n\n",k+1);
  fprintf(fp,"(2)潮流計算結果(節點電壓)：\n Node\t\tV\t\tθ\t\tP\t\tQ\n");
  NodeDataOutput(fp,NodeVol,Generator,N,GenPower,NodePow,Load,Nl);
  fprintf(fp,"\n(3)潮流計算結果(支路功率):\n    Branch\t\tP\t\tQ\t\tQc\t\tQr\n");
  BranchDataOutput(fp,Nb,Nc,Nr,Branch,Compensation,Reactance,NodeVol);
  fprintf(fp,"\n本系統為%d個節點、%d條支路系統，含有%d個補償電容和%d個并聯電抗\n",N,Nb,Nc,Nr);
  now=time(NULL);
  fprintf(fp,"\n\t\t\t處理時間： %s\n",ctime(&now));

DATA:
  fprintf(fp,"\n\n附原始數據：\n");
  fprintf(fp,"節點數=%d 個,支路數=%d 條,補償電容=%d 個,并聯電抗=%d 個，發電機=%d 個,負荷=%d 個,系統平均電壓=%lf,誤差=%lf\n",N,Nb,Nc,Nr,Ng,Nl,V0,epsilon);
  for(n=1;n<=Nb;n++)
  {
    if((Branch[n].i>0)&&(Branch[n].j>0))
    {
        fprintf(fp,"節點%3d和節點%3d之間的電力線參數:",Branch[n].i,Branch[n].j);
        fprintf(fp,"R+jX=%e+j%e   Y0=%e\n",Branch[n].R,Branch[n].X,Branch[n].YK);
    }
    else
    {
        fprintf(fp,"節點%3d和節點%3d之間的變壓器參數:",abs(Branch[n].i),abs(Branch[n].j));
        fprintf(fp,"R+jX=%e+j%e   K=%lf,",Branch[n].R,Branch[n].X,fabs(Branch[n].YK));
    if(Branch[n].i<0)
        i=abs(Branch[n].i);
    else
        i=abs(Branch[n].j);
    fprintf(fp,"非標準變比位于節點%3d側\n",i);
    }
  }
  for(n=1;n<=Nc;n++)
  {
    fprintf(fp,"位于節點%3d的補償電容:",Compensation[n].i);
    fprintf(fp,"Y=j%e\n",Compensation[n].Y);
  }
  for(n=1;n<=Nr;n++)
  {
    fprintf(fp,"位于節點%3d的并聯電抗:",Reactance[n].i);
    fprintf(fp,"X=j%e\n",Reactance[n].X);
  }
  for(n=1;n<=Nl;n++)
  {
    fprintf(fp,"節點%3d的負荷  :",Load[n].i);
    fprintf(fp,"P+jQ=%lf-j%lf   V=%lf \n",Load[n].P,fabs(Load[n].Q),Load[n].V);
  }
  for(n=1;n<=Ng;n++)
  {
    fprintf(fp,"節點%3d的發電機:",Generator[n].i);
    fprintf(fp,"P+jQ=%lf+j%lf   V=%lf \n",Generator[n].P,Generator[n].Q,Generator[n].V);
  }
  for(n=1;n<Npv;n++)
  {
    if(n==1)
    fprintf(fp,"PV節點有%3d個,分別為：",Npv-1);
    fprintf(fp,"節點%3d,",PVNode[n].i); 
  }

  fprintf(fp,"\n\n\n");

  fclose(fp);
  free(Branch);
  free(Compensation);
  free(Reactance);
  free(Generator);
  free(Load);
  free(PVNode);
  free(Yii);
  free(Yiil);
  free(Yij);
  free(Yijl);
  free(NYseq);
  free(U1);
  free(U2);
  free(D1);
  free(D2);
  free(NUsum1);
  free(NUsum2);
  free(DI1);
  free(DI2);
  free(NodePow);
  free(NodeVol);
  free(GenPower);
}



void Datain(int Nb,int Nl,int Nc,int Nr,int Ng,FILE *fp,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance,struct Load_Type *Load,struct Generator_Type *Generator,struct PVNode_Type *PVNode,int *Npv)
{
  int n;
  /**********下面從數據文件中讀取支路、發電機、負荷的數據***********/

  for(n=1;n<=Nb;n++)
  {
    fscanf(fp,"%d,%d,%lf,%lf,%lf",&(Branch[n].i),&(Branch[n].j),&(Branch[n].R),&(Branch[n].X),&(Branch[n].YK));
  }
  for(n=1;n<=Nc;n++)
  {
    fscanf(fp,"%d,%lf",&(Compensation[n].i),&(Compensation[n].Y));
  }
  for(n=1;n<=Nr;n++)
  {
    fscanf(fp,"%d,%lf",&(Reactance[n].i),&(Reactance[n].X));
  }
  for(n=1;n<=Ng;n++)
  {
    fscanf(fp,"%d,%lf,%lf,%lf",&(Generator[n].i),&(Generator[n].P),&(Generator[n].Q),&(Generator[n].V));
    if((Generator[n].V)<0)
    {
        (*Npv)++;
        PVNode[(*Npv)].i=Generator[n].i;
        PVNode[(*Npv)].V=-(Generator[n].V);
    }
  }
  for(n=1;n<=Nl;n++)
  {
    fscanf(fp,"%d,%lf,%lf,%lf",&(Load[n].i),&(Load[n].P),&(Load[n].Q),&(Load[n].V));
    if((Load[n].V)<0)
    {
        (*Npv)++;
        PVNode[(*Npv)].i=Load[n].i;
        PVNode[(*Npv)].V=-(Load[n].V);
    }
  }



  fclose(fp);/*關閉數據輸入流*/
/******************讀取支路、發電機、負荷的數據完成，并關閉系統數據文件********************/



}







void AdmittanceMatrix(int N,int Nb,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,int *NYseq1,int *NYsum1)
{
  int i,n;


  int *NYsum=malloc((N+1)*sizeof(int));
  int *NYseq=malloc((N+1)*sizeof(int));
  for(i=1;i<=N;i++)
  {
    Yii[i].G=0.0;
    Yii[i].B=0.0;
    Yiil[i].G=0.0;
    Yiil[i].B=0.0;
    NYsum[i]=0;
  }

/************下面形成不考慮接地支路時的導納矩陣*******************/
  for(n=1;n<=Nb;n++)
  {
    int i=abs(Branch[n].i);
    int j=abs(Branch[n].j);
    double R=(Branch[n].R);
    double X=(Branch[n].X);
    double YK=(Branch[n].YK);

    double Zmag2=R*R+X*X;
    double Gij=R/Zmag2;
    double Bij=-X/Zmag2;

    double b_ij=-1.0/X;

    if((Branch[n].i<0)||(Branch[n].j<0))
    {
        Yij[n].G=-Gij/YK;
        Yij[n].B=-Bij/YK;
        Yijl[n].G=0.0;
        Yijl[n].B=-b_ij/YK;
    }
    else
    {
        Yij[n].G=-Gij;
        Yij[n].B=-Bij;
        Yijl[n].G=0.0;
        Yijl[n].B=-b_ij;

        /*這幾句在最后一次執行時，執行ps.dat時NYsum[1]、[2]變為0，NYsum[3]=-21845,NYsum[4]=-21845,但是[1]、[2]、[3]均應為2，[4]應為0;在執行ps-5.dat時則沒有這問題*//*現在通過在此子程序內定義新的NYsum和NYseq已經得到了解決*/

    }
    Yij[n].j=j;
    Yijl[n].j=j;

    if((Branch[n].i<0)||(Branch[n].j<0))
    {
        Yii[i].G=Yii[i].G+Gij/YK;
        Yii[i].B=Yii[i].B+Bij/YK;
        Yii[j].G=Yii[j].G+Gij/YK;
        Yii[j].B=Yii[j].B+Bij/YK;

        Yiil[i].B=Yiil[i].B+b_ij/YK;
        Yiil[j].B=Yiil[j].B+b_ij/YK;
    }
    else
    {
        Yii[i].G=Yii[i].G+Gij;
        Yii[i].B=Yii[i].B+Bij;
        Yii[j].G=Yii[j].G+Gij;
        Yii[j].B=Yii[j].B+Bij;

        Yiil[i].B=Yiil[i].B+b_ij;
        Yiil[j].B=Yiil[j].B+b_ij;
    }
    NYsum[i]+=1;

  }

  NYseq[1]=1;

/****************這里給各行的首元素定位********************/

  for(i=1;i<=N-1;i++)
  {
    NYseq[i+1]=NYseq[i]+NYsum[i];
  }

  for(n=1;n<=N;n++)
 {
  NYseq1[n]=NYseq[n];
 }




}


void AdmittanceMatrixAdd(int Nb,int Nc,int Nr,struct Yii_Type *Yii,struct Yii_Type *Yiil,struct Yij_Type *Yij,struct Yij_Type *Yijl,struct Branch_Type *Branch,struct Compensation_Type *Compensation,struct Reactance_Type *Reactance)
{
  int n;
  /*char ch;*/




/*******************下面在導納矩陣中追加接地支路*********************/

  for(n=1;n<=Nb;n++)
  {
    int i=Branch[n].i;
    int j=Branch[n].j;
    double YK=Branch[n].YK;

    if(!(i<0||j<0))
    {
        double Bij=YK/2.0;
        double b_ij=YK/2.0;

        Yii[i].B=Yii[i].B+Bij;
        Yii[j].B=Yii[j].B+Bij;

        Yiil[i].B=Yiil[i].B+b_ij;
        Yiil[j].B=Yiil[j].B+b_ij;
    }
    else
    {
        double Gij;
        double Bij;
        double b_ij;
        if(i<0)
        {
            i=abs(i);
            Gij=Yij[n].G;
            Bij=Yij[n].B;

            b_ij=Yijl[n].B;

            Yii[i].G=Yii[i].G+(1.0-1.0/YK)*Gij;
            Yii[i].B=Yii[i].B+(1.0-1.0/YK)*Bij;
            Yiil[i].B=Yiil[i].B+(1.0-1.0/YK)*b_ij;


            Yii[j].G=Yii[j].G+(1.0-YK)*Gij;
            Yii[j].B=Yii[j].B+(1.0-YK)*Bij;
            Yiil[j].B=Yiil[j].B+(1.0-YK)*b_ij;
        }
        else
        {
            j=abs(j);
            Gij=Yij[n].G;
            Bij=Yij[n].B;

            b_ij=Yijl[n].B;

            Yii[j].G=Yii[j].G+(1.0-1.0/YK)*Gij;
            Yii[j].B=Yii[j].B+(1.0-1.0/YK)*Bij;
            Yiil[j].B=Yiil[j].B+(1.0-1.0/YK)*b_ij;


            Yii[i].G=Yii[i].G+(1.0-YK)*Gij;
            Yii[i].B=Yii[i].B+(1.0-YK)*Bij;
            Yiil[i].B=Yiil[i].B+(1.0-YK)*b_ij;
        }
    }
  }
/********************追加接地支路完成**********************/
/********************追加對地補償電容***********************/
  for(n=1;n<=Nc;n++)
  {
    int i=Compensation[n].i;
    double Bij=Compensation[n].Y;
    double b_ij=Compensation[n].Y;

    Yii[i].B=Yii[i].B+Bij;
    
    Yiil[i].B=Yiil[i].B+b_ij;
   }
/********************追加對地并聯電感***********************/
  for(n=1;n<=Nr;n++)
  {
    int i=Reactance[n].i;
    double Bij=-1/Reactance[n].X;
    double b_ij=-1/Reactance[n].X;

    Yii[i].B=Yii[i].B+Bij;
    
    Yiil[i].B=Yiil[i].B+b_ij;
    
  }



}



void Factorial(int flag,int N,int Npv,struct PVNode_Type *PVNode,int *NUsum,struct Yii_Type *Yii,struct Yij_Type *Yij,int *NYseq,double *D,struct U_Type *U)/*因子表形成函數*/
{
  int n_pv,i_pv,j,n_u,i_above;
  int i,count;
  double *B=malloc((N+1)*sizeof(double));
  double Btemp;
  
  n_pv=1;
  i_pv=PVNode[1].i;
  
  for(i=1;i<N;i++)
  {
    if((flag==2)&&(i==i_pv))
    {
        n_pv++;
        i_pv=PVNode[n_pv].i;
        NUsum[i]=0;
        D[i]=0.0;
        continue;
    }
    else
    {
        for(count=i+1;count<N;count++)
            B[count]=0.0;
        
        B[i]=Yii[i].B;
        
        for(count=NYseq[i];count<NYseq[i+1];count++)
        {
            j=Yij[count].j;
            B[j]=Yij[count].B;
        }
        
        if(flag==2)
        {
            for(count=1;count<=Npv;count++)
            {
                j=PVNode[count].i;
                B[j]=0.0;
            }
            
        }

        
        n_u=1;
        i_above=1;
        
        while(i_above<=(i-1))
        {
            count=1;
            
            while(count<=NUsum[i_above])
            {
                if(U[n_u].j==i)
                {   
                    
                    Btemp=U[n_u].value/D[i_above];