From da1253f4ea69a31162602def5084363b4538d880 Mon Sep 17 00:00:00 2001
From: Lew Rossman <LRossman@cinci.rr.com>
Date: Thu, 26 Sep 2019 11:01:01 -0400
Subject: [PATCH] Fixes bugs in pump and demand head loss gradients

---
 src/hydcoeffs.c | 20 ++++++++++----------
 src/project.c   |  3 ++-
 2 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c
index 07da0dc..ffdf98c 100644
--- a/src/hydcoeffs.c
+++ b/src/hydcoeffs.c
@@ -505,20 +505,19 @@ void demandheadloss(Project *pr, int i, double dp, double n,
         *hloss = (*hgrad) * d;
     }
 
-    // Otherwise use power head loss function
-    else
+    // Use power head loss function for demand less than full
+    else if (r < 1.0)
     {
         *hgrad = n * dp * pow(r, n - 1.0) / dfull;
         *hloss = (*hgrad) * d / n;
-      
-        // Add on barrier function for any demand above full value
-        if (r >= 1.0)
-        {
-            *hgrad += CBIG;
-            *hloss += CBIG * (d - dfull);
-        }
     }
     
+    // Use upper barrier function for demand above full value
+    else
+    {
+        *hgrad = CBIG;
+        *hloss = dp + CBIG * (d - dfull);
+    }
 }
 
 
@@ -761,7 +760,8 @@ void  pumpcoeff(Project *pr, int k)
         if (n != 1.0)
         {
              // ... use linear approx. to pump curve for small flows
-            qa = pow(hyd->RQtol / n / r, 1.0 / (n - 1.0));
+            if (pump->Ptype == CONST_HP) qa = -CBIG;
+            else qa = pow(hyd->RQtol / n / r, 1.0 / (n - 1.0));
             if (q <= qa)
             {
                 hgrad = hyd->RQtol;
diff --git a/src/project.c b/src/project.c
index 93ab40b..b585233 100644
--- a/src/project.c
+++ b/src/project.c
@@ -14,11 +14,12 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
-#include <unistd.h>
 
 //*** For the Windows SDK _tempnam function ***//
 #ifdef _WIN32
 #include <windows.h>
+#else
+#include <unistd.h >
 #endif
 
 #include "types.h"