/********************************************************************/
/*                                                                  */
/*      QHCSIM.C -- Flow, Pressure, & Water Quality Simulation      */
/*                  for Water Distribution Pipe Networks            */
/*                                                                  */
/*      Written by: L. Rossman                                      */
/*                  US EPA - RREL                                   */
/*                  10/16/91                                        */
/********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*  MACROS.H  include file from
 *
 *  COMMON C FUNCTIONS
 *
 *  several conversion macros to be used in place of functions
 *
 */

#define INT(x)   ((int)(x))                   /* integer portion of x  */
#define FRAC(x)  ((x)-(int)(x))               /* fractional part of x  */
#define ABS(x)   (((x)<0) ? -(x) : (x))       /* absolute value of x   */
#define MIN(x,y) (((x)<=(y)) ? (x) : (y))     /* minimum of x and y    */
#define MAX(x,y) (((x)>=(y)) ? (x) : (y))     /* maximum of x and y    */
#define ROUND(x) (((x)>=0) ? (x)+.5 : (x)-.5) /* round-off of x        */
#define MOD(x,y) ((x)%(y))                    /* x modulus y           */
#define SQR(x)   ((x)*(x))                    /* x-squared             */
#define SGN(x)   (((x)<0) ? (-1) : (1))       /* sign of x             */

/*--------------------------------------------------------------------*/
/*      Global Constants                                              */
/*--------------------------------------------------------------------*/
#define  MAXN      100          /* Max. dimension of Jacobian matrix */
#define  MAXNODES  100          /* Max. # nodes */
#define  MAXPIPES  200          /* Max. # pipes */
#define  MAXTANKS  10           /* Max. # tanks */
#define  MAXTEXTS  50           /* Max. # text lines on map */
#define  MAXPUMPS  50           /* Max. # pump control rules */
#define  MAXID     10000        /* Highest node/pipe ID number */
#define  MAXITER   20           /* Max. # Newton-Raphson iterations */
#define  MAXHRS    24           /* Max. hrs. per demand pattern */
#define  MAXSEGS   50           /* Max. segments per pipe */
#define  ACCEL     1.85         /* Factor used in Newton-Raphson */
#define  FLOW0     1            /* Flow used to initialize network */
#define  HTOL      .01          /* Default zero head tolerance (ft) */
#define  QTOL      5            /* Default zero flow tolerance (gpm) */
#define  TTOL      5            /* Default zero time tolerance (min) */
#define  UNKNOWN   -999999
#define  HUGE      1.E10
#define  HUGE_INT  999999999
#define  PI        3.14
#define  GPMperCFS 448.4
#define  LperFT3   28.3

/*--------------------------------------------------------------------*/
/*      Global Data Structures                                        */
/*--------------------------------------------------------------------*/
typedef float BIG;

typedef struct
{
   int    ID;                   /* Node ID number */
   int    Tank;                 /* Tank index if node is fixed grade */
   float  X,Y;                  /* Map coordinates (not used ) */
   float  E;                    /* Elevation of node */
   float  C0;                   /* Initial concen. at node */
   float  Qbase;                /* Base demand flow from node */
   float  Cs;                   /* Source concen. */
   float  Q;                    /* Actual demand flow from node */
   float  H;                    /* Head at node */
   float  C;                    /* Concentration at node */
}  Snode;

typedef struct
{
   int    ID;                   /* Pipe ID number */
   int    N1;                   /* Start node index of pipe */
   int    N2;                   /* End node index of pipe */
   float  D;                    /* Diam. of pipe (in) */
   float  L;                    /* Length of pipe (ft) */
   float  C;                    /* Hazen Williams Coeff. */
   float  K1;                   /* Growth coeff. within pipe */
   float  K;                    /* Head loss coeff. for pipe */
   float  KK;                   /* Discharge coeff. for pipe */
   float  V;                    /* Volume of pipe */
   float  Q;                    /* Flow in pipe */
}  Spipe;

typedef struct
{
   int    Node;                 /* Node index of tank */
   float  A;                    /* Tank area */
   float  Hmin;                 /* Minimum tank height */
   float  Hmax;                 /* Maximum tank height */
   float  H0;                   /* Initial height of water in tank */
   float  K1;                   /* Growth coeff. within tank */
   float  V;                    /* Tank volume */
   float  C;                    /* Tank concentration */
}  Stank;

typedef struct
{
   int    Pump;                 /* Pump being controlled */
   int    Pipe;                 /* Pipe index of pump */
   int    Ontime;               /* Time period turned on */
   int    Offtime;              /* Time period turned off */
   int    Tank;                 /* Tank controlling pump */
   float  Onlevel;              /* Level turning pump on */
   float  Offlevel;             /* Level turning pump off */
}  Spump;

/*--------------------------------------------------------------------*/
/*        Global Variables                                            */
/*--------------------------------------------------------------------*/
FILE     *DataFile;             /* Data file handle */
FILE     *OutFile;              /* Output file handle */
char     Title[81];             /* Problem title */
int      Nnodes,                /* Number of network nodes */
         Npipes,                /* Number of network pipes */
         Npumps,                /* Number of pump control rules */
         Ntanks,                /* Number of network tanks */
         Nhsteps,               /* Number of hydraulic time steps */
         Htime,                 /* Hydraulic time period */
         Hdur;                  /* Simulation duration */
Spipe    Pipe[MAXPIPES];        /* Pipe data */
Snode    Node[MAXNODES];        /* Node data */
Stank    Tank[MAXTANKS];        /* Tank data */
Spump    Pump[MAXPUMPS];        /* Pump control rules */
float    Demand[MAXHRS];        /* Demand schedule */
float    Htol;                  /* Head tolerance (ft) */
float    Qtol;                  /* Flow tolerance (cfs) */
float    Ttol;                  /* Time tolerance (sec) */
float    Krate;                 /* Universal reaction rate coeff. */
BIG      Hstep,                 /* Hydraulic time step (sec) */
         Rstep,                 /* Reporting time step (sec) */
         Cstep;                 /* Water quality time step (sec) */
BIG      FlowVol[MAXPIPES];     /* Flow volume per time step */
BIG      SegVol[MAXPIPES];      /* Pipe segment volume */
BIG      *SegMass[MAXPIPES];    /* Pointer to pipe segment masses */
int      Nseg[MAXPIPES];        /* Number of segments per pipe */
char     Dir[MAXPIPES];         /* Flow direction (+/-) */

