/*  PQ 潮流程序 */
  /********************************************************************/               
  /* 信息輸入                                                         */
  /*   線路數，變壓器數，發電機數，負荷數，負荷靜特性標志，最大節點號 */
  /*   最大迭代次數，節點電壓上限，節點電壓下限，最大誤差             */
  /*                                                                  */ 
  /*     負荷靜特性標志：如不考慮負荷的靜特性添0，負荷靜特性可以不輸入*/
  /*     考慮時，有幾種類型添幾，必須輸入負荷靜特性的系數             */
  /*  支路數據                                                        */
  /*     I      J      R       X      B/2(添正)                       */
  /*     .                                                            */  
  /*     .                                                            */
  /*     .                                                            */
  /*     接地支路   I，J一致                                          */
  /*  變壓器數據
        I      J      RT      XT      KT
        .
        .
        .
         J為KT所在側
   	發電機數據
        I      標志      PG       QG     V
        .        
        .
        .
      平衡機次序任意，PG處添平衡機電壓相角，平衡機的QG 任意
      標志：平衡機  0    
            PQ節點   1
            PV節點  -1
	負荷數據
        I      特性      PL       QL       V
        .        1
        .
      負荷數  .
	有10種考慮負荷靜特性的方案,添1為不考慮負荷靜特性，一般添1
          負荷PL、QL均添正
	負荷靜特性
             PA      PB      PC       QA      QB      QC
         1.  0       0       1        0       0       1
             .
             .
            (程序中默認第一行為恒定功率)
          原始數據存放在pqcin.dat文件中
          輸出數據存放在pq                                   */

# include <stdio.h>
# include <string.h>
# include <math.h>
# include <time.h>
# include <process.h>
# include <stdlib.h>
# include <malloc.h>
# define sita 57.29577951
# define max(a,b)   ((a)>=(b)? a:b)
# define min(a,b)   ((a)<=(b)? a:b)


double **TwoArrayAlloc_D(int,int);
void TwoArrayFree_D(double **);
int **TwoArrayAlloc_I(int,int);
void TwoArrayFree_I(int **);
void data_read(int,int,int,int,int,struct data *,struct data *,struct data *,
     struct data *,double **);
void order(int,int,int,int,int,int *,int *,struct data *,struct data *);
void matrix_adm(int,int,int,struct data *,struct data *,int *,int *,double *,
     double *,double *,double *);
void fdlf(int,int,int,int,int,struct data *,struct data *,struct data *,
     struct data *,double **,int *,int *,double *,double *,double *,double *,
     int,double,int *,double *,double *,double **,double **,int *);
void fact(int,int,int,int,int *,int *,double **,struct data *,struct data *,
     struct data *,int *,int *,double *,double *,int *,int *,double *,
     double *,int,int *);
void cbr(int,int,struct data *,struct data *,int *,double *,double *);
void power_node(int,int,double **,double *,double *,int *,int *,double *,
     double *,double *,double *);
void mismach(int,int,int,int,double **,double **,double **,struct data *,
     struct data *,double *,double *,double *,int *);
void czm(int,int,int *,int *,int *);
void zero(double *,int,int);
void revise(int,double *,double *,int *,int *,double *,double *);
void eliminate(int,double *,double *,int,int *,int *,int *,double *);

struct data
{
    int i;
    int j;
    double a;
    double b;
    double c;
};

double **TwoArrayAlloc_D(int m,int n)
{
   double *x,**y;
   int i;
   x=(double *)calloc(m*n,sizeof(double));
   if(x==NULL)  {  printf("\nno memory in xd");   exit(0);  }
   y=(double **)calloc(m,sizeof(double *));
   if(y==NULL)  {  printf("\nno memory in yd");   exit(0);  }
   for(i=0;i<m;i++)
     y[i]=&x[n*i];
   return(y);
}

void TwoArrayFree_D(double **x)
{
   free(x[0]);
   free(x);
}

int **TwoArrayAlloc_I(int m,int n)
{
   int i,*x,**y;
   x=(int *)calloc(m*n,sizeof(int));
   if(x==NULL)  {  printf("\nno memory in xi");   exit(0);  }
   y=(int **)calloc(m,sizeof(int *));
   if(y==NULL)  {  printf("\nno memory in yi");   exit(0);  }
   for(i=0;i<m;i++)
     y[i]=&x[n*i];
   return(y);
}

