Skybright.cpp   Skybright.cpp 
skipping to change at line 39 skipping to change at line 39
{ {
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
void Skybright::setDate(const int year, const int month, const float moonPh ase) void Skybright::setDate(const int year, const int month, const float moonPh ase)
{ {
magMoon = -12.73f + 1.4896903f * std::fabs(moonPhase) + 0.04310727f * std::pow(moonPhase, 4.f); magMoon = -12.73f + 1.4896903f * fabsf(moonPhase) + 0.04310727f * po wf(moonPhase, 4.f);
// GZ: Bah, a very crude estimate for the solar position... // GZ: Bah, a very crude estimate for the solar position...
RA = (month - 3.f) * 0.52359878f; RA = (month - 3.f) * 0.52359878f;
// Term for dark sky brightness computation. // Term for dark sky brightness computation.
// GZ: This works for a few 11-year solar cycles around 1992... // GZ: This works for a few 11-year solar cycles around 1992... ...
bNightTerm = 1.0e-13 + 0.3e-13 * std::cos(0.57118f * (year-1992.f)); cos((y-1992)/11 * 2pi)
bNightTerm = 1.0e-13 + 0.3e-13 * cosf(0.57118f * (year-1992.f));
} }
void Skybright::setLocation(const float latitude, const float altitude, con st float temperature, const float relativeHumidity) void Skybright::setLocation(const float latitude, const float altitude, con st float temperature, const float relativeHumidity)
{ {
float sign_latitude = (latitude>=0.f) * 2.f - 1.f; float sign_latitude = (latitude>=0.f) * 2.f - 1.f;
// extinction Coefficient for V band // extinction Coefficient for V band
float KR = 0.1066f * std::exp(-altitude/8200.f); // Rayleigh // GZ TODO: re-create UBVRI for colored extinction, and get RGB exti
float KA = 0.1f * std::exp(-altitude/1500.f) * std::pow(1.f - 0.32f/ nction factors from SkyBright!
std::log(relativeHumidity/100.f) ,1.33f) * float KR = 0.1066f * expf(-altitude/8200.f); // Rayleigh
(1.f + 0.33f * sign_latitude * std::sin(RA)); // Aerosol float KA = 0.1f * expf(-altitude/1500.f) * powf(1.f - 0.32f/logf(rel
float KO = 0.031f * std::exp(-altitude/8200.f) * ( 3.f + 0.4f * (lat ativeHumidity/100.f) ,1.33f) *
itude * std::cos(RA) - std::cos(3.f*latitude)) )/3.f; // Ozone (1.f + 0.33f * sign_latitude * sinf(RA)); // Aerosol
float KW = 0.031f * 0.94f * (relativeHumidity/100.f) * std::exp(temp float KO = 0.031f * expf(-altitude/8200.f) * ( 3.f + 0.4f * (latitud
erature/15.f) * std::exp(-altitude/8200.f); // Water e * cosf(RA) - cosf(3.f*latitude)) )/3.f; // Ozone
float KW = 0.031f * 0.94f * (relativeHumidity/100.f) * expf(temperat
ure/15.f) * expf(-altitude/8200.f); // Water
K = KR + KA + KO + KW; // Total extinction coefficient K = KR + KA + KO + KW; // Total extinction coefficient
} }
// Set the moon and sun zenith angular distance (cosin given) // Set the moon and sun zenith angular distance (cosin given)
// and precompute what can be // and precompute what can be
void Skybright::setSunMoon(const float cosDistMoonZenith, const float cosDi stSunZenith) void Skybright::setSunMoon(const float cosDistMoonZenith, const float cosDi stSunZenith)
{ {
// Air mass for Moon // Air mass for Moon
if (cosDistMoonZenith<0) airMassMoon = 40.f; if (cosDistMoonZenith<0) airMassMoon = 40.f;
else airMassMoon = 1.f / (cosDistMoonZenith+0.025f*std::exp(-11.f*co sDistMoonZenith)); else airMassMoon = 1.f / (cosDistMoonZenith+0.025f*expf(-11.f*cosDis tMoonZenith));
// Air mass for Sun // Air mass for Sun
if (cosDistSunZenith<0) airMassSun = 40.f; if (cosDistSunZenith<0) airMassSun = 40.f;
else airMassSun = 1.f / (cosDistSunZenith+0.025f*std::exp(-11.f*cosD istSunZenith)); else airMassSun = 1.f / (cosDistSunZenith+0.025f*expf(-11.f*cosDistS unZenith));
bMoonTerm1 = stelpow10f(-0.4f * (magMoon + 54.32f)); bMoonTerm1 = stelpow10f(-0.4f * (magMoon + 54.32f));
// Moon should have no impact if below the horizon // Moon should have no impact if below the horizon
// .05 is ad hoc fadeout range - Rob // .05 is ad hoc fadeout range - Rob
if( cosDistMoonZenith < 0.f ) bMoonTerm1 *= 1.f + cosDistMoonZenith/ 0.05f; if( cosDistMoonZenith < 0.f ) bMoonTerm1 *= 1.f + cosDistMoonZenith/ 0.05f;
if(cosDistMoonZenith < -0.05f) bMoonTerm1 = 0.f; if(cosDistMoonZenith < -0.05f) bMoonTerm1 = 0.f;
C3 = stelpow10f(-0.4f*K*airMassMoon); // Term for moon brightness computation C3 = stelpow10f(-0.4f*K*airMassMoon); // Term for moon brightness computation
bTwilightTerm = -6.724f + 22.918312f * (M_PI_2-std::acos(cosDistSunZ enith)); bTwilightTerm = -6.724f + 22.918312f * (M_PI_2-acosf(cosDistSunZenit h));
C4 = stelpow10f(-0.4f*K*airMassSun); // Term for sky brightness c omputation C4 = stelpow10f(-0.4f*K*airMassSun); // Term for sky brightness c omputation
} }
// 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,
const float cosDistSun, const float cosDistSun,
skipping to change at line 125 skipping to change at line 126
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.1341274269f * C3 + 440000. f * (1.f - C3)))/b_total>0.01f) if ((bMoonTerm1 * (1.f - bKX) * (28860205.1341274269f * 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.99f ? std::acos(cosDistM oon) : StelUtils::fastAcos(cosDistMoon); dist_moon = cosDistMoon > 0.99f ? acosf(cosDistMoon) : 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
+ stelpow10f(6.15f - dist_moon * 1.43239f) + stelpow10f(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)
{ {
b_total += (0.4f + 0.6f / std::sqrt(0.04f + 0.96f * cosDistZ enith*cosDistZenith)) * bNightTerm * bKX; b_total += (0.4f + 0.6f / sqrtf(0.04f + 0.96f * cosDistZenit h*cosDistZenith)) * bNightTerm * bKX;
} }
return (b_total<0.f) ? 0.f : b_total * (900900.9f * static_cast<floa t>(M_PI) * 1e-4f * 3239389.f*2.f *1.5f); return (b_total<0.f) ? 0.f : b_total * (900900.9f * static_cast<floa t>(M_PI) * 1e-4f * 3239389.f*2.f *1.5f);
//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...
} }
 End of changes. 8 change blocks. 
16 lines changed or deleted 19 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/