/*--------------------------------------------------------------------*/
/*       Function Prototypes                                          */
/*--------------------------------------------------------------------*/
void    setdefaults(void);              /* Initializes default values */
void    getdata(void);                  /* Gets network data          */
int     qhcsim(void);                   /* Performs total simulation  */
void    initparams(void);               /* Initializes nodes & pipes  */
void    demands(void);                  /* Computes current demands   */
int     controlpumps(void);             /* Turns pumps on/off         */
void    output();                       /* Writes output report       */

int     hsim(int, int *);               /* Performs hydraulic simulation */
void    head0(void);                    /* Initializes heads          */
int     hsolve(int *);                  /* Solves network for heads   */
float   getqc(int,int,int);             /* Gets a discharge coeff.    */
void    tankheights(void);              /* Computes new tank levels   */

/* These functions solve a set of simultaneous linear equations */
void    gefa(float [][], int, int [], int *);
void    gesl(float [][], int, int [], float []);
void    simul(float [][], int, int [], float [], int *);

int     csim(void);                     /* Performs quality simulation */
BIG     mintraveltime(void);            /* Finds minimum travel time   */
int     segmentpipes(BIG);              /* Divides pipe into segments  */
void    interpmass(int);                /* Finds initial mass distribution */
void    shiftdown(int, int,             /* Redistributes mass within a */
                  int, int,             /* pipe when the number of     */
                  BIG, BIG,             /* segments decreases          */
                  BIG *, BIG *);
void    shiftup(int, int,               /* Redistributes mass within a */
                int, int,               /* pipe when the number of     */
                BIG, BIG,               /* segments increases          */
                BIG *, BIG *);
void    reactmass(BIG);                 /* Reacts mass within pipes    */
void    advectmass(BIG, BIG);           /* Computes mass flow along pipes */


void  main(int argc, char *argv[])
/*------------------------------------------------------------------*/
/*      Main Program Segment                                        */
/*------------------------------------------------------------------*/
{
   int i;
   if (argc < 3)
   {
      printf("\nCorrect syntax is:");
      printf(" qhcsim <input file> <output file>\n");
      exit(0);
   }
   if ((DataFile = fopen(argv[1],"r")) == NULL)
   {
      printf("\nCannot open input file %s.\n",argv[1]);
      exit(0);
   }
   if ((OutFile = fopen(argv[2],"w")) == NULL)
   {
      printf("\nCannot open output file %s.\n",argv[2]);
      exit(0);
   }
   setdefaults();
   getdata();
   do
   {
      printf("\nSimulating time period %d",Htime+1);
      i = qhcsim();
      if (i == 1) output();
   }  while (i == 1 && Htime < Hdur);
   fclose(OutFile);
   printf("\n");
}


void  setdefaults()
/*----------------------------------------------------------------*/
/*      Assigns default values to global variables                */
/*----------------------------------------------------------------*/
{
   int i;
   strcpy(Title,"Flow, Pressure, Quality Simulation");
   Htol  = HTOL;                /* Zero head tolerance */
   Qtol  = QTOL/GPMperCFS;      /* Zero flow tolerance */
   Ttol  = TTOL*60;             /* Zero travel time tolerance */
   Krate = 0;                   /* No universal reaction rate */
   Hstep = 3600;                /* Hydraulic time step = 1 hr */
   Cstep = 0;                   /* No preset qual. time step */
   Rstep = Hstep;
   Nhsteps = MAXHRS;            /* Demand schedule is full */
   for (i=0; i<MAXHRS; i++)
      Demand[i] = 1;            /* All demands = base demands */
}


