Eric and others,<div>Thanks for the insights. I totally neglected the issue that there could be problems taking non-integer powers of negative numbers. However, more importantly, I discovered there was an error in my calculation that I fixed. Since then, the problem has been fixed for the most part.</div>
<div><br></div><div>Jeff<br><br><div class="gmail_quote">On Fri, Jan 18, 2013 at 7:26 PM, Eric Altshuler <span dir="ltr"><<a href="mailto:ela@cola.iges.org" target="_blank">ela@cola.iges.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><div style="font-size:12pt;font-family:arial,helvetica,sans-serif">
There seem to be two separate issues here. One is the pow function returning 'nan' instead of the missing data value. As Jennifer says, this is a bug. However, the other issue is that lev-e should never be negative. The vapor pressure, e, should be much smaller than the pressure, lev (since lev is in millibars, e must also be in mb). In the calculation of e, I suspect 'rhprs' is given as a percentage, in which case e will be a factor of 100 too large. Multiply by 0.0611 instead of 6.11 to fix this problem. (Since I'm not familiar with your other formulas, I can't vouch for their correctness.)</div>
<div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><br></div><div style="font-size:12pt;font-family:arial,helvetica,sans-serif">Using maskout in this situation will get rid of the nan's but it won't fix the underlying problem. You'll get missing values where lev-e is (erroneously) negative, and non-missing but spurious values where lev-e is positive (but still wrong).</div>
<div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><br></div><div style="font-size:12pt;font-family:arial,helvetica,sans-serif">Finally, your calculator's order of operations may perform the ^ before the negation: -(3.55595^0.2854) = -1.43629. No calculator should give a real result for (-3.55595)^0.2854; this is mathematically impossible. Only if the power is the exact reciprocal of an odd integer will the raising of a negative number to a (decimal) non-integer power result in a real value. As far as I can tell, this is only true for reciprocals of powers of 5 (1/5, 1/25, 1/125, etc.), i.e. 0.2, 0.04, 0.008, etc. Note that the power 1/3 (cube root) of a negative number is real, but the decimal 0.333333333... doesn't represent 1/3 exactly, so some calculators will give an error. (On calculators that accept fractional notation, the operation (-8)^(1/3) does give -2, as expected, when 1/3 is given as a fraction instead of a decimal. This is also true for other odd roots, e.g. (-128)^(1/7)=-2.)</div>
<div style="font-size:12pt;font-family:arial,helvetica,sans-serif"><br></div><div style="font-size:12pt;font-family:arial,helvetica,sans-serif">Best regards,</div><div style="font-size:12pt;font-family:arial,helvetica,sans-serif">
<br><br><div><span></span>Eric L. Altshuler<br>Research Scientist<br>Center for Ocean-Land-Atmosphere Studies<div class="im"><br>4041 Powder Mill Road, Suite 302<br></div>Calverton, MD 20705-3106<br>USA<br><br>E-mail: <a href="mailto:ela@cola.iges.org" target="_blank">ela@cola.iges.org</a><br>
Phone: (301) 902-1257<br>Fax: (301) 595-9793<span></span><br></div><br><hr><b>From: </b>"Jennifer Adams" <<a href="mailto:jma@cola.iges.org" target="_blank">jma@cola.iges.org</a>><br><b>To: </b>"GrADS Users Forum" <<a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a>><br>
<b>Sent: </b>Friday, January 18, 2013 2:09:24 PM<br><b>Subject: </b>Re: [gradsusr] odd behavior of pow() function<div><div class="h5"><br><br>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.<div>
<div><br></div><div>Here are some GrADS calls that illustrate this functionality of 'pow':</div><div><div>ga-> d pow(2,0.5) </div><div>Result value = 1.41421 </div><div>ga-> d pow(-2,2) </div><div>Result value = 4 </div>
<div>ga-> d pow(-2,-1)</div><div>Result value = -0.5 </div><div>ga-> d pow(-2,0.9999999999)</div><div>Result value = nan </div><div>ga-> d pow(-2,1.1) </div><div>Result value = nan </div><div><br></div></div>
<div>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'. </div>
<div>--Jennifer</div><div><br></div><div><br></div><div><br><div><div>On Jan 18, 2013, at 1:21 PM, Jeff Duda wrote:</div><br><blockquote>Hello,<br>I am trying to compute some thermodynamic parameters, shown below.<br clear="all">
<br> 'define e = rhprs*6.11*exp((2.5e6/461.5)*((1/273)-(1/tmpprs)))'<br> 'define r = 622*e/(lev-e)'<br> 'define condtemp = 55 + 1 / (1/(tmpprs-55) - log(rhprs/100)/2840)'<br>
'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'<br> 'define thetae = thetaDL*exp(((3.036/condtemp)-0.00178)*r*(1+0.448+r/1000))'<br> 'define thetaw = 45.114 - 51.489*pow(thetae/273.15,-3.504)'<br>
'define LSI = max(thetaw,lev=1000,lev=500) - max(thetaw,lev=1000,lev=900)'<br><br>I get a lot of NaNs. They are coming from this line:<br> 'define thetaDL = tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'<br>
<br>And in fact, from this piece:<br>pow(1000/(lev-e),0.2854)<br><br>I know this because I tested this for a model grid column (from 1000 - 500 mb). <br>ga-> d 1000/(lev-e) <br>Printing Grid -- 13 Values -- Undef = -9.99e+08<br>
-3.55595 -4.04332 -4.69364 -9.98891 39.3747 -5.3341 -11.5849 4.19346 <br>2.07427 2.42358 2.11214 2.1287 2.24423 <br>ga-> d pow(1000/(lev-e),0.2854)<br>Printing Grid -- 13 Values -- Undef = -9.99e+08<br>-nan -nan -nan -nan 2.85282 -nan -nan 1.5055 <br>
1.2315 1.28743 1.23787 1.24063 1.25949 <br><br>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?<br>
<br>Jeff Duda<br><br>-- <br>Jeff Duda<br>Graduate research assistant<br>University of Oklahoma School of Meteorology<br>Center for Analysis and Prediction of Storms<br>
_______________________________________________<br>gradsusr mailing list<br><a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br><a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
</blockquote></div><br><div>
<span style="font-size:12px"><div>--</div><div>Jennifer M. Adams</div><div>IGES/COLA</div><div>4041 Powder Mill Road, Suite 302</div><div>Calverton, MD 20705</div><div><a href="mailto:jma@cola.iges.org" target="_blank">jma@cola.iges.org</a></div>
<div><br></div><br></span>
</div>
<br></div></div><br>_______________________________________________<br>gradsusr mailing list<br><a href="mailto:gradsusr@gradsusr.org" target="_blank">gradsusr@gradsusr.org</a><br><a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
</div></div></div></div></div></div></div><br>_______________________________________________<br>
gradsusr mailing list<br>
<a href="mailto:gradsusr@gradsusr.org">gradsusr@gradsusr.org</a><br>
<a href="http://gradsusr.org/mailman/listinfo/gradsusr" target="_blank">http://gradsusr.org/mailman/listinfo/gradsusr</a><br>
<br></blockquote></div><br><br clear="all"><div><br></div>-- <br>Jeff Duda<br>Graduate research assistant<br>University of Oklahoma School of Meteorology<br>Center for Analysis and Prediction of Storms<br>
</div>