I have a complicated equation with lots of pow's and was getting a NaN result. So I broke it into pieces, called temp1 to temp4. temp4
which is pow(temp3, 0.25)
is where I get a nan which then results in a NaN being returned.
#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;
const double kVal = 273.15;
// [[Rcpp::export]]
double thermalRadiance(double tas, double wind, double Tg) {
double temp1 = pow(Tg kVal, 4.0);
double temp2 = pow(wind, 0.6) * (Tg - tas);
double temp3 = (temp1 2.5e 8 * temp2);
double temp4 = pow(temp3, 0.25);
double tr = temp4 - kVal;
Rcpp::Rcout << "tas " << tas << " wind " << wind << " Tg " << Tg << " temp1 " << temp1 << " temp2 " << temp2 << " temp3 " << temp3 << " temp4 " << temp4 << " tr " << tr << std::endl;
return tr;
}
/*** R
thermalRadiance(29., 4.5, 17.4)
# test of temp4 in R
-2.37018e 07^.25
*/
tas 29 wind 4.5 Tg 17.4 temp1 7.12662e 09 temp2 -28.6013 temp3 -2.37018e 07 temp4 nan tr nan
[1] NaN
> # test of temp4 in R
> -2.37018e 07^.25
[1] -69.77427
Seems really straightforward until I get to the temp4 value. Any assistance greatly appreciated!
CodePudding user response:
You are trying to take the 4th root of a negative number in C . Perhaps you should replace pow(temp3, 0.25)
with -pow(abs(temp3), 0.25)
#include <Rcpp.h>
#include <cmath>
using namespace Rcpp;
const double kVal = 273.15;
// [[Rcpp::export]]
double thermalRadiance(double tas, double wind, double Tg) {
double temp1 = pow(Tg kVal, 4.0);
double temp2 = pow(wind, 0.6) * (Tg - tas);
double temp3 = (temp1 2.5e 8 * temp2);
double temp4 = -pow(abs(temp3), 0.25);
double tr = temp4 - kVal;
Rcpp::Rcout << "tas " << tas << " wind " << wind << " Tg " << Tg << " temp1 "
<< temp1 << " temp2 " << temp2 << " temp3 " << temp3 << " temp4 "
<< temp4 << " tr " << tr << std::endl;
return tr;
}
Which in R gives you:
thermalRadiance(29., 4.5, 17.4)
tas 29 wind 4.5 Tg 17.4 temp1 7.12662e 09 temp2 -28.6013 temp3 -2.37018e 07 temp4 -69.7743 tr -342.924
[1] -342.9243