From 4200b67c8bfca32b5bf6d6c3ff540b378707d8a8 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Wed, 16 Aug 2023 15:30:18 -0600
Subject: [PATCH] Bugfix #2644 main_v11.0 percentile (#2645)

---
 src/basic/vx_math/ptile.cc | 49 ++++++++++++++++++++++++++++++++++----
 1 file changed, 44 insertions(+), 5 deletions(-)

diff --git a/src/basic/vx_math/ptile.cc b/src/basic/vx_math/ptile.cc
index 1869320d33..6bcc178e52 100644
--- a/src/basic/vx_math/ptile.cc
+++ b/src/basic/vx_math/ptile.cc
@@ -23,6 +23,7 @@ using namespace std;
 #include "is_bad_data.h"
 
 #include "nint.h"
+#include "vx_log.h"
 
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -63,10 +64,31 @@ int index;
 double delta; 
 double p = bad_data_double;
 
+// Range check
+if ( t < 0.0 || t > 1.0 ) {
+
+   mlog << Error << "\npercentile() -> "
+        << "requested percentile value (" << t
+        << ") must be between 0 and 1!\n\n";
+
+   exit ( 1 );
+
+}
+
 if ( n > 0 ) {
+
    index = nint(floor((n - 1)*t));
-   delta = (n-1)*t - index;
-   p = (1 - delta)*ordered_array[index] + delta*ordered_array[index+1];
+
+   // Use the last value
+   if ( index == (n - 1) ) {
+      p = ordered_array[index];
+   }
+   // Interpolate linearly between two values
+   else if ( index >= 0 && index < (n - 1) ) {
+      delta = (n - 1)*t - index;
+      p = (1 - delta)*ordered_array[index] + delta*ordered_array[index + 1];
+   }
+
 }
 
 return ( p );
@@ -101,13 +123,30 @@ int index;
 float delta; 
 float p = bad_data_float;
 
+// Range check
+if ( t < 0.0 || t > 1.0 ) {
+
+   mlog << Error << "\npercentile_f() -> "
+        << "requested percentile value (" << t
+        << ") must be between 0 and 1!\n\n";
+
+   exit ( 1 );
+
+}
+
 if ( n > 0 ) {
 
    index = nint(floor((n - 1)*t));
 
-   delta = (n - 1)*t - index;
-
-   p = (1 - delta)*(ordered_array[index]) + delta*(ordered_array[index + 1]);
+   // Use the last value
+   if ( index == (n - 1) ) {
+      p = ordered_array[index];
+   }
+   // Interpolate linearly between two values
+   else if ( index >= 0 && index < (n - 1) ) {
+      delta = (n - 1)*t - index;
+      p = (1 - delta)*ordered_array[index] + delta*ordered_array[index + 1];
+   }
 
 }