void TwoArrayFree_I(int **x)
{
   free(x[0]);
   free(x);
}

FILE *fin,*fou,*fst;

void main()
{

   clock_t st,en;
   char pinf[20],pouf[20],ptef[20];
   float a1,a2,a3,a4,a5,a6;
   int num_line,num_tran,num_gene,num_load,num_loadp,no_max,iter_max;
   int schem,num_node,num_bran,conv,i,j;
   int *otn,*nto,*iy,*jy;
   double vmax,vmin,error_max,c1,c2;
   double **property,*dyr,*dyi,*uyr,*uyi,*vm,*va,**pq_node,**pq_load;
   struct data *line,*tran,*gene,*load,*p,*end;
//   printf("the top of the heap is: %ld bytes\n",(unsigned long) coreleft());
  
   //fin=fopen("pqcin.dat","r");
   fin=fopen("Edit2.txt","r");
   if(fin==NULL)
   { printf("cannot open file\n");
      exit(0);
   }  
   fou=fopen("pqcout.dat","w");
   if(fou==NULL)
   { printf("cannot open outfile\n");
      exit(0);
   }
//  fin=fopen(pinf,"r");   fou=fopen(pouf,"w");   
//   fst=fopen(ptef,"w");
//   從文件中讀原始數據   
   fscanf(fin,"%d %d %d %d %d %d %d %lf %lf %lf",&num_line,&num_tran,
	  &num_gene,&num_load,&num_loadp,&no_max,&iter_max,
	  &vmax,&vmin,&error_max);
//   printf("%d,%d",num_line,num_tran);
//   exit(0);
//    從文件中讀線路、變壓器、發電機、負荷原始數據    
   line=(struct data *)calloc(num_line+1,sizeof(struct data));
   if(line==NULL)  
   {  
	   printf("\nno memory in line");
	   exit(0);  
   }
   tran=(struct data *)calloc(num_tran+1,sizeof(struct data));
   if(tran==NULL)  
   {  
	   printf("\nno memory in tran");  
	   exit(0);  
   }
   gene=(struct data *)calloc(num_gene+1,sizeof(struct data));
   if(gene==NULL)  
   {  
	   printf("\nno memory in gene");   
	   exit(0);  
   }
   load=(struct data *)calloc(num_load+1,sizeof(struct data));
   if(load==NULL)  
   {  
	   printf("\nno memory in load");  
	   exit(0);  
   }
   property=TwoArrayAlloc_D(num_loadp+1,6);
   
   data_read(num_line,num_tran,num_gene,num_load,num_loadp,line,tran,gene,
	     load,property);
   en=clock();   
   a1=(en-st)/CLK_TCK;
   
// 節點優化  
   otn=(int *)calloc(no_max+1,sizeof(int));
   if(otn==NULL) 
   {  
	   printf("\nno memory in otn");  
	   exit(0);  
   }
   for(end=line+num_line,p=line+1;p<=end;p++)  
   { 
	   otn[p->i]=1;  
	   otn[p->j]=1; 
   }
   for(end=tran+num_tran,p=tran+1;p<=end;p++)  
   { 
	   otn[p->i]=1;  
	   otn[p->j]=1; 
   }
   for(num_node=0,i=1;i<=no_max;i++)  
	   if(otn[i]==1)  otn[i]=++num_node;
    nto=(int *)calloc(num_node+1,sizeof(int));
   if(nto==NULL)  
   {  
	   printf("\nno memory in nto");   
	   exit(0); 
   }
   st=clock();
   order(schem,num_node,num_line,num_tran,no_max,otn,nto,line,tran);
   en=clock();   
   a2=(en-st)/CLK_TCK;
   for(end=gene+num_gene,p=gene+1;p<=end;p++)  
	   p->i=otn[p->i];
   for(end=load+num_load,p=load+1;p<=end;p++)  
	   p->i=otn[p->i];

//  形成導納矩陣   
   num_bran=num_line+num_tran;
   iy=(int *)calloc(num_node+2,sizeof(int));
   if(iy==NULL)  {  printf("\nno memory in iy");   exit(0);  }
   jy=(int *)calloc(num_bran+1,sizeof(int));
   if(jy==NULL)  {  printf("\nno memory in jy");   exit(0);  }
   dyr=(double *)calloc(num_node+1,sizeof(double));
   if(dyr==NULL)  {  printf("\nno memory in dyr");   exit(0);  }
   dyi=(double *)calloc(num_node+1,sizeof(double));
   if(dyi==NULL)  {  printf("\nno memory in dyi");   exit(0);  }
   uyr=(double *)calloc(num_bran+1,sizeof(double));
   if(uyr==NULL)  {  printf("\nno memory in uyr");   exit(0);  }
   uyi=(double *)calloc(num_bran+1,sizeof(double));
   if(uyi==NULL)  {  printf("\nno memory in uyi");   exit(0);  }
   st=clock();
   matrix_adm(num_node,num_line,num_tran,line,tran,iy,jy,dyr,dyi,uyr,uyi);
   en=clock();   a3=(en-st)/CLK_TCK;
   
//   開始計算   
   vm=(double *)calloc(num_node+1,sizeof(double));
   if(vm==NULL)  {  printf("\nno memory in vm");   exit(0);  }
   va=(double *)calloc(num_node+1,sizeof(double));
   if(va==NULL)  {  printf("\nno memory in va");   exit(0);  }
   pq_node=TwoArrayAlloc_D(num_node+1,2);
   pq_load=TwoArrayAlloc_D(num_node+1,2);
//   printf("the top of the heap is: %ld bytes\n",(unsigned long) coreleft());
  
   {
      st=clock();
      fdlf(num_node,num_line,num_tran,num_gene,num_load,line,tran,gene,load,
	       property,iy,jy,dyr,dyi,uyr,uyi,iter_max,error_max,nto,vm,va,
	       pq_node,pq_load,&conv);
      en=clock();    
	  a4=(en-st)/CLK_TCK;
   }
 //  輸出節點功率  
   if(conv==1)
   {
      fprintf(fou,"\n                      ********* load flow *********\n\n");
      fprintf(fou,"node  voltage   phase(deg.)    gen.p       gen.q       load.p      load.q\n");
      for(j=1;j<=no_max;j++)
      {
	 i=otn[j];  if(i==0)  continue;
	 c1=pq_load[i][0];  c2=pq_load[i][1];
	 fprintf(fou,"%4d%9.4f%12.4f%12.4f%12.4f%12.4f%12.4f\n",j,vm[i],
		 va[i]*sita,pq_node[i][0]+c1,pq_node[i][1]+c2,c1,c2);
      }
   }
   cbr(num_line,num_tran,line,tran,nto,vm,va);
   en=clock();    a6=(en-st)/CLK_TCK;
//   計算和輸出支路功率   
  /* fprintf(fou,"\n\n            |     read      data      | %f",a1);
   fprintf(fou,"\n            |    node    ordering     | %f",a2);
   fprintf(fou,"\nthe time of |forming admittance matrix| %f",a3);
   fprintf(fou,"\n            |          fdlp           | %f",a4);
   fprintf(fou,"\n            |  calculation and print  | %f",a6);
  // for(i=1;i<=num_node;i++)  fprintf(fst,"%f   %f  \n",vm[i],va[i]);*/
}

