Changeset 13063


Ignore:
Timestamp:
08/16/12 09:42:18 (13 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixed problem with cubic root

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/shared/Numerics/cubic.cpp

    r13061 r13063  
    11#include <math.h>
    22#include "./numerics.h"
     3
     4IssmDouble CBRT(IssmDouble Z){
     5
     6        IssmDouble ret;
     7
     8        if (Z> 0.0){
     9                ret = fabs(pow(fabs(Z),1./3.));
     10        }
     11        else if(Z< 0.0){
     12                ret = - fabs(pow(fabs(Z),1./3.));
     13        }
     14        else{
     15                ret = 0.;
     16        }
     17        return ret;
     18}
    319
    420int cubic(IssmDouble a,IssmDouble b,IssmDouble c,IssmDouble d, IssmDouble x[3], int* num){
     
    3854                        // only one real solution!
    3955                        delta = sqrt(delta);
    40                         x[0] = pow(Q+delta,third)+pow(Q-delta,third)-W;
     56                        x[0] = CBRT(Q+delta)+CBRT(Q-delta)-W;
    4157                        *num=1;
    4258                }
Note: See TracChangeset for help on using the changeset viewer.