WSL/SLF GitLab Repository

Commit 9293031a authored by Mathias Bavay's avatar Mathias Bavay
Browse files

New optimized pow() replacement have been implemented. The Sun object can not...

New optimized pow() replacement have been implemented. The Sun object can not use them, though (error accumulation leading to a dramatic loss of accuracy). The Sun object now checks that a given julian_date is different from the already stored one (to avoid recalculating all if not necessary).
parent 8eba7de6
......@@ -114,6 +114,54 @@ namespace Optim {
inline double pow2(const double& val) {return (val*val);}
inline double pow3(const double& val) {return (val*val*val);}
inline double pow4(const double& val) {return (val*val*val*val);}
//please do not use this method directly, call fastPow() instead!
inline double fastPowInternal(double a, double b) {
//see http://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
// calculate approximation with fraction of the exponent
int e = (int) b;
union {
double d;
int x[2];
} u = { a };
u.x[1] = (int)((b - e) * (u.x[1] - 1072632447) + 1072632447);
u.x[0] = 0;
// exponentiation by squaring with the exponent's integer part
// double r = u.d makes everything much slower, not sure why
double r = 1.0;
while (e) {
if (e & 1) {
r *= a;
}
a *= a;
e >>= 1;
}
return r * u.d;
}
/**
* @brief Optimized version of c++ pow()
* This version works with positive and negative exponents and handles exponents bigger than 1.
* The relative error remains less than 6% for the benachmarks that we ran (argument between 0 and 500
* and exponent between -10 and +10). It is ~3.3 times faster than cmath's pow().
* Source: http://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp/
*
* Please benchmark your code before deciding to use this!!
* @param a argument
* @param b exponent
* @return a^b
*/
inline double fastPow(double a, double b) {
if(b>0.) {
return fastPowInternal(a,b);
} else {
const double tmp = fastPowInternal(a,-b);
return 1./tmp;
}
}
}
} //end namespace
......
......@@ -54,11 +54,18 @@ SunObject::SunObject(const double& i_latitude, const double& i_longitude, const
}
void SunObject::setDate(const double& i_julian, const double& TZ) {
position.reset();
julian_gmt = i_julian - TZ/24.;
beam_toa = beam_direct = beam_diffuse = IOUtils::nodata;
if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata && altitude!=IOUtils::nodata) {
update();
const double i_julian_gmt = i_julian - TZ/24.;
if(i_julian_gmt!=julian_gmt) { //invalidate all fields if receiving a new date
position.reset();
julian_gmt = i_julian_gmt;
beam_toa = beam_direct = beam_diffuse = IOUtils::nodata;
}
//if the date was new or if the previous date had not lead to an update -> update now
if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata && altitude!=IOUtils::nodata &&
(beam_toa==IOUtils::nodata || beam_direct==IOUtils::nodata || beam_diffuse==IOUtils::nodata)) {
update();
}
}
......@@ -110,7 +117,7 @@ void SunObject::getBeamPotential(const double& sun_elevation, const double& Ecce
const double& ta, const double& rh, const double& pressure, const double& mean_albedo,
double& R_toa, double& R_direct, double& R_diffuse)
{
//TODO: find a suitable replacement for pow(): these pow cost us a lot here
//these pow cost us a lot here, but replacing them by fastPow() has a large impact o accuracy (because of the exp())
const double olt = 0.32; // ozone layer thickness (cm) U.S.standard = 0.34 cm
const double w0 = 0.9; // fraction of energy scattered to total attenuation by aerosols (Bird and Hulstrom(1981))
const double fc = 0.84; // fraction of forward scattering to total scattering (Bird and Hulstrom(1981))
......@@ -148,12 +155,12 @@ void SunObject::getBeamPotential(const double& sun_elevation, const double& Ecce
// on the quoted pages
// broadband transmittance by Rayleigh scattering (Iqbal (1983), p.189)
const double taur = exp( -0.0903 * pow( ma,0.84 ) * (1 + ma - pow( ma,1.01 )) );
const double taur = exp( -0.0903 * pow( ma,0.84 ) * (1. + ma - pow( ma,1.01 )) );
// broadband transmittance by ozone (Iqbal (1983), p.189)
const double u3 = olt * mr; // ozone relative optical path length
const double tauoz = 1. - (0.1611 * u3 * pow( 1. + 139.48 * u3,-0.3035 ) -
0.002715 * u3 * pow( 1. + 0.044 * u3 + 0.0003 * u3 * u3,-1 ));
0.002715 * u3 * 1./( 1. + 0.044 * u3 + 0.0003 * u3 * u3));
// broadband transmittance by uniformly mixed gases (Iqbal (1983), p.189)
const double taug = exp( -0.0127 * pow( ma,0.26 ) );
......@@ -169,7 +176,7 @@ void SunObject::getBeamPotential(const double& sun_elevation, const double& Ecce
const double u1 = precw * mr;
// broadband transmittance by water vapor (in Iqbal (1983), p.189)
const double tauw = 1 - 2.4959 * u1 * (1. / (pow( 1.0 + 79.034 * u1,0.6828 ) + 6.385 * u1));
const double tauw = 1. - 2.4959 * u1 * (1. / (pow( 1.0 + 79.034 * u1,0.6828 ) + 6.385 * u1));
// broadband total transmittance by aerosols (in Iqbal (1983), pp.189-190)
// using Angstroem's turbidity formula Angstroem (1929, 1930) for the aerosol thickness
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment