From 3585997ac0a51d10008307f16c86182b93433af2 Mon Sep 17 00:00:00 2001
From: Sondre Sortland <sonso@equinor.com>
Date: Wed, 15 Jan 2020 15:17:36 +0100
Subject: [PATCH] Fix error in triangular distribution

---
 lib/enkf/tests/trans_func.cpp | 15 +++++++++++++++
 lib/enkf/trans_func.cpp       |  7 ++++---
 2 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/lib/enkf/tests/trans_func.cpp b/lib/enkf/tests/trans_func.cpp
index 9590108a37..1594fd556f 100644
--- a/lib/enkf/tests/trans_func.cpp
+++ b/lib/enkf/tests/trans_func.cpp
@@ -40,6 +40,20 @@ void test_triangular() {
   stringlist_free(args);
 }
 
+void test_triangular_assymetric() {
+  stringlist_type * args = stringlist_alloc_new();
+  stringlist_append_copy(args , "TRIANGULAR");
+  stringlist_append_copy(args, "0");
+  stringlist_append_copy(args,"1.0");
+  stringlist_append_copy(args, "4.0");
+
+  trans_func_type * trans_func = trans_func_alloc(args);
+  test_assert_double_equal( trans_func_eval(trans_func, -1.0), 0.7966310411513150456286);
+  test_assert_double_equal( trans_func_eval(trans_func, 1.1), 2.72407181575270778882286);
+  trans_func_free( trans_func );
+  stringlist_free(args);
+}
+
 void test_create() {
   {
     stringlist_type * args = stringlist_alloc_new();
@@ -79,5 +93,6 @@ void test_create() {
 int main(int argc , char ** argv) {
   test_create();
   test_triangular();
+  test_triangular_assymetric();
 }
 
diff --git a/lib/enkf/trans_func.cpp b/lib/enkf/trans_func.cpp
index f2567d61ae..7a0a1753df 100644
--- a/lib/enkf/trans_func.cpp
+++ b/lib/enkf/trans_func.cpp
@@ -182,14 +182,15 @@ static double trans_triangular(double x, const double_vector_type * arg) {
   double xmode = double_vector_iget(arg, 1);
   double xmax = double_vector_iget(arg,2);
 
-  double inv_norm = (xmax - xmin) * (xmode- xmin);
+  double inv_norm_left = (xmax - xmin) * (xmode - xmin);
+  double inv_norm_right = (xmax - xmin) * (xmax - xmode);
   double ymode = (xmode - xmin) / (xmax - xmin);
   double y = 0.5*(1 + erf(x/sqrt(2.0)));           /* 0 - 1 */
 
   if (y < ymode)
-    return xmin + sqrt(y * inv_norm);
+    return xmin + sqrt(y * inv_norm_left);
   else
-    return xmax - sqrt((1 - y)*inv_norm);
+    return xmax - sqrt((1 - y)*inv_norm_right);
 }