/* Power Flow. */

#include <stdlib.h>
//#ifndef WINDOWS
//#include <stdio.h>
//#else
#include "pfwstdio.h"
//#endif
#include <math.h>
#include <string.h>
#include "constant.h"
#include "param.h"
#include "sparse.h"
#include "pflow.h"

#ifdef ANSIPROTO
void ErrorStop(char *Msg);
AreaData *ACFunJac(SparseMatrix *Mptr,int *val,BOOLEAN flagF,BOOLEAN flagJ,BOOLEAN flagFirst);
BOOLEAN DCFunJac(SparseMatrix *Mptr,BOOLEAN flagF,BOOLEAN flagJ);
void UpdateSVCvar(VALUETYPE cons,INDEX j);                               /* FACTS */
void SVCFunJac(SparseMatrix *Mptr,BOOLEAN flagF,BOOLEAN flagJ);          /* FACTS */
void UpdateTCSCvar(VALUETYPE cons,INDEX j);                              /* FACTS */
void TCSCFunJac(SparseMatrix *Mptr,BOOLEAN flagF,BOOLEAN flagJ);         /* FACTS */
void UpdateSTATCOMvar(VALUETYPE cons,INDEX j);                           /* FACTS */
void STATCOMFunJac(SparseMatrix *Mptr,BOOLEAN flagF,BOOLEAN flagJ);      /* FACTS */
int HFunJac(BOOLEAN flagF,BOOLEAN flagJ,AreaData *Aptr,VALUETYPE *vec);
void ACFunHes(BOOLEAN flagF,BOOLEAN flagJ);
BOOLEAN DCFunHes(BOOLEAN flagF,BOOLEAN flagJ);
void SVCFunHes(BOOLEAN flagF,BOOLEAN flagJ);       /* FACTS */
void TCSCFunHes(BOOLEAN flagF,BOOLEAN flagJ);      /* FACTS */
void STATCOMFunHes(BOOLEAN flagF,BOOLEAN flagJ);   /* FACTS */
int factorns(SparseMatrix *Mptr,double Param,IntegerVector *PartRow,IntegerVector *PartCol,
             IntegerVector *P1Row,IntegerVector *P1Col,IntegerVector *P2Row,IntegerVector *P2Col);
void repsolp(SparseMatrix *Mptr,VALUETYPE *Vptr,
             IntegerVector *PermR,IntegerVector *PermC);
VALUETYPE Norm(VALUETYPE *Vptr,INDEX N,INDEX *N1);
void WriteSolution(INDEX Iter,char *File1,char *str);
int factor(SparseMatrix *Mptr);
void UpdateACvar(VALUETYPE cons,INDEX j,BOOLEAN Limits,BOOLEAN Recover);
void UpdateDCvar(VALUETYPE cons,INDEX j,BOOLEAN Limits);
void UpdateEvector(VALUETYPE cons);
BOOLEAN CheckRlimits(void);
BOOLEAN CheckVlimits(void);
BOOLEAN CheckQlimits(void);
BOOLEAN CheckDClimits(void);
BOOLEAN ChangeSVCmode(void);     /* FACTS */
BOOLEAN ChangeTCSCmode(void);    /* FACTS */
BOOLEAN ChangeSTATCOMmode(void); /* FACTS */
BOOLEAN ChangeDCmode(void);
ACbusData *GetACbus(INDEX N);
void PrintMismatch(VALUETYPE val,INDEX j,INDEX N1);
void DeleteJac(SparseMatrix *Mptr,IntegerVector *P1Row,IntegerVector *P1Col,
               IntegerVector *P2Row,IntegerVector *P2Col);
void WriteJac(void);
int Pflow(int iter,BOOLEAN flagF,BOOLEAN flagD,BOOLEAN flagFirst);
void InitializeLoad(void);

#else
void ErrorStop();
AreaData *ACFunJac();
BOOLEAN DCFunJac();
void UpdateSVCvar();      // FACTS 
void SVCFunJac();         // FACTS 
void UpdateTCSCvar();     // FACTS 
void TCSCFunJac();        // FACTS 
void UpdateSTATCOMvar();  // FACTS 
void STATCOMFunJac();     // FACTS 
int HFunJac();
void ACFunHes();
BOOLEAN DCFunHes();
void SVCFunHes();       // FACTS 
void TCSCFunHes();      // FACTS 
void STATCOMFunHes();   // FACTS 
int factorns();
void repsolp();
VALUETYPE Norm();
void WriteSolution();
int factor();
void UpdateACvar();
void UpdateDCvar();
void UpdateEvector();
BOOLEAN CheckRlimits();
BOOLEAN CheckVlimits();
BOOLEAN CheckQlimits();
BOOLEAN CheckDClimits();
BOOLEAN ChangeSVCmode();     // FACTS 
BOOLEAN ChangeTCSCmode();    // FACTS 
BOOLEAN ChangeSTATCOMmode(); // FACTS 
BOOLEAN ChangeDCmode();
ACbusData *GetACbus();
void PrintMismatch();
void DeleteJac();
void WriteJac();
int Pflow();
void InitializeLoad();
#endif

/* ------- Global Variables ------ */
extern Data *dataPtr;
extern SparseMatrix *Jac;
extern INDEX MaxIter,Nac,NacEl,NregPQ,NregV,Ndc,Nslack,Nvolt,Narea,NacVar,Bl,
             Nsvc,Ntcsc,NtcscVar,Nstatcom; /* FACTS */
