Changeset 13063

Timestamp:

08/16/12 09:42:18 (13 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed problem with cubic root

File:

• issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp

@@ -1 +1 @@
 #include <math.h>
 #include "./numerics.h"
+
+IssmDouble CBRT(IssmDouble Z){
+
+        IssmDouble ret;
+
+        if (Z> 0.0){
+                ret = fabs(pow(fabs(Z),1./3.));
+        }
+        else if(Z< 0.0){
+                ret = - fabs(pow(fabs(Z),1./3.));
+        }
+        else{
+                ret = 0.;
+        }
+        return ret;
+}
 
 int cubic(IssmDouble a,IssmDouble b,IssmDouble c,IssmDouble d, IssmDouble x[3], int* num){
@@ -38 +54 @@
                         // only one real solution!
                         delta = sqrt(delta);
-                        x[0] = pow(Q+delta,third)+pow(Q-delta,third)-W;
+                        x[0] = CBRT(Q+delta)+CBRT(Q-delta)-W;
                         *num=1;
                 }