allows changing model units on-the-fly

replaces PR #140
and with special thanks to @eladsal

include/epanet2.h

@@ -513,6 +513,13 @@ extern "C" {
    */
   int  DLLEXPORT ENgetflowunits(int *code);
   
+  /**
+   @brief Sets the flow units
+   @param code Code of flow units to use
+   @return Error code
+   */
+  int  DLLEXPORT ENsetflowunits(int code);
+  
   /**
    @brief Retrieves the index of the time pattern with specified ID
    @param id String ID of the time pattern
@@ -1149,6 +1156,7 @@ extern "C" {
   int DLLEXPORT EN_getoption(EN_Project *p, EN_Option opt, EN_API_FLOAT_TYPE *value);
   int DLLEXPORT EN_gettimeparam(EN_Project *p, int code, long *value);
   int DLLEXPORT EN_getflowunits(EN_Project *p, int *code);
+  int DLLEXPORT EN_setflowunits(EN_Project *p, int code);
   int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index);
   int DLLEXPORT EN_getpatternid(EN_Project *p, int index, char *id);
   int DLLEXPORT EN_getpatternlen(EN_Project *p, int index, int *len);

src/epanet.c

@@ -223,6 +223,9 @@ int DLLEXPORT ENgettimeparam(int code, long *value) {
 int DLLEXPORT ENgetflowunits(int *code) {
   return EN_getflowunits(_defaultModel, code);
 }
+int DLLEXPORT ENsetflowunits(int code) {
+  return EN_setflowunits(_defaultModel, code);
+}
 int DLLEXPORT ENgetpatternindex(char *id, int *index) {
   return EN_getpatternindex(_defaultModel, id, index);
 }
@@ -1328,6 +1331,80 @@ int DLLEXPORT EN_getflowunits(EN_Project *p, int *code) {
   return (0);
 }
 
+int DLLEXPORT EN_setflowunits(EN_Project *p, int code) {
+  int i, j;
+  double qfactor, vfactor, hfactor, efactor, xfactor, yfactor;
+  double *Ucf = p->Ucf;
+  EN_Network *net = &p->network;
+  
+  if (!p->Openflag) { 
+    return(102);
+  }
+  
+  /* Determine unit system based on flow units */
+  qfactor = Ucf[FLOW];
+  vfactor = Ucf[VOLUME];
+  hfactor = Ucf[HEAD];
+  efactor = Ucf[ELEV];
+  
+  p->parser.Flowflag = code;
+  switch (code)
+  {
+    case LPS:          /* Liters/sec */
+    case LPM:          /* Liters/min */
+    case MLD:          /* megaliters/day  */
+    case CMH:          /* cubic meters/hr */
+    case CMD:          /* cubic meters/day */
+      p->parser.Unitsflag = SI;
+      break;
+    default:
+      p->parser.Unitsflag = US;
+      break;
+  }
+  
+  /* Revise pressure units depending on flow units */
+  if (p->parser.Unitsflag != SI) { 
+    p->parser.Pressflag = PSI;
+  }
+  else if (p->parser.Pressflag == PSI) { 
+    p->parser.Pressflag = METERS;
+  }
+  
+  initunits(p);
+  
+  //update curves
+  for (i=1; i <= net->Ncurves; i++)
+  {
+    switch (net->Curve[i].Type)
+    {
+      case V_CURVE:
+        xfactor = efactor/Ucf[ELEV];
+        yfactor = vfactor/Ucf[VOLUME];
+        break;
+      case H_CURVE:
+      case P_CURVE:
+        xfactor = qfactor/Ucf[FLOW];
+        yfactor = hfactor/Ucf[HEAD];
+        break;
+      case E_CURVE:
+        xfactor = qfactor/Ucf[FLOW];
+        yfactor = 1;
+        break;
+      default:
+        xfactor = 1;
+        yfactor = 1;
+    }
+    
+    for (j=0; j < net->Curve[i].Npts; j++)
+    {
+      net->Curve[i].X[j] = net->Curve[i].X[j]/xfactor;
+      net->Curve[i].Y[j] = net->Curve[i].Y[j]/yfactor;
+    }
+  }
+  
+  return(0);
+}
+
 int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index) {
   int i;
   *index = 0;
@@ -2757,6 +2834,7 @@ int DLLEXPORT EN_addcurve(EN_Project *p, char *id) {
 
   strcpy(tmpCur[n].ID, id);
   tmpCur[n].Npts = 1;
+  tmpCur[n].Type = -1;
   tmpCur[n].X = (double *)calloc(tmpCur[n].Npts, sizeof(double));
   tmpCur[n].Y = (double *)calloc(tmpCur[n].Npts, sizeof(double));
   if (tmpCur[n].X == NULL)

src/input2.c

@@ -400,6 +400,7 @@ int getpumpparams(EN_Project *pr)
         return (200);
       }
       curve = &net->Curve[j];
+      curve->Type = P_CURVE;
       n = curve->Npts;
       if (n == 1) {  /* Only a single h-q point supplied so use generic */
         pump->Ptype = POWER_FUNC; /* power function curve.   */

src/input3.c

@@ -198,9 +198,11 @@ int tankdata(EN_Project *pr)
     /* If volume curve supplied check it exists */
     if (n == 8) {
       t = findID(par->Tok[7], par->Curvelist);
-      if (t == NULL)
+      if (t == NULL) {
         return (202);
+      }
       vcurve = t->i;
+      net->Curve[t->i].Type = V_CURVE;
     }
   }
 
@@ -529,16 +531,16 @@ int valvedata(EN_Project *pr)
     return (202);
   if (diam <= 0.0)
     return (202);     /* Illegal diameter.*/
-  if (type == EN_GPV) /* Headloss curve for GPV */
-  {
+  if (type == EN_GPV) { /* Headloss curve for GPV */
     t = findID(par->Tok[5], par->Curvelist);
-    if (t == NULL)
+    if (t == NULL) {
       return (206);
+    }
     setting = t->i;
+    net->Curve[t->i].Type = H_CURVE;
 
     /*** Updated 9/7/00 ***/
     status = OPEN;
-
   } else if (!getfloat(par->Tok[5], &setting))
     return (202);
   if (n >= 7 && !getfloat(par->Tok[6], &lcoeff))
@@ -1436,6 +1438,7 @@ int energydata(EN_Project *pr)
       if (t == NULL)
         return (217);
       Pump[j].Ecurve = t->i;
+      net->Curve[t->i].Type = E_CURVE;
     }
     return (0);
   }