Skybright.cpp   Skybright.cpp 
skipping to change at line 31 skipping to change at line 31
#include "StelUtils.hpp" #include "StelUtils.hpp"
#include <config.h> #include <config.h>
#ifndef HAVE_POW10 #ifndef HAVE_POW10
# define HAVE_POW10 1 # define HAVE_POW10 1
# 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"
#include "StelUtils.hpp"
//! Compute exp(x) for small x
inline float fastExp(float x)
{
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 / (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)
//! The taylor serie is not accurate around x=1 and x=-1
inline float fastAcos(float x)
{
return M_PI_2 - (x + x*x*x * (1.f/6.f + x*x * (3.f/40.f + 5.f/112.f
* 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
// moonPhase in radian 0=Full Moon, PI/2=First Quadrant/Last Quadran, PI=No Moon // moonPhase in radian 0=Full Moon, PI/2=First Quadrant/Last Quadran, PI=No Moon
skipping to change at line 114 skipping to change at line 100
// 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
{ {
// Air mass // Air mass
const float bKX = pow10(-0.4f * K * (1.f / (cosDistZenith + 0.025f*f astExp(-11.f*cosDistZenith)))); const float bKX = pow10(-0.4f * K * (1.f / (cosDistZenith + 0.025f*S telUtils::fastExp(-11.f*cosDistZenith))));
// Daylight brightness // Daylight brightness
const float distSun = fastAcos(cosDistSun); const float distSun = StelUtils::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 b_twilight = pow10(bTwilightTerm + 0.063661977f * fastAc os(cosDistZenith)/(K> 0.05f ? K : 0.05f)) * (1.7453293f / distSun) * (1.f-b KX); const float b_twilight = pow10(bTwilightTerm + 0.063661977f * StelUt ils::fastAcos(cosDistZenith)/(K> 0.05f ? K : 0.05f)) * (1.7453293f / distSu n) * (1.f-bKX);
// 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 = 0.f;} 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) : StelUtils::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 );
b_total += bMoonTerm1 * (1.f - bKX) * (FM * C3 + 440000.f * (1.f - C3)); b_total += bMoonTerm1 * (1.f - bKX) * (FM * C3 + 440000.f * (1.f - C3));
} }
// 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)
 End of changes. 5 change blocks. 
21 lines changed or deleted 5 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/