void getdata(void)
/*----------------------------------------------------------------*/
/*      Reads in network data from disk file                      */
/*----------------------------------------------------------------*/
{
   char  line[80];
   char  label[40];
   char  key[15];
   char  *item;
   int   i1,i2,i3,i4,mpumps;
   float x1,x2,x3,x4,x5,x6,t1,t2;
   unsigned int *ID;

   ID = (unsigned int *) calloc(MAXID,sizeof(unsigned int));
   if (ID == NULL) return;

   Nnodes = 0;
   Npipes = 0;
   mpumps = 0;
   Npumps = 0;
   Ntanks = 0;
   Hdur   = 1;
   for(;;)
   {
      if (fgets(line,80,DataFile) == NULL) break;
      if (strstr(line,"TITLE") != NULL)
         strcpy(Title,line);
      else if (strstr(line,"NODE") != NULL && Nnodes < MAXNODES)
      {
         if (sscanf(line,"%s %d %f %f %f %f %f %f",
             key,&i1,&x1,&x2,&x3,&x4,&x5,&x6) == 8)
         {
            Nnodes++;
            Node[Nnodes].ID     = i1;
            Node[Nnodes].X      = x1;
            Node[Nnodes].Y      = x2;
            Node[Nnodes].E      = x3;
            Node[Nnodes].C0     = x4*LperFT3;
            Node[Nnodes].Qbase  = x5/GPMperCFS;
            Node[Nnodes].Cs     = x6*LperFT3;
            Node[Nnodes].Tank   = 0;
            *(ID+i1) = Nnodes;
         }
      }
      else if (strstr(line,"PIPE") != NULL && Npipes < MAXPIPES)
      {
         if (sscanf(line,"%s %d %d %d %f %f %f %f",
                    key,&i1,&i2,&i3,&x1,&x2,&x3,&x4) == 8)
         {
            Npipes++;
            Pipe[Npipes].ID  = i1;
            Pipe[Npipes].N1  = *(ID+i2);
            Pipe[Npipes].N2  = *(ID+i3);
            Pipe[Npipes].D   = x1;
            Pipe[Npipes].L   = x2;
            Pipe[Npipes].C   = x3;
            Pipe[Npipes].K1  = x4;
         }
      }
      else if (strstr(line,"PUMP") != NULL && Npipes < MAXPIPES)
      {
         if (sscanf(line,"%s %d %d %d %d %f %f %f",
                    key,&i1,&i2,&i3,&i4,&x1,&x2,&x3) == 8)
         {
            Npipes++;
            mpumps++;
            Pipe[Npipes].ID   = i1;
            Pipe[Npipes].N1   = *(ID+i2);
            Pipe[Npipes].N2   = *(ID+i3);
            Pipe[Npipes].K1   = 0;
            Pipe[Npipes].K    = i4;
            Pipe[Npipes].KK   = 1.e-10;
            Pipe[Npipes].V    = 0;
            Pipe[Npipes].D    = -x2/2/x1;
            Pipe[Npipes].L    = (x2*x2 - 4*x1*x3)/(4*x1*x1);
            Pipe[Npipes].C    = 1/x1;
         }
      }
      else if (strstr(line,"CONTROL") != NULL && Npumps < MAXPUMPS)
      {
         if (sscanf(line,"%s %d %d %d",key,&i1,&i2,&i3) == 4)
         {
            Npumps++;
            Pump[Npumps].Pump    = i1;
            Pump[Npumps].Pipe    = 0;
            Pump[Npumps].Ontime  = i2;
            Pump[Npumps].Offtime = i3;
            Pump[Npumps].Tank     = 0;
            Pump[Npumps].Onlevel  = 0;
            Pump[Npumps].Offlevel = 0;
         }
         else if (sscanf(line,"%s %d %d %f %f",key,&i1,&i2,&x1,&x2) == 5)
         {
            Npumps++;
            Pump[Npumps].Pump     = i1;
            Pump[Npumps].Pipe     = 0;
            Pump[Npumps].Ontime   = 0;
            Pump[Npumps].Offtime  = 0;
            Pump[Npumps].Tank     = i2;
            Pump[Npumps].Onlevel  = x1;
            Pump[Npumps].Offlevel = x2;
         }
      }
      else if (strstr(line,"TANK") != NULL && Ntanks < MAXTANKS)
      {
         if (sscanf(line,"%s %d %f %f %f %f %f",
             key,&i1,&x1,&x2,&x3,&x4,&x5) == 7)
         {
            Ntanks++;
            Tank[Ntanks].Node   = *(ID+i1);
            Tank[Ntanks].A      = x1;
            Tank[Ntanks].Hmin   = x2;
            Tank[Ntanks].Hmax   = x3;
            Tank[Ntanks].H0     = x4;
            Tank[Ntanks].K1     = x5;
            Node[Tank[Ntanks].Node].Tank = Ntanks;
         }
      }
      else if (strstr(line,"USAGE") != NULL)
      {
         item = strtok(line," ");        /* skip key word */
         item = strtok(NULL," ");        /* get start period */
         if (item != NULL                /* check its validity */
            && (i1 = atoi(item)) != 0
            &&  i1 <= MAXHRS)
         {
            do
            {
               item = strtok(NULL," ");  /* get next usage factor */
               if (item != NULL          /* check its validity */
                  && (x1 = atof(item)) != 0
                  &&  i1 <= MAXHRS)
                  Demand[i1-1] = x1;     /* record usage factor */
               i1++;                     /* move to next period */
            }  while (item != NULL);
         }
      }
      else if (strstr(line,"DURATION") != NULL)
      {
         if (sscanf(line,"%s %d",key,&i1) == 2) Hdur = i1;
      }
      else if (strstr(line,"REACT") != NULL)
      {
         if (sscanf(line,"%s %f",key,&x1) == 2) Krate = x1;
      }
      else if (strstr(line,"TOLERANCE") != NULL)
      {
         if (sscanf(line,"%s %f %f %f",key,&x1,&x2,&x3) == 4)
         {
            Qtol = x1/GPMperCFS;
            Htol = x2;
            Ttol = x3*60;
         }
      }
   }
   free(ID);
   fclose(DataFile);
   Htime = 0;
   printf("\n%s",Title);
   printf("\n   Number of Pipes = %3d",Npipes);
   printf("\n   Number of Nodes = %3d",Nnodes);
   printf("\n   Number of Tanks = %3d",Ntanks);
   printf("\n   Number of Pumps = %3d",mpumps);
   printf("\n   Duration (hrs)  = %3d",Hdur);
   printf("\n");
}


int  qhcsim()
/*--------------------------------------------------------------------*/
/*      Main Simulation Procedure                                     */
/*--------------------------------------------------------------------*/
{
   int   iter,                          /* Iteration counter */
         reset;                         /* Network reset indicator */

   if (Nnodes == 0 || Npipes == 0)
   {
      printf("\nThere are no network data");
      printf(" -- simulation ended.");
      return(0);
   }
   if (Htime == 0) initparams();        /* Initialize parameters */
   demands();                           /* Find demands in time period */
   reset = controlpumps();              /* Turn pumps on/off */
   if (Htime == 0) reset = 1;
   if (!hsim(reset,&iter))               /* Solve network equations */
   {
      printf("\n   Cannot solve network hydraulic equations");
      printf(" -- simulation ended.");
      return(0);
   }
   printf("\n   Hydraulic solution reached in %2d iterations.",iter);
   if (!csim())                          /* Transport mass thru network */
   {
      printf("\n   Cannot solve mass transport equations");
      printf(" -- simulation ended.");
      return(0);
   }
   printf("\n   Mass transport completed.");
   Htime = Htime + 1;
   return(1);
}


void  initparams()
/*----------------------------------------------------------------*/
/*      Initialize pipe & tank parameters                         */
/*----------------------------------------------------------------*/
{
   int   i,k,n;
   float t1,t2;

   for (k=1; k<=Npipes; k++)
   {
      if (Pipe[k].D < 0)
      {
         if (Npumps > 0)
           for (i=1; i<=Npumps; i++)
            if (Pump[i].Pump == Pipe[k].ID) Pump[i].Pipe = k;
      }
      else
      {
         t1 = pow(Pipe[k].C,1.85);
         t2 = Pipe[k].D/12.;
         t2 = pow(t2,4.87);
         Pipe[k].K = 4.72*Pipe[k].L/t1/t2;
         Pipe[k].KK = pow((1/Pipe[k].K),.54);
         Pipe[k].V = PI*Pipe[k].D*Pipe[k].D/4/144*Pipe[k].L;
      }
   }
   for (i=1; i<=Nnodes; i++)
   {
      Node[i].C = Node[i].C0;
      n = Node[i].Tank;
      if (n > 0)
      {
         Node[i].C = Node[i].Cs;
         Tank[n].C = Node[i].C;
         Tank[n].V = Tank[n].H0*Tank[n].A;
      }
   }
}


