Skybright.cpp   Skybright.cpp 
skipping to change at line 36 skipping to change at line 36
# include <math.h> # include <math.h>
# define pow10(x) std::exp((x) * 2.3025850930f) # define pow10(x) std::exp((x) * 2.3025850930f)
#endif #endif
#include "Skybright.hpp" #include "Skybright.hpp"
//! Compute exp(x) for small x //! Compute exp(x) for small x
inline float fastExp(float x) inline float fastExp(float x)
{ {
return (x>=0)? return (x>=0)?
1.f + x*(1.f+ x/2.f*(1.f+ x/3.f*(1.f+x/4.f*(1.f+x/5.f)))): (1.f + x*(1.f+ x/2.f*(1.f+ x/3.f*(1.f+x/4.f*(1.f+x/5.f))))):
1.f / fastExp(-x); 1.f / (1.f -x*(1.f -x/2.f*(1.f- x/3.f*(1.f-x/4.f*(1.
f-x/5.f)))));
} }
//! Compute acos(x) //! Compute acos(x)
//! The taylor serie is not accurate around x=1 and x=-1 //! The taylor serie is not accurate around x=1 and x=-1
inline float fastAcos(float x) inline float fastAcos(float x)
{ {
return M_PI/2.f - ( return M_PI_2 - (x + x*x*x * (1.f/6.f + x*x * (3.f/40.f + 5.f/112.f
x + 1.f/6.f * x*x*x + 3.f/40.f * x*x*x*x*x + * x*x)) );
5.f/112.f * x*x*x*x*x*x*x
//105.f/3456.f * x*x*x*x*x*x*x*x*x
);
} }
Skybright::Skybright() : SN(1.f) Skybright::Skybright() : SN(1.f)
{ {
setDate(2003, 8, 0); setDate(2003, 8, 0);
setLocation(M_PI_4, 1000., 25.f, 40.f); setLocation(M_PI_4, 1000., 25.f, 40.f);
setSunMoon(0.5, 0.5); setSunMoon(0.5, 0.5);
} }
// month : 1=Jan, 12=Dec // month : 1=Jan, 12=Dec
skipping to change at line 117 skipping to change at line 113
} }
// Compute the luminance at the given position // Compute the luminance at the given position
// Inputs : cosDistMoon = cos(angular distance between moon and the positio n) // Inputs : cosDistMoon = cos(angular distance between moon and the positio n)
// cosDistSun = cos(angular distance between sun and the position) // cosDistSun = cos(angular distance between sun and the position)
// cosDistZenith = cos(angular distance between zenith and the position) // cosDistZenith = cos(angular distance between zenith and the position)
float Skybright::getLuminance(float cosDistMoon, float Skybright::getLuminance(float cosDistMoon,
float cosDistSun, float cosDistSun,
float cosDistZenith) const float cosDistZenith) const
{ {
float distSun,dist_zenith;
// catch rounding errors here or end up with white flashes in some case
s
if (cosDistSun <= -1.f) {cosDistSun = -1.f;distSun = M_PI;}
else if (cosDistSun >= 1.f) {cosDistSun = 1.f;distSun = 0.f;}
else distSun = fastAcos(cosDistSun);
if (cosDistZenith <= -1.f) {cosDistZenith = -1.f;dist_zenith = M_PI;
}
else if (cosDistZenith >= 1.f) {cosDistZenith = 1.f;dist_zenith = 0.
f;}
else dist_zenith = fastAcos(cosDistZenith);
// Air mass // Air mass
// const float X = 1.f / (cosDistZenith + 0.025f*std::exp(-11.f*cosD const float bKX = pow10(-0.4f * K * (1.f / (cosDistZenith + 0.025f*f
istZenith)); astExp(-11.f*cosDistZenith))));
const float X = 1.f / (cosDistZenith + 0.025f*fastExp(-11.f*cosDistZ
enith));
const float bKX = pow10(-0.4f * K * X);
// Daylight brightness // Daylight brightness
const float distSun = fastAcos(cosDistSun);
const float FS = 18886.28f / (distSun*distSun + 0.0007f) const float FS = 18886.28f / (distSun*distSun + 0.0007f)
+ pow10(6.15f - (distSun+0.001f)* 1.43239f) + pow10(6.15f - (distSun+0.001f)* 1.43239f)
+ 229086.77f * ( 1.06f + cosDistSun*cosDistSun ); + 229086.77f * ( 1.06f + cosDistSun*cosDistSun );
const float b_daylight = 9.289663e-12 * (1.f - bKX) * (FS * C4 + 440 000.f * (1.f - C4)); const float b_daylight = 9.289663e-12 * (1.f - bKX) * (FS * C4 + 440 000.f * (1.f - C4));
//Twilight brightness //Twilight brightness
const float Ktrimed = K> 0.05f ? K : 0.05f; const float b_twilight = pow10(bTwilightTerm + 0.063661977f * fastAc
const float b_twilight = pow10(bTwilightTerm + 0.063661977f * dist_z os(cosDistZenith)/(K> 0.05f ? K : 0.05f)) * (1.7453293f / distSun) * (1.f-b
enith/Ktrimed) * (1.7453293f / distSun) * (1.f-bKX); KX);
// Total sky brightness // Total sky brightness
float b_total = ((b_twilight<b_daylight) ? b_twilight : b_daylight); float b_total = ((b_twilight<b_daylight) ? b_twilight : b_daylight);
// Moonlight brightness, don't compute if less than 1% daylight // Moonlight brightness, don't compute if less than 1% daylight
if ((bMoonTerm1 * (1.f - bKX) * (28860205.1341274269 * C3 + 440000.f * (1.f - C3)))/b_total>0.01f) if ((bMoonTerm1 * (1.f - bKX) * (28860205.1341274269 * C3 + 440000.f * (1.f - C3)))/b_total>0.01f)
{ {
float dist_moon; float dist_moon;
if (cosDistMoon <= -1.f) {cosDistMoon = -1.f;dist_moon = M_P if (cosDistMoon >= 1.f) {cosDistMoon = 1.f;dist_moon = 0.f;}
I;}
else if (cosDistMoon >= 1.f) {cosDistMoon = 1.f;dist_moon =
0.f;}
else else
{ {
// Because the accuracy of our power serie is bad ar ound 1, call the real acos if it's the case // Because the accuracy of our power serie is bad ar ound 1, call the real acos if it's the case
dist_moon = cosDistMoon > 0.99 ? std::acos(cosDistMo on) : fastAcos(cosDistMoon); dist_moon = cosDistMoon > 0.99 ? std::acos(cosDistMo on) : fastAcos(cosDistMoon);
} }
const float FM = 18886.28f / (dist_moon*dist_moon + 0.0005f) // The last 0.0005 should be 0, but it causes too fast brightness chang e const float FM = 18886.28f / (dist_moon*dist_moon + 0.0005f) // The last 0.0005 should be 0, but it causes too fast brightness chang e
+ pow10(6.15f - dist_moon * 1.43239f) + pow10(6.15f - dist_moon * 1.43239f)
+ 229086.77f * ( 1.06f + cosDistMoon*cosDistMoon ); + 229086.77f * ( 1.06f + cosDistMoon*cosDistMoon );
float b_moon = bMoonTerm1 * (1.f - bKX) * (FM * C3 + 440000. b_total += bMoonTerm1 * (1.f - bKX) * (FM * C3 + 440000.f *
f * (1.f - C3)); (1.f - C3));
b_total += b_moon;
} }
// Dark night sky brightness, don't compute if less than 1% daylight // Dark night sky brightness, don't compute if less than 1% daylight
if ((bNightTerm*bKX)/b_total>0.01f) if ((bNightTerm*bKX)/b_total>0.01f)
{ {
const float b_night = (0.4f + 0.6f / std::sqrt(0.04f + 0.96f b_total += (0.4f + 0.6f / std::sqrt(0.04f + 0.96f * cosDistZ
* cosDistZenith*cosDistZenith)) * bNightTerm * bKX; enith*cosDistZenith)) * bNightTerm * bKX;
b_total += b_night;
} }
return (b_total<0.f) ? 0.f : b_total * 900900.9f * M_PI * 1e-4 * 323 9389.*2. *1.5; return (b_total<0.f) ? 0.f : b_total * (900900.9f * M_PI * 1e-4 * 32 39389.*2. *1.5);
//5; // In cd/m^2 : the 32393895 is empirical term because the //5; // In cd/m^2 : the 32393895 is empirical term because the
// lambert -> cd/m^2 formula seems to be wrong... // lambert -> cd/m^2 formula seems to be wrong...
} }
/* /*
250 REM Visual limiting magnitude 250 REM Visual limiting magnitude
260 BL=B(3)/1.11E-15 : REM in nanolamberts*/ 260 BL=B(3)/1.11E-15 : REM in nanolamberts*/
// Airmass for each component // Airmass for each component
//cosDistZenith = std::cos(dist_zenith); //cosDistZenith = std::cos(dist_zenith);
//float gaz_mass = 1.f / ( cosDistZenith + 0.0286f * std::exp(-10.5f * cosDistZenith) ); //float gaz_mass = 1.f / ( cosDistZenith + 0.0286f * std::exp(-10.5f * cosDistZenith) );
 End of changes. 10 change blocks. 
40 lines changed or deleted 17 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/