VSP efficiency + retrieving efficiency directly (#81)

* VSP efficiency + retrieving efficiency directly

Changes to the getenergy function (in hydraulic.c) to compute the
correct efficiency of variable speed pumps according to the affinity
laws.
Changes to epanet2.h and epanet.c are for retrieving the efficiency of
pumps throught the ENgetlinkvalue function (code 16 for efficiency)

* Update epanet.c

* VSP efficiency (incl. bas and vb file)

- bas and vb file have been updated with code 16 for efficiency.
- Retrieving efficiency case has been improved.

include/epanet2.bas

@@ -52,6 +52,7 @@ Global Const EN_SETTING = 12
 Global Const EN_ENERGY = 13
 Global Const EN_LINKQUAL = 14   'ES
 Global Const EN_LINKPATTERN = 15
+Global Const EN_EFFICIENCY = 16
 
 Global Const EN_DURATION = 0      ' Time parameters
 Global Const EN_HYDSTEP = 1

include/epanet2.h

@@ -110,7 +110,8 @@ typedef enum {
   EN_SETTING      = 12,
   EN_ENERGY       = 13,
   EN_LINKQUAL     = 14,
-  EN_LINKPATTERN  = 15
+  EN_LINKPATTERN  = 15,
+  EN_EFFICIENCY   = 16
 } EN_LinkProperty;
 
 /// Time parameter codes

include/epanet2.vb

@@ -56,6 +56,7 @@ Public Const EN_SETTING = 12
 Public Const EN_ENERGY = 13
 Public Const EN_LINKQUAL = 14   'ES
 Public Const EN_LINKPATTERN = 15
+Public Const EN_EFFICIENCY = 16
 
 Public Const EN_DURATION = 0      ' Time parameters
 Public Const EN_HYDSTEP = 1

src/epanet.c

@@ -1458,7 +1458,11 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value)
          if (Link[index].Type == PUMP)
             v = (double)Pump[PUMPINDEX(index)].Upat;
          break;
-         
+
+      case EN_EFFICIENCY: 
+         getenergy(index, &a, &v);  
+         break; 
+
       default: return(251);
    }
    *value = (EN_API_FLOAT_TYPE)v;

src/hydraul.c

@@ -998,6 +998,7 @@ void  getenergy(int k, double *kw, double *eff)
 {
    int   i,j;
    double dh, q, e;
+   double q4eff; //q4eff=flow at nominal speed to compute efficiency
 
 /*** Updated 6/24/02 ***/
    /* No energy if link is closed */
@@ -1019,7 +1020,11 @@ void  getenergy(int k, double *kw, double *eff)
       j = PUMPINDEX(k);
       e = Epump;
       if ( (i = Pump[j].Ecurve) > 0)
-         e = interp(Curve[i].Npts,Curve[i].X,Curve[i].Y,q*Ucf[FLOW]);
+      { 
+         q4eff = q/LinkSetting[k]; 
+         e = interp(Curve[i].Npts,Curve[i].X, Curve[i].Y, q4eff*Ucf[FLOW]); 
+      } 
+      //   e = interp(Curve[i].Npts,Curve[i].X,Curve[i].Y,q*Ucf[FLOW]); //old line of code
       e = MIN(e, 100.0);
       e = MAX(e, 1.0);
       e /= 100.0;