void  demands()
/*--------------------------------------------------------------------*/
/*      Computes demands at nodes during current time period          */
/*--------------------------------------------------------------------*/
{
   int i;
   for (i=1; i<=Nnodes; i++)
   {
      if (Node[i].Qbase <= 0) Node[i].Q = Node[i].Qbase;
      else
         Node[i].Q = Node[i].Qbase*Demand[Htime%Nhsteps];
   }
}


int  controlpumps()
/*-----------------------------------------------------------------*/
/*      Determines which pumps are on or off                       */
/*-----------------------------------------------------------------*/
{
   int   i, k, n, t, status, reset = 0;
   float z;

   if (Npumps == 0) return(0);
   t = Htime + 1;

/* Examine each pump control condition */

   for (i=1; i<=Npumps; i++)
   {

   /* Make sure that pipe is a pump */

      k = Pump[i].Pipe;
      if (Pipe[k].D >= 0) continue;
      status = (int) Pipe[k].K;

   /* Pump controlled by tank level */

      if ((n = Pump[i].Tank) > 0)
      {
         z = Node[n].H - Node[n].E;
         if (Node[n].Tank == 0) z = 2.308*z;
         if (z <= Pump[i].Onlevel)  status = 1;
         if (z >= Pump[i].Offlevel) status = 0;
      }

   /* Pump is time-controlled */

      else
      {
         if (t == Pump[i].Ontime)  status = 1;
         if (t == Pump[i].Offtime) status = 0;
      }

   /* Update pump status */
      if (status != (int) Pipe[k].K)
      {
         Pipe[k].K = status;
         if (status == 1)  printf("\n   Pump %3d turned on.",k);
         else              printf("\n   Pump %3d turned off.",k);
         reset = 1;
      }
   }
   return(reset);
}


void  output()
/*-----------------------------------------------------------------*/
/*      Writes out simulation results for current time step        */
/*-----------------------------------------------------------------*/
{
   static char ititle = 0;
   int    i,j,n1,n2;
   float  c,p,q,t,v;

   if (ititle == 0)
   {
      fprintf(OutFile,"\n    %s",Title);
      ititle = 1;
   }
   fprintf(OutFile,"\n\n    Node Data at Hour %3.0f",Htime*Hstep/3600);
   fprintf(OutFile,
   "\n    ---------------------------------------------------------------");
   fprintf(OutFile,
   "\n    Node   Elev.    Demand     Head   Pressure   Concen.");
   fprintf(OutFile,
   "\n     No.    Ft.       GPM       Ft.      PSI       M/L  ");
   fprintf(OutFile,
   "\n    ---------------------------------------------------------------");
   for (i=1; i<=Nnodes; i++)
   {
      q = Node[i].Q*GPMperCFS;
      p = (Node[i].H-Node[i].E)*.433;
      c = Node[i].C/LperFT3;
      fprintf(OutFile,"\n    %4d  %6.1f    %6.0f ",
      Node[i].ID,Node[i].E,q);
      if ((j=Node[i].Tank) > 0)
      {
         if (Tank[j].A == 0)
         fprintf(OutFile,"                     %7.1f   Reservoir",c);
         else
         fprintf(OutFile,"  %6.1f    %6.1f   %7.1f   Tank",
         Node[i].H,p,c);
      }
      else
      fprintf(OutFile,"  %6.1f    %6.1f   %7.1f",Node[i].H,p,c);
   }
   fprintf(OutFile,
   "\n    ---------------------------------------------------------------");
   fprintf(OutFile,"\n\n    Pipe Data at Hour %3.0f",Htime*Hstep/3600);
   fprintf(OutFile,
   "\n    ---------------------------------------------------------------");
   fprintf(OutFile,
   "\n    Pipe    Nodes     Diam.   Length    Flow  Velocity  Travel Time");
   fprintf(OutFile,
   "\n     No.   From  To    In.      Ft.      GPM   ft/sec        Min   ");
   fprintf(OutFile,
   "\n    ---------------------------------------------------------------");
   for (i=1; i<=Npipes; i++)
   {
      if (Pipe[i].D < 0)
      fprintf(OutFile,"\n    %4d   %3d  %3d    PUMP            %5.0f",
      Pipe[i].ID,Node[Pipe[i].N1].ID,Node[Pipe[i].N2].ID,
      Pipe[i].Q*GPMperCFS);
      else
      {
         n1 = Pipe[i].N1;
         n2 = Pipe[i].N2;
         if (Pipe[i].Q < 0)
         {
            n1 = Pipe[i].N2;
            n2 = Pipe[i].N1;
         }
         q = ABS(Pipe[i].Q);
         v = q/Pipe[i].D/Pipe[i].D/3.14*4*144;
         t = 0;
         if (v > 0) t = Pipe[i].L/v/60;
         q = q*GPMperCFS;
         fprintf(OutFile,
         "\n    %4d   %3d  %3d   %5.1f   %6.0f   %5.0f    %4.1f     %7.0f",
         Pipe[i].ID,Node[n1].ID,Node[n2].ID,Pipe[i].D,Pipe[i].L,q,v,t);
      }
   }
fprintf(OutFile,
"\n    ---------------------------------------------------------------\n");
}

/*----------------------------------------------------------------------*/
/*                                                                      */
/*                Network Hydraulic Simulator                           */
/*                                                                      */
/*----------------------------------------------------------------------*/

int   hsim(int reset, int *iter)
/*-------------------------------------------------------------------*/
/*      Hydraulic Simulator                                          */
/*-------------------------------------------------------------------*/
{
   if (reset == 1)  head0();            /* Compute initial heads */
   if (hsolve(iter) > 0)                /* Solve network equations */
   {                                    /* If ill-conditioned matrix, */
      head0();                          /* re-compute initial heads */
      if (hsolve(iter) > 0) return(0);  /* & try to solve once again */
   }
   tankheights();                       /* Adjust tank heights */
   return(1);
}