//*    讀原始數據    
void data_read(int num_line,int num_tran,int num_gene,int num_load,
     int num_loadp,struct data *line,struct data *tran,struct data *gene,
     struct data *load,double **property)
{
   int i,j;
   for(i=1;i<=num_line;i++)  fscanf(fin,"%d %d %lf %lf %lf",
       &line[i].i,&line[i].j,&line[i].a,&line[i].b,&line[i].c);
   for(i=1;i<=num_tran;i++)  fscanf(fin,"%d %d %lf %lf %lf",
       &tran[i].i,&tran[i].j,&tran[i].a,&tran[i].b,&tran[i].c);
   for(i=1;i<=num_gene;i++)  fscanf(fin,"%d %d %lf %lf %lf",
       &gene[i].i,&gene[i].j,&gene[i].a,&gene[i].b,&gene[i].c);
   for(i=1;i<=num_load;i++)  fscanf(fin,"%d %d %lf %lf %lf",
       &load[i].i,&load[i].j,&load[i].a,&load[i].b,&load[i].c);
   for(i=1;i<=num_loadp;i++)
	   for(j=0;j<6;j++)  fscanf(fin,"%lf",&property[i][j]);
}

//       節點優化      
void order(int schem,int num_node,int num_line,int num_tran,int no_max,
	   int *otn,int *nto, struct data *line,struct data *tran)
{
   int num_adj,i,j,k,ii,jj,kk,id,jd,ip,sel,mid;
   int **adjacent,*degree;
   struct data *p,*end;

   num_adj=20;
   adjacent=TwoArrayAlloc_I(num_node+1,num_adj);
   degree=(int *)calloc(num_node+1,sizeof(int));
   if(degree==NULL) {  printf("\nno memory in degree");   exit(0);  }
/*    forming adjacent list and achieving minimum degree of each node    */
   for(end=line+num_line,p=line+1;p<=end;p++)
   {
      i=p->i;  
	  j=p->j;
      if(i==j) continue;
      i=otn[i];  
	  j=otn[j];
      id=degree[i];  
	  jd=degree[j];
      adjacent[i][id]=j;  
	  adjacent[j][jd]=i;
      degree[i]=id+1;  
	  degree[j]=jd+1;
   }
   for(end=tran+num_tran,p=tran+1;p<=end;p++)
   {
      i=otn[p->i];  
	  j=otn[p->j];
      id=degree[i];  
	  jd=degree[j];
      adjacent[i][id]=j;  
	  adjacent[j][jd]=i;
      degree[i]=id+1;  
	  degree[j]=jd+1;
   }
/*    starting order    */
   ip=1;
   for(;;)
   {
      mid=num_node;
      for(i=1;i<=num_node;i++)
      {
         id=degree[i];
         if((nto[i]==0)&&(id<mid)) 
		 {  
			 mid=id;  
			 sel=i;
		 }
      }
      nto[sel]=ip;
      if(schem==2)
      {
	   for(i=0;i<mid;i++)
       {
	       ii=adjacent[sel][i];  
		   id=degree[ii];
            for(j=0;j<id;j++)
            {
	         if(adjacent[ii][j]==sel) adjacent[ii][j]=adjacent[ii][id-1];
            }
            degree[ii]=id-1;
	   }
         for(i=0;i<mid-1;i++)
         {
            ii=adjacent[sel][i];  
			id=degree[ii];
            for(j=i+1;j<mid;j++)
            {
               jj=adjacent[sel][j];  
			   jd=degree[jj]; 
			   kk=0;
	       for(k=0;k<jd;k++)
               {
		  if(adjacent[jj][k]==ii)  {kk=1;  break;}
               }
	       if(kk==0)
               {
		         adjacent[ii][id]=jj;  
				 adjacent[jj][jd]=ii;
                 degree[ii]=id+1;  
				 degree[jj]=jd+1;
               }
			}
         }
      }
      if(ip==num_node) break;
      else {  ip++;  continue;}
   }
   for(i=1;i<=no_max;i++)
   {
      j=otn[i];
      if(j!=0) otn[i]=nto[j];
   }
   for(i=1;i<=no_max;i++)
   {
      j=otn[i];
      if(j!=0) nto[j]=i;
   }
/*    old number i, new number otn[i]    */
/*    new number i, old number nto[i]    */
   for(end=line+num_line,p=line+1;p<=end;p++)
   { 
	   p->i=otn[p->i];  
	   p->j=otn[p->j]; 
   }
   for(end=tran+num_tran,p=tran+1;p<=end;p++)
   {  
	   p->i=otn[p->i]; 
	   p->j=otn[p->j]; 
   }
   TwoArrayFree_I(adjacent);  free(degree);
}

//         形成導納矩陣         
void matrix_adm(int num_node,int num_line,int num_tran,struct data *line,
     struct data *tran,int *iy,int *jy,double *dyr,double *dyi,double *uyr,
     double *uyi)
{
   int ii,jj,count,i,j;
   double r,x,b;
   struct data *p,*end;
   count=1;
   for(i=1;i<=num_node;i++)
   {
      iy[i]=count;
      end=line+num_line;
      for(p=line+1;p<=end;p++)
      {
         ii=p->i; 
		 jj=p->j;
	    if(min(ii,jj)!=i) continue;
         r=p->a; 
		 x=p->b; 
		 b=r*r+x*x;
         r/=b;  
		 x/=-b;
         if(ii==jj)   
		 {  
			 dyr[ii]+=r; 
			 dyi[ii]+=x;  
			 continue;  
		 }