Rcpp Unexpected NaN results

Question

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!

Answer 1

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