void head0()
/*--------------------------------------------------------------------*/
/*      Determines initial heads in network when each pipe            */
/*      carries a flow of FLOW0                                       */
/*--------------------------------------------------------------------*/
{
   int   i,k,n1,n2,nh;
   float a,b,c,havg,dh;

/* Initialize heads & compute average fixed grade node elevation */

   havg = 0.;
   nh = 0;
   for (i=1; i<=Nnodes; i++)
   {
      Node[i].H = UNKNOWN;              /* Head unknown if node not at */
      if (Node[i].Tank > 0)             /*   a fixed grade elevation */
      {
         Node[i].H = Node[i].E + Tank[Node[i].Tank].H0;
         havg += Node[i].H;
         nh++;
      }
   }
   havg = havg/nh;

/* For each pipe in network */

   for (k=1; k<=Npipes; k++)
   {
      n1 = Pipe[k].N1;
      n2 = Pipe[k].N2;
      if (Node[n1].H != UNKNOWN &&              /* Skip if heads at ends */
          Node[n2].H != UNKNOWN) continue;      /*  are known */
      if (Pipe[k].D < 0)                        /* Pipe is a pump */
      {
         if (Pipe[k].K == 0) continue;          /* Skip if pump is offline */
         else
         {
            a = 1/Pipe[k].C;
            b = -2*a*Pipe[k].D;
            c = (b*b-4*a*a*Pipe[k].L)/(4*a);
            dh = a*FLOW0*FLOW0 + b*FLOW0 + c;
         }
      }
      else
         dh = Pipe[k].K*pow(FLOW0,1.85);        /* dh = Head loss */
      if (Node[n1].H != UNKNOWN)                /* Adjust head at end of */
      {                                         /*  pipe with unknown head */
         Node[n2].H = Node[n1].H - SGN(FLOW0)*dh;
      }
      else if (Node[n2].H != UNKNOWN)
      {
         Node[n1].H = Node[n2].H + SGN(FLOW0)*dh;
      }
      else
      {
         Node[n1].H = havg;                       /* Both ends have unknown */
         Node[n2].H = Node[n1].H - SGN(FLOW0)*dh; /*  head so assign one end */
      }                                           /*  the value havg */
   }
}



void  tankheights()
/*----------------------------------------------------------------------*/
/*      Computes new storage tank heights at the end of the             */
/*              current time period                                     */
/*----------------------------------------------------------------------*/
{
   int i,n;
   if (Ntanks == 0) return;
   for (i=1; i<=Ntanks; i++)
   {
      if (Tank[i].A == 0) continue;
      n = Tank[i].Node;
      Node[n].H += Node[n].Q*Hstep/Tank[i].A;
      Node[n].H = MAX(Node[n].H, Node[n].E+Tank[i].Hmin);
      Node[n].H = MIN(Node[n].H, Node[n].E+Tank[i].Hmax);
   }
}


int  hsolve(int *iter)
/*----------------------------------------------------------------------*/
/*      Solves network nodal continuity equations in terms of heads.    */
/*      iter = iteration counter                                        */
/*----------------------------------------------------------------------*/
{
   float F[MAXN],               /* Value of continuity equation */
         J[MAXN][MAXN];         /* Jacobian of continuity equation */
   int   ipvt[MAXN];            /* Vector needed for simul routine */
   float a,b,f,                 /* Temporary variables */
         dh,                    /* Head loss */
         qc,                    /* Discharge coefficient (q = qc*dh^.54) */
         maxdiff;               /* Max. change in heads */
   int   i,j,k,n1,n2,           /* Node & pipe indices */
         err;                   /* Error flag from simul routine */

/* Repeat until convergence acheived */

   *iter = 1;
   for (;;)
   {

/* Find flow in each pipe & its contribution towards continuity  */
/* eqn. for nodes at each end of pipe.                           */
/* NOTE: Continuity Eqn. at node i is:                           */
/*       F[i] = Sum(flow out of i) - Sum(flow into i) + Demand   */

      for (i=1; i<=Nnodes; i++) F[i] = 0;
      for (k=1; k<=Npipes; k++)
      {
         n1 = Pipe[k].N1;
         n2 = Pipe[k].N2;
         qc = getqc(k,n1,n2);
         dh = Node[n1].H - Node[n2].H;
         if (Pipe[k].D < 0)                     /* Pipe is a pump */
         {
            if (Pipe[k].K == 1 && dh < 0)
            Pipe[k].Q = Pipe[k].D + sqrt(Pipe[k].L - Pipe[k].C*dh);
            else Pipe[k].Q = 0;
         }
         else                                   /* Pipe not a pump */
            Pipe[k].Q = SGN(dh)*qc*pow(ABS(dh),0.54);
         if (ABS(Pipe[k].Q) < Qtol) Pipe[k].Q = 0;
         F[n1] += Pipe[k].Q;
         F[n2] -= Pipe[k].Q;
      }

/* For each node, add demand to F & reverse sign of F.  */
/* (NOTE: nodes with fixed grades (tanks) have F = 0.)  */

      for (i=1; i<=Nnodes; i++)
      {
         if (Node[i].Tank > 0)                  /* Node is a tank */
         {
            Node[i].Q = -F[i];
            F[i] = 0;
         }
         else  F[i] = -(F[i] + Node[i].Q);      /* Node not a tank */
      }

/* Initialize entries of Jacobian matrix (dF/dH) */

      for (i=1; i<=Nnodes; i++)
      for (j=1; j<=Nnodes; j++) J[i][j] = 0;

/* For each pipe, determine its contribution to the Jacobian */
/* elements associated with the nodes at each end  of pipe.  */

      for (k=1; k<=Npipes; k++)
      {
         n1 = Pipe[k].N1;
         n2 = Pipe[k].N2;
         qc = getqc(k,n1,n2);
         dh = Node[n1].H-Node[n2].H;
         if (dh != 0)
         {
            if (Pipe[k].D < 0)                          /* Pipe is a pump */
            {
               if (Pipe[k].K == 1 && dh < 0)
               f = -Pipe[k].C/sqrt(Pipe[k].L - Pipe[k].C*dh)/2;
               else f = 0;
            }
            else  f = -qc*pow(ABS(dh),-.46)/.54;        /* Pipe not a pump */
            J[n1][n2] = f;
            J[n2][n1] = f;
            J[n1][n1] -= f;
            J[n2][n2] -= f;
         }
      }

/* For fixed grade nodes (tanks), set off-diagonal */
/* elements of J to 0, & diagonal element to 1.    */

      for (i=1; i<=Nnodes; i++)
      {
         if (Node[i].Tank > 0)
         {
            for (j=1; j<=Nnodes; j++) J[i][j] = 0;
            J[i][i] = 1;
         }
         else if (J[i][i] == 0) J[i][i] = 1;
      }

/* Invoke simul routine to solve Jh = -F     */
/* where solution h is returned in vector F. */

      simul(J,Nnodes,ipvt,F,&err);
      if (err > 0) return(err);

/* Update heads at each node using formula:    */
/*   H = H + ACCEL*F                           */
/* where ACCEL = overrelaxation parameter,     */

      maxdiff = 0;
      if (*iter > 1 && *iter < 11) a = ACCEL;
      else           a = 1;
      for (i=1; i<=Nnodes; i++)
      {
         if (ABS(F[i]) > maxdiff) maxdiff = ABS(F[i]);
         if (Node[i].Tank == 0) Node[i].H += a*F[i];
      }

/* Increment iteration counter & repeat the process */

      if (maxdiff <= Htol) return(0);
      if (*iter == MAXITER) return(0);
      *iter = *iter + 1;
   }
}


