[gradsusr] odd behavior of pow() function

[Jennifer Adams]

Interesting question! From the man page on 'pow': pow(x, y) returns a NaN and raises the "invalid" floating-point exception for finite x < 0 and finite non-integer y.

Here are some GrADS calls that illustrate this functionality of 'pow':
ga-> d pow(2,0.5) 
Result value = 1.41421 
ga-> d pow(-2,2)  
Result value = 4 
ga-> d pow(-2,-1)
Result value = -0.5 
ga-> d pow(-2,0.9999999999)
Result value = nan 
ga-> d pow(-2,1.1)         
Result value = nan 

The problem is that GrADS is returning the result of the call to 'pow' without checking whether it is a NaN and turning it into a proper missing data value. I will add this to my list of bugs to fix. Meanwhile, wrap a maskout around the expression 1000/(lev-e) and that should any negative values for the first argument to 'pow'. 
--Jennifer



On Jan 18, 2013, at 1:21 PM, Jeff Duda wrote:

> Hello,
> I am trying to compute some thermodynamic parameters, shown below.
> 
>  'define e = rhprs*6.11*exp((2.5e6/461.5)*((1/273)-(1/tmpprs)))'
>  'define r = 622*e/(lev-e)'
>  'define condtemp = 55 + 1 / (1/(tmpprs-55) - log(rhprs/100)/2840)'
>  'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'
>  'define thetae = thetaDL*exp(((3.036/condtemp)-0.00178)*r*(1+0.448+r/1000))'
>  'define thetaw = 45.114 - 51.489*pow(thetae/273.15,-3.504)'
>  'define LSI = max(thetaw,lev=1000,lev=500) - max(thetaw,lev=1000,lev=900)'
> 
> I get a lot of NaNs.  They are coming from this line:
>  'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'
> 
> And in fact, from this piece:
> pow(1000/(lev-e),0.2854)
> 
> I know this because I tested this for a model grid column (from 1000 - 500 mb). 
> ga-> d 1000/(lev-e)            
> Printing Grid -- 13 Values -- Undef = -9.99e+08
> -3.55595 -4.04332 -4.69364 -9.98891 39.3747 -5.3341 -11.5849 4.19346 
> 2.07427 2.42358 2.11214 2.1287 2.24423 
> ga-> d pow(1000/(lev-e),0.2854)
> Printing Grid -- 13 Values -- Undef = -9.99e+08
> -nan -nan -nan -nan 2.85282 -nan -nan 1.5055 
> 1.2315 1.28743 1.23787 1.24063 1.25949 
> 
> For whatever reason, the pow() function doesn't like negative bases, because for example, I can use a calculator to compute -3.55595^0.2854 as -1.43629.  Why, then, does the pow() function give a -nan as a result?
> 
> Jeff Duda
> 
> -- 
> Jeff Duda
> Graduate research assistant
> University of Oklahoma School of Meteorology
> Center for Analysis and Prediction of Storms
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr

--
Jennifer M. Adams
IGES/COLA
4041 Powder Mill Road, Suite 302
Calverton, MD 20705
jma at cola.iges.org



-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20130118/e4e4c383/attachment-0003.html