extern INDEX *ACvar;
extern VALUETYPE *dx,*dF,tol,Tol,Sn,lambda,*x0,*Dx;
extern VALUETYPE K1,K2,MaxdFi,alpha;
extern IntegerVector *NewRow,*OldRow,*NewCol,*OldCol,*RowPartition,*ColPartition;
extern IntegerVector *RowPer,*ColPer;
extern BOOLEAN Acont,PQcont,QRcont,Rcont,PQlim,Tlim,Qlim,Vlim,flagH,flagPoC,flagL,flagR;
extern INDEX *InvRowPerm,*InvColPerm;
extern BOOLEAN *MarkRowPerm,*MarkColPerm;
extern int SD0;

/* --------------- Norm ---------------------- */
#ifdef ANSIPROTO
VALUETYPE Norm(VALUETYPE *Vptr,INDEX N,INDEX *N1)
#else
VALUETYPE Norm(Vptr,N,N1)
VALUETYPE *Vptr;
INDEX N,*N1;
#endif
/* Find the norm (max. value) of a vector. */
{
  VALUETYPE val;
  INDEX i;

  val=-0.1;
  for (i=1;i<=N;i++) if (fabs(Vptr[i])>val) {val=fabs(Vptr[i]); *N1=i;}
  return(val);
}


/* ----------------- GetACbus ------------------------ */
#ifdef ANSIPROTO
ACbusData *GetACbus(INDEX N)
#else
ACbusData *GetACbus(N)
INDEX N;
#endif
{
  ACbusData *ACptr;

  for(ACptr=dataPtr->ACbus;ACptr!=NULL;ACptr=ACptr->Next) {
     if (ACptr->N==N) return(ACptr);
  }
  return(NULL);
}


/* ----------------- PrintMismatch ------------------------ */
#ifdef ANSIPROTO
void PrintMismatch(VALUETYPE val,INDEX j,INDEX N1)
#else
void PrintMismatch(val,j,N1)
VALUETYPE val;
INDEX j,N1;
#endif
{
  INDEX k,l,N,N2,N3; 
  INDEX m,n,o,N4,N5,N6; /* FACTS */
  ACbusData *ACptr;

  N2=N3=N4=N5=N6=0; l=NacVar;
  m=NacVar+11*Ndc/2;
  n=m+3*Nsvc;      /* FACTS */
  o=n+NtcscVar;    /* FACTS */
  N=o+7*Nstatcom;  /* FACTS */
  if (((!flagPoC || Jac->n1==N) && !flagH) || N1!=Jac->n1) {
     if (flagPoC && N1>N) N1=N1-N;
     else N=0;
     if (N1>NacVar && N1<=NacVar+11*Ndc/2) for (k=1; k<=Ndc/2; k++){
        N3=k;
        if (N1<=l+11) {N1=N1-l; break;}
        l=l+11;
     }
                              /* FACTS */
     else if (N1>NacVar+11*Ndc/2 && N1<=NacVar+11*Ndc/2+3*Nsvc) for (k=1;k<=Nsvc;k++){
        N4=k;
        if (N1<=m+3){N1=N1-m;break;}
        m=m+3;
     }
     else if (N1>NacVar+11*Ndc/2+3*Nsvc && N1<=NacVar+11*Ndc/2+3*Nsvc+NtcscVar) for (k=1;k<=Ntcsc;k++){
        N5=k;
        if (N1<=n+7){N1=N1-n;break;}
        n=n+7;
     }
     else if (N1>NacVar+11*Ndc/2+3*Nsvc+NtcscVar && N1<=NacVar+11*Ndc/2+3*Nsvc+NtcscVar+7*Nstatcom) for (k=1;k<=Nstatcom;k++){
        N6=k;
        if (N1<=o+7){N1=N1-o;break;}
        o=o+7;
     }
                              /* END OF FACTS */
     else for(k=1;k<=Nac;k++){
        N2=k;
        if (k+1>Nac || N1<ACvar[k+1]) {N1=N1-ACvar[k]+1; break;}
     }
  } else N=0;
  if (j!=0) fCustomPrint(stderr,"%15s","");
  fCustomPrint(stderr,"Maximum mismatch: %8.4lg  ",val);
  if (N) fCustomPrint(stderr,"PoC-");
  fCustomPrint(stderr,"Equation: %d  ",N1);
  if (N2) {
    ACptr=(ACbusData *) GetACbus(N2);
    if (ACptr!=NULL) fCustomPrint(stderr,"AC bus: %d\n",ACptr->Num);
  }
  else if(N3) fCustomPrint(stderr,"DC link: %d\n",N3);
  else if(N4) fCustomPrint(stderr,"SVC: %d\n",N4);          /* FACTS */
  else if(N5) fCustomPrint(stderr,"TCSC: %d\n",N5);         /* FACTS */
  else if(N6) fCustomPrint(stderr,"STATCOM: %d\n",N6);      /* FACTS */
  else if (flagH) fCustomPrint(stderr,"Continuation Equation\n");
  else fCustomPrint(stderr,"PoC Eigenvector Equation\n");
}

/* ------------------------ DeleteJac --------------------------------- */
#ifdef ANSIPROTO
void DeleteJac(SparseMatrix *Mptr,IntegerVector *P1Row,IntegerVector *P1Col,
               IntegerVector *P2Row,IntegerVector *P2Col)
#else
void DeleteJac(Mptr,P1Row,P1Col,P2Row,P2Col)
SparseMatrix *Mptr;
IntegerVector *P1Row,*P1Col,*P2Row,*P2Col;
#endif
{
  INDEX k;
  SparseMatrixElement *Jptr,*Jptrp;

  for (k=1;k<=Mptr->n1;k++) {
    Jptr=Mptr->RowHead[k];
    while (Jptr!=NULL) {
      Jptrp=Jptr->RowNext;