float  getqc(int k, int n1, int n2)
/*----------------------------------------------------------------------*/
/*   Gets discharge coefficient (qc) for pipe k between nodes n1 & n2   */
/*----------------------------------------------------------------------*/
{
   int i;

/* Check if start node is a tank */

   if ((i=Node[n1].Tank) > 0 && Tank[i].A > 0)
   {

   /* qc = 0 if tank full & end node at higher head */

      if ((Node[n1].H-Node[n1].E) >= Tank[i].Hmax
          && Node[n2].H > Node[n1].H) return(0);

   /* qc = 0 if tank empty & end node at lower head */

      if ((Node[n1].H-Node[n1].E) <= Tank[i].Hmin
          && Node[n1].H > Node[n2].H) return(0);
   }

/* Check if end node is a tank */

   if ((i=Node[n2].Tank) > 0 && Tank[i].A > 0)
   {

   /* qc = 0 if tank full & start node at higher head */

      if ((Node[n2].H-Node[n2].E) >= Tank[i].Hmax
          && Node[n1].H > Node[n2].H) return(0);

   /* qc = 0 if tank empty & start node at lower head */

      if ((Node[n2].H-Node[n2].E) <= Tank[i].Hmin
          && Node[n2].H > Node[n1].H) return(0);
   }

/* Otherwise return with pipe's normal qc */

   return(Pipe[k].KK);
}


/*----------------------------------------------------------------------*/
/*                                                                      */
/*                Simultaneous Equation Solver                          */
/*                                                                      */
/*----------------------------------------------------------------------*/

void  gefa(a,n,ipvt,err)
/*-----------------------------------------------------------------------
     Gaussian Elimination FActorization of a matrix
--------------------------------------------------------------------------*/
float a[][MAXN];
int    n;
int    ipvt[];
int    *err;
{
   int i,j,k,l;
   float t;

   *err = 0;
   if (n == 1) goto RETURN;
   for (k=1; k<n; k++)               /* find l = pivot index */
   {
      l = k;
      t = a[k][k];
      for (j=k+1; j<=n; j++)
      {
         if (ABS(a[j][k]) > t)
         {
            l = j;
            t = ABS(a[j][k]);
         }
      }
      ipvt[k] = l;
      if (a[l][k] == 0.0) *err = k;
      else
      {                             /* Interchange pivot element */
         t = a[l][k];
         a[l][k] = a[k][k];
         a[k][k] = t;
         t = -1./a[k][k];           /* Compute multiplier */
         for (j=k+1; j<=n; j++) a[j][k] = t*a[j][k];
         for (j=k+1; j<=n; j++)     /* Row elimination with */
         {                          /*  column indexing */
            t = a[l][j];
            if (l != k)
            {
               a[l][j] = a[k][j];
               a[k][j] = t;
            }
            for (i=k+1; i<=n; i++) a[i][j] = a[i][j] + t*a[i][k];
         }  /* next j */
      }
   }        /* next k */
   RETURN:
   ipvt[n] = n;
   if (a[n][n] == 0.0) *err = n;
}  /* end gefa */


void  gesl(a,n,ipvt,b)
/*------------------------------------------------------------------
   Gaussian Elimination solution of Simultaneous Linear equations
--------------------------------------------------------------------*/
float a[][MAXN];
int    n;
int    ipvt[];
float b[];
{
   int j,k,kb,l;
   float t;

   if (n > 1)                /* First solve Ly = b */
   {
      for (k=1; k<n; k++)
      {
         l = ipvt[k];
         t = b[l];
         b[l] = b[k];
         b[k] = t;
         for (j=k+1; j<=n; j++) b[j] = b[j] + t*a[j][k];
      }
   }
                             /* Now solve Ux = y */
   for (kb=1; kb<=n; kb++)
   {
      k = n+1-kb;
      b[k] = b[k]/a[k][k];
      t = -b[k];
      for (j=1; j<k; j++) b[j] = b[j] + t*a[j][k];
   }
}     /* end gesl */


void  simul(a,n,ipvt,b,err)
/*------------------------------------------------------------------------
    Simultaneous Linear Equation Solver
-------------------------------------------------------------------------*/
float a[][MAXN];
int    n;
int    ipvt[];
float b[];
int    *err;
{
   gefa(a,n,ipvt,err);   /* L-U factorization of a */
   if (*err > 0) return;
   gesl(a,n,ipvt,b);     /* Uses factorized a to solve ax = b */
}   /* end simul */


/*----------------------------------------------------------------------*/
/*                                                                      */
/*                      Water Quality Simulator                         */
/*                                                                      */
/*----------------------------------------------------------------------*/

int  csim()
/*-------------------------------------------------------------------*/
/*      Simulates water quality conditions in pipe network           */
/*-------------------------------------------------------------------*/
{
   BIG  time,         /* Overall elapsed time  */
        ctime,        /* Elapsed time for hydraulic time step */
        rtime;        /* Elapsed time for reporting time step */
   BIG  cstep,        /* Computed water quality time step */
        dt;           /* Modified water quality time step */
   int  i,n;

   rtime = 0;
   ctime = 0;
   time  = Htime*Hstep;

/* If not set by user then water quality time step is the shortest */
/* time for water to flow through a pipe in the network            */

   cstep = mintraveltime();
   if (cstep == 0) return(1);
   if (Cstep > 0 && Cstep < cstep) cstep = Cstep;

/* Based on this time step, divide pipes into segments */
/* and assign  contaminant masses to each segment      */

   if (!segmentpipes(cstep)) return(0);

/* Repeat until elapsed time equals hydraulic time step */

   do
   {

   /* Modify water quality time step to end at hydraulic time step */
      dt = MIN(cstep,Hstep-ctime);

   /* React mass within pipes */
      reactmass(dt);

   /* Advect mass through pipes */
      advectmass(dt,cstep);

   /* Update elapsed times */
      time += dt;
      ctime += dt;
      rtime += dt;
   }  while (Hstep - ctime > Ttol);
   return(1);
}


BIG    mintraveltime(void)
/*--------------------------------------------------------------------*/
/*      Finds smallest time for water to travel through a pipe        */
/*      in the network.                                               */
/*--------------------------------------------------------------------*/
{
   BIG     t,
           cstep;               /* Min. travel time */
   int     k;
   float   q;

   cstep = HUGE_INT;

/* Find travel time in each pipe */

   for (k=1; k<=Npipes; k++)
   {
      if (Pipe[k].Q == 0) continue;
      q = ABS(Pipe[k].Q);
      t = Pipe[k].V/q;

   /* Update min. travel time */

      if (t < Ttol) continue;
      if (t < cstep) cstep = t;
   }
   if (cstep == HUGE_INT) cstep = 0;
   return(cstep);
}


int  segmentpipes(BIG cstep)
/*---------------------------------------------------------------------*/
/*      Divides each pipe into segments and determines how much        */
/*      contaminant mass is in each segment.                           */
/*---------------------------------------------------------------------*/
{
   int    i,k;
   float  qv;                   /* Flow volume during time step */
   char   dir1;                 /* Previous flow direction */
   int    nseg1;                /* Previous number of segments */
   BIG    segvol1;              /* Previous segment volume */
   BIG    *segmass1;            /* Ptr. to previous segment masses */

/* Examine each pipe in turn */

   for (k=1; k<=Npipes; k++)
   {

   /* Find volume of flow during time step */

      qv = ABS(Pipe[k].Q);
      if (qv < Qtol) qv = 0;
      qv = qv*cstep;
      qv = MIN(qv,Pipe[k].V);
      FlowVol[k] = qv;

   /* Save pipe segment data from previous time step */

      if (Htime > 0)
      {
         dir1 = Dir[k];
         nseg1 = Nseg[k];
         segvol1 = SegVol[k];
         segmass1 = (BIG *)
                    calloc(nseg1,sizeof(BIG));
         if (segmass1 == NULL) return(0);
         for (i=0; i<nseg1; i++)
           *(segmass1+i) = *(SegMass[k]+i);
         free(SegMass[k]);
      }

   /* Compute new set of pipe segments */

      Dir[k] = SGN(Pipe[k].Q);
      Nseg[k] = 0;
      if (FlowVol[k] > 0)
         Nseg[k] = MIN(MAXSEGS, Pipe[k].V/FlowVol[k]);
      if (Nseg[k] == 0) Nseg[k] = 1;
      SegVol[k] = Pipe[k].V/Nseg[k];

   /* Allocate memory for mass of each segment */

      SegMass[k] = (BIG *)
                   calloc(Nseg[k],sizeof(BIG));
      if (SegMass[k] == NULL) return(0);

   /* For the very first hydraulic time step, interpolate between */
   /* nodal concentrations to find mass within each segment.      */

      if (Htime == 0) interpmass(k);

   /* Otherwise, shift mass from old pipe segments into new ones */

      else
      {
         if (nseg1 == Nseg[k])
            for (i=0; i<nseg1; i++)
            *(SegMass[k]+i) = *(segmass1+i);
         else if (nseg1 >= Nseg[k])
            shiftdown(dir1,Dir[k],nseg1,Nseg[k],
            segvol1,SegVol[k],segmass1,SegMass[k]);
         else
            shiftup(dir1,Dir[k],nseg1,Nseg[k],
            segvol1,SegVol[k],segmass1,SegMass[k]);
         free(segmass1);
      }
   }
   return(1);
}


void  interpmass(int k)
/*-------------------------------------------------------------------*/
/*      Interpolates between concentrations at end nodes of a        */
/*      pipe to find mass contained within each volume segment.      */
/*-------------------------------------------------------------------*/
{
   int    i,n1,n2;
   BIG    m1, m2, dm;

/* Distinguish between start and end nodes of the pipe, */
/* depending on flow direction                          */

   if (Dir[k] > 0)
   {
      n1 = Pipe[k].N1;
      n2 = Pipe[k].N2;
   }
   else
   {
      n1 = Pipe[k].N2;
      n2 = Pipe[k].N1;
   }

/* Convert the nodal concentration to an equivalent pipe segment mass */

   m1 = Node[n1].C0*SegVol[k];
   m2 = Node[n2].C0*SegVol[k];

/* Interpolate between these masses to find mass within each pipe segment */

   if (Nseg[k] == 1) *(SegMass[k]) = (m1+m2)/2;
   else
   {
      dm = (m2-m1)/(Nseg[k]-1);
      for (i=0; i<Nseg[k]; i++) *(SegMass[k]+i) = m1+i*dm;
   }
}


void  shiftdown(int dir1, int dir2, int nseg1, int nseg2,
                BIG segvol1, BIG segvol2, BIG *segmass1,
                BIG *segmass2)
/*---------------------------------------------------------------------*/
/*      Re-distributes mass within a pipe going from                   */
/*      a higher to lower number of pipe segments.                     */
/*---------------------------------------------------------------------*/
{
   int    i1,i2,dir;
   BIG    segvol, f1;

   dir = dir1*dir2;
   i2 = 0;
   if (dir < 0) i2 = nseg2-1;
   segvol = 0;
   *(segmass2+i2) = 0;
   for (i1=0; i1<nseg1; i1++)
   {
      if (segvol+segvol1 < segvol2)
      {
         *(segmass2+i2) += *(segmass1+i1);
         segvol += segvol1;
      }
      else
      {
         f1 = segvol2-segvol;
         *(segmass2+i2) += *(segmass1+i1)*f1/segvol1;
         if (dir > 0)
         {
            if (i2 < nseg2-1)
            {
               i2++;
               *(segmass2+i2) = 0;
            }
         }
         else
         {
            if (i2 > 0)
            {
               i2--;
               *(segmass2+i2) = 0;
            }
         }
         segvol = segvol1-f1;
         *(segmass2+i2) += *(segmass1+i1)*segvol/segvol1;
      }
   }
}


void  shiftup(int dir1, int dir2, int nseg1, int nseg2,
              BIG segvol1, BIG segvol2,
              BIG *segmass1, BIG *segmass2)
/*---------------------------------------------------------------------*/
/*      Re-distributes mass within a pipe going from                   */
/*      a lower to higher number of pipe segments.                     */
/*---------------------------------------------------------------------*/
{
   int    i1,i2,dir;
   BIG    segvol, f1;

   dir = dir1*dir2;
   i1 = 0;
   if (dir < 0) i1 = nseg1-1;
   segvol = 0;
   for (i2=0; i2<nseg2; i2++)
   {
      if (segvol+segvol2 < segvol1)
      {
         *(segmass2+i2) = *(segmass1+i1)*segvol2/segvol1;
         segvol += segvol2;
      }
      else
      {
         f1 = segvol1-segvol;
         *(segmass2+i2) = *(segmass1+i1)*f1/segvol1;
         if (dir > 0)
              if (i1 < nseg1-1) i1++;
         else if (i1 > 0)       i1--;
         segvol = segvol2-f1;
         *(segmass2+i2) += *(segmass1+i1)*segvol/segvol1;
      }
   }
}


void  reactmass(BIG dt)
/*---------------------------------------------------------------------*/
/*      Reacts mass within each pipe segment and tank                  */
/*---------------------------------------------------------------------*/
{
   int   i,j,k;
   float m,rate,growth;

   for (k=1; k<=Npipes; k++)
   {
      rate = Pipe[k].K1;
      if (Krate != 0) rate = Krate;
      if (rate == 0) continue;
      growth = exp(rate/3600*dt);
      for (j=0; j<Nseg[k]; j++)
      {
         m = *(SegMass[k]+j);
         m = m*growth;
         *(SegMass[k]+j) = m;
      }
   }
   if (Ntanks == 0) return;
   for (i=1; i<=Ntanks; i++)
   {
      if (Tank[i].A == 0)  continue;
      rate = Tank[i].K1;
      if (Krate != 0) rate = Krate;
      if (rate == 0) continue;
      Tank[i].C *= exp(rate/3600*dt);
   }
}


void  advectmass(BIG dt, BIG cstep)
/*----------------------------------------------------------------------*/
/*      Advects mass from one pipe segment to the next, mixing          */
/*      mass & volumes at pipe junctions.                               */
/*----------------------------------------------------------------------*/
{
   int   i,j,k,n,n1,n2,nseg1;
   float qvol,f;
   BIG   vol,mass;
   BIG   summass[MAXNODES] = {0},
         sumvol[MAXNODES]  = {0};

/* Advect mass along pipe segments, accumulating mass from */
/* last segment in the node it enters                      */

   for (k=1; k<=Npipes; k++)
   {
      if (FlowVol[k] == 0) continue;
      vol = FlowVol[k];
      if (dt < cstep) vol = (vol*dt)/cstep;
      f = (float) vol/ (float) SegVol[k];
      nseg1 = Nseg[k]-1;
      for (j=nseg1; j>=0; j--)
      {
         mass = *(SegMass[k]+j)*f;
         *(SegMass[k]+j) -= mass;
         *(SegMass[k]+j) = MAX(0, *(SegMass[k]+j));

      /* Advect last segment's mass into end node */

         if (j == nseg1)
         {
            i = Pipe[k].N2;
            if (Dir[k] < 0) i = Pipe[k].N1;
            summass[i] += mass;
            sumvol[i]  += vol;
         }

      /* Advect segment's mass into next segment */

         else *(SegMass[k]+j+1) += mass;

      }  /* Next segment */
   }  /* Next pipe */

/* Compute external inputs/outputs at each node */

   for (i=1; i<=Nnodes; i++)
   {

   /* If node has input from external source then */
   /* add this to node's mass and volume          */

      if (Node[i].Q < 0)
      {
         qvol = dt*(-Node[i].Q);
         vol = qvol;
         sumvol[i]  += vol;
         if ((n=Node[i].Tank) > 0) mass = vol*Tank[n].C;
         else                      mass = vol*Node[i].Cs;
         summass[i] += mass;
      }

   }  /* Next node */

/* Transfer mass from upstream nodes to downstream nodes */
/* of pipes with zero volumes (e.g., pumps)              */

   for (k=1; k<=Npipes; k++)
   {
      if (FlowVol[k] > 0) continue;
      if (Pipe[k].Q == 0) continue;
      if (Pipe[k].Q > 0)
      {
         n1 = Pipe[k].N1;
         n2 = Pipe[k].N2;
      }
      else
      {
         n2 = Pipe[k].N1;
         n1 = Pipe[k].N2;
      }
      sumvol[n2] += sumvol[n1];
      summass[n2] += summass[n1];
   }

/*  Compute concentration at each node */

   for (i=1; i<=Nnodes; i++)
   {
      if (sumvol[i] == 0) continue;
      Node[i].C = summass[i]/sumvol[i];
   }

/* For each pipe, add mass from start node to first segment of pipe */

   for (k=1; k<=Npipes; k++)
   {
      if (FlowVol[k] == 0) continue;
      i = Pipe[k].N1;
      if (Dir[k] < 0) i = Pipe[k].N2;
      vol = FlowVol[k];
      if (dt < cstep) vol = (vol*dt)/cstep;
      *(SegMass[k]) += Node[i].C*vol;
   }  /* Next pipe */

/* For each tank, adjust volume and concentration */

   if (Ntanks == 0) return;
   for (i=1; i<=Ntanks; i++)
   {

   /* No adjustment needed if tank is a reservoir */

      if (Tank[i].A == 0) continue;

   /* If tank emptying, then reduce its volume */

      n = Tank[i].Node;
      qvol = dt*Node[n].Q;
      if (qvol < 0) Tank[i].V += qvol;

   /* If filling, then increase volume & adjust concen. */

      else
      {
         Tank[i].C = (Node[n].C*qvol + Tank[i].C*Tank[i].V)/
                     (Tank[i].V + qvol);
         Tank[i].V += qvol;
      }

   /* Replace previous nodal concen. with tank's concen. */

      Node[n].C = Tank[i].C;

   } /* Next tank */
}


**********************************************************
* To sign off, email to:   listserv@listserv.uoguelph.ca *
* In the body of the message type: signoff epanet-users  *
**********************************************************