From 435110cef9b8b1d4ee4e71b8595d1a162af5a956 Mon Sep 17 00:00:00 2001 From: Hirek193 Date: Sat, 17 Oct 2026 22:47:24 +0200 Subject: [PATCH] Place the sun and the moon where they are and light the night by the phase of the moon --- docs/opengl46-renderer.md | 2 + environment/moon.cpp | 128 ++++++++++++++---- environment/moon.h | 15 ++- environment/sun.cpp | 35 ++++- rendering/opengl46/opengl46renderer.cpp | 171 +++++++++++++----------- shaders/celestial.frag | 64 ++++++++- shaders/celestial.vert | 4 + simulation/simulationenvironment.cpp | 24 +++- utilities/Globals.cpp | 4 + utilities/Globals.h | 4 + 10 files changed, 338 insertions(+), 113 deletions(-) diff --git a/docs/opengl46-renderer.md b/docs/opengl46-renderer.md index 21f72faf..92f25434 100644 --- a/docs/opengl46-renderer.md +++ b/docs/opengl46-renderer.md @@ -255,6 +255,8 @@ are no longer loaded. | `gfx.clouds.shadows.range` | `16000` | 1000..100000: the side of the cloud shadow map, metres; farther the clouds cast no shadows | | `gfx.stars yes/no` | `yes` | the stars of the catalogue | | `gfx.stars.brightness` | `1.0` | 0..16: their brightness | +| `gfx.sun.disc` | `1.0` | 0..4: brightness of the sun disc against the radiance its sunlight gives it; 1 is physical (held within the range of the image), the glare around it comes from the bloom | +| `gfx.moon.disc` | `1.0` | 0..16: brightness of the moon disc | ### Models drawn from the GPU diff --git a/environment/moon.cpp b/environment/moon.cpp index 302e5d8e..9c81304b 100644 --- a/environment/moon.cpp +++ b/environment/moon.cpp @@ -6,7 +6,30 @@ #include "simulation/simulationtime.h" ////////////////////////////////////////////////////////////////////////////////////////// -// cSun -- class responsible for dynamic calculation of position and intensity of the Sun, +// cMoon -- class responsible for dynamic calculation of position, phase and intensity of the Moon + +namespace { + +// direction to a body on the sky of the observer (x to the west, y up, z to the north) from its hour angle and declination, +// at the given (refracted) elevation. degrees +glm::vec3 +horizontal_direction( double const Hourangle, double const Declination, double const Latitude, double const Elevation ) { + + auto const hourangle { glm::radians( Hourangle ) }; + auto const declination { glm::radians( Declination ) }; + auto const latitude { glm::radians( Latitude ) }; + auto const elevation { glm::radians( Elevation ) }; + auto const west { std::cos( declination ) * std::sin( hourangle ) }; + auto const north { std::cos( latitude ) * std::sin( declination ) - std::sin( latitude ) * std::cos( declination ) * std::cos( hourangle ) }; + auto const horizontal { std::hypot( west, north ) }; + if( horizontal < 1e-9 ) { + return { 0.f, ( Elevation >= 0.0 ? 1.f : -1.f ), 0.f }; + } + auto const scale { std::cos( elevation ) / horizontal }; + return glm::normalize( glm::vec3( static_cast( west * scale ), static_cast( std::sin( elevation ) ), static_cast( north * scale ) ) ); +} + +} // namespace cMoon::cMoon() { @@ -21,23 +44,20 @@ void cMoon::init() { m_observer.timezone = -1.0 * simulation::Time.zone_bias(); - // NOTE: we're calculating phase just once, because it's unlikely simulation will last a few days, - // plus a sudden texture change would be pretty jarring + move(); phase(); } void -cMoon::update( bool const Includephase ) { +cMoon::update( [[maybe_unused]] bool const Includephase ) { m_observer.temp = Global.AirTemperature; move(); - if( Includephase ) { phase(); } - glm::vec3 position( 0.f, 0.f, -1.f ); - position = glm::rotateX( position, glm::radians( static_cast( m_body.elevref ) ) ); - position = glm::rotateY( position, glm::radians( static_cast( -m_body.hrang ) ) ); - - m_position = glm::normalize( position ); + // the phase follows the positions of the moon and the sun, worked out in move(); the moon is drawn lit from the side + // of the sun, so it changes smoothly and is kept up to date on every update + phase(); + m_position = horizontal_direction( m_body.hrang, m_body.declin, m_observer.latitude, m_body.elevref ); } void @@ -65,8 +85,31 @@ float cMoon::getIntensity() { // calculating intensity of the sun instead, and returning 15% of the value, // which roughly matches how much sunlight is reflected by the moon // We alter the intensity further based on current phase of the moon - auto const phasefactor = 1.0f - std::abs( m_phase - 29.53f * 0.5f ) / ( 29.53f * 0.5f ); - return static_cast( m_body.etr / 1399.0 * phasefactor * 0.15 ); // arbitrary scaling factor taken from etrn value + return static_cast( m_body.etr / 1399.0 * getBrightness() * 0.15 ); // arbitrary scaling factor taken from etrn value +} + +// brightness of the moon relative to the full moon. it falls with the phase much faster than the lit part of the disc +// shrinks, the quarter moon gives about a tenth of the light of the full one (phase law of the lunar magnitude, Allen 1976) +float +cMoon::getBrightness() const { + + auto const angle { std::abs( m_phaseangle ) }; + return std::pow( 10.f, -0.4f * ( 0.026f * angle + 4.0e-9f * angle * angle * angle * angle ) ); +} + +// lit fraction of the disc, 0 at the new moon, 1 at the full moon +float +cMoon::getIllumination() const { + + return 0.5f * ( 1.f + std::cos( glm::radians( m_phaseangle ) ) ); +} + +// angular radius of the disc seen from the earth, radians +float +cMoon::getAngularRadius() const { + + // radius of the moon in earth radii, over the distance in earth radii + return static_cast( std::asin( 0.2725 / std::max( m_body.distance, 50.0 ) ) ); } void cMoon::setLocation( float const Longitude, float const Latitude ) { @@ -118,7 +161,7 @@ void cMoon::move() { + 275 * localtime.wMonth / 9 + localtime.wDay - 730530 - + localut / 24.0; + + ( localut - m_observer.timezone ) / 24.0; // the orbital elements count days of universal time // Universal Coordinated (Greenwich standard) time m_observer.utime = localut - m_observer.timezone; @@ -175,19 +218,57 @@ void cMoon::move() { double const pertevection = -1.274 * std::sin( degtorad * ( m_body.mnanom - 2.0 * mnelong ) ); // Evection double const pertvariation = +0.658 * std::sin( degtorad * ( 2.0 * mnelong ) ); // Variation double const pertyearlyeqt = -0.186 * std::sin( degtorad * sunmnanom ); // Yearly equation + // smaller longitude perturbations, all of them together up to about a quarter of a degree + // Moon's argument of latitude: F = Lm - N + double const latarg = clamp_circular( m_body.mnlong - longascnode ); + auto const sindeg = [&]( double const Angle ) { return std::sin( degtorad * Angle ); }; + double const pertsmall = + - 0.059 * sindeg( 2.0 * m_body.mnanom - 2.0 * mnelong ) + - 0.057 * sindeg( m_body.mnanom - 2.0 * mnelong + sunmnanom ) + + 0.053 * sindeg( m_body.mnanom + 2.0 * mnelong ) + + 0.046 * sindeg( 2.0 * mnelong - sunmnanom ) + + 0.041 * sindeg( m_body.mnanom - sunmnanom ) + - 0.035 * sindeg( mnelong ) // Parallactic equation + - 0.031 * sindeg( m_body.mnanom + sunmnanom ) + - 0.015 * sindeg( 2.0 * latarg - 2.0 * mnelong ) + + 0.011 * sindeg( m_body.mnanom - 4.0 * mnelong ); // latitude perturbations + double const pertlatitude = + - 0.173 * sindeg( latarg - 2.0 * mnelong ) + - 0.055 * sindeg( m_body.mnanom - latarg - 2.0 * mnelong ) + - 0.046 * sindeg( m_body.mnanom + latarg - 2.0 * mnelong ) + + 0.033 * sindeg( latarg + 2.0 * mnelong ) + + 0.017 * sindeg( 2.0 * m_body.mnanom + latarg ); + // distance perturbations, earth radii + m_body.distance += + - 0.58 * std::cos( degtorad * ( m_body.mnanom - 2.0 * mnelong ) ) + - 0.46 * std::cos( degtorad * ( 2.0 * mnelong ) ); - m_body.eclong += pertevection + pertvariation + pertyearlyeqt; + m_body.eclong = clamp_circular( m_body.eclong + pertevection + pertvariation + pertyearlyeqt + pertsmall ); + // ecliptic latitude: the orbit is tilted by about 5 degrees, which moves the moon as much off the ecliptic + m_body.eclat = radtodeg * std::atan2( zeclip, std::sqrt( xeclip * xeclip + yeclip * yeclip ) ) + pertlatitude; + + // equatorial coordinates + double const xequat = std::cos( degtorad * m_body.eclong ) * std::cos( degtorad * m_body.eclat ); + double const yeclon = std::sin( degtorad * m_body.eclong ) * std::cos( degtorad * m_body.eclat ); + double const zeclat = std::sin( degtorad * m_body.eclat ); + double const yequat = yeclon * std::cos( degtorad * m_body.oblecl ) - zeclat * std::sin( degtorad * m_body.oblecl ); + double const zequat = yeclon * std::sin( degtorad * m_body.oblecl ) + zeclat * std::cos( degtorad * m_body.oblecl ); // declination - m_body.declin = radtodeg * std::asin( std::sin (m_body.oblecl * degtorad) * std::sin(m_body.eclong * degtorad) ); - + m_body.declin = radtodeg * std::atan2( zequat, std::sqrt( xequat * xequat + yequat * yequat ) ); // right ascension - double top = std::cos( degtorad * m_body.oblecl ) * std::sin( degtorad * m_body.eclong ); - double bottom = std::cos( degtorad * m_body.eclong ); - m_body.rascen = clamp_circular( radtodeg * std::atan2( top, bottom ) ); + m_body.rascen = clamp_circular( radtodeg * std::atan2( yequat, xequat ) ); + + // phase: the angle between the sun and the moon seen from the earth (elongation), from the true longitude of the sun + double const sunlong = sunmnlong + 1.915 * sindeg( sunmnanom ) + 0.020 * sindeg( 2.0 * sunmnanom ); + m_elongation = clamp_circular( m_body.eclong - sunlong ); + double const coselongation = std::cos( degtorad * m_body.eclat ) * std::cos( degtorad * m_elongation ); + // the phase angle (sun - moon - earth) is the supplement of the elongation, closely enough at the distance of the sun + m_phaseangle = static_cast( 180.0 - radtodeg * std::acos( std::clamp( coselongation, -1.0, 1.0 ) ) ); // Greenwich mean sidereal time - m_observer.gmst = 6.697375 + 0.0657098242 * daynumber + m_observer.utime; + // (Schlyter: GMST0 = sun's mean longitude + 180 degrees, for the day number counted from 2000 Jan 0.0 UT) + m_observer.gmst = 6.5991604 + 0.0657098242 * daynumber + m_observer.utime; m_observer.gmst -= 24.0 * (int)( m_observer.gmst / 24.0 ); if( m_observer.gmst < 0.0 ) m_observer.gmst += 24.0; @@ -219,6 +300,8 @@ void cMoon::move() { m_body.zenetr = std::acos( cz ) * radtodeg; m_body.elevetr = 90.0 - m_body.zenetr; + // topocentric parallax: the moon is close enough to be seen about a degree lower than from the centre of the earth + m_body.elevetr -= radtodeg * std::asin( 1.0 / m_body.distance ) * std::cos( degtorad * m_body.elevetr ); refract(); } @@ -290,9 +373,8 @@ cMoon::phase() { if (lt.wMonth == 5 && lt.wDay == 4) //May the forth be with you! m_phase = 50; else { - // calculate moon's age in days from new moon - float ip = normalize( ( static_cast(simulation::Time.julian_day()) - 2451550.1f ) / 29.530588853f ); - m_phase = ip * 29.53f; + // moon's age in days from new moon, from the elongation worked out in move() + m_phase = normalize( static_cast( m_elongation / 360.0 ) ) * 29.53f; } } diff --git a/environment/moon.h b/environment/moon.h index ca8ff29c..f51cb926 100644 --- a/environment/moon.h +++ b/environment/moon.h @@ -19,8 +19,16 @@ public: float getAngle() const; // returns current intensity of the sun float getIntensity(); - // returns current phase of the moon + // returns current phase of the moon: age in days from the new moon (50 on May 4th) float getPhase() const { return m_phase; } + // returns current phase angle (sun - moon - earth), degrees: 0 at the full moon, 180 at the new moon + float getPhaseAngle() const { return m_phaseangle; } + // returns lit fraction of the disc, 0...1 + float getIllumination() const; + // returns brightness of the moon relative to the full moon, 0...1 + float getBrightness() const; + // returns angular radius of the disc, radians + float getAngularRadius() const; // sets current time, overriding one acquired from the system clock void setTime( int const Hour, int const Minute, int const Second ); // sets current geographic location @@ -62,6 +70,7 @@ protected: double mnanom; // mean anomaly, degrees double tranom; // true anomaly, degrees double eclong; // ecliptic longitude, degrees. + double eclat{ 0.0 }; // ecliptic latitude, degrees. double oblecl; // obliquity of ecliptic. double declin; // declination--zenith angle of solar noon at equator, degrees NORTH. double rascen; // right ascension, degrees @@ -95,5 +104,7 @@ protected: celestialbody m_body; observer m_observer; glm::vec3 m_position; - float m_phase; + float m_phase{ 0.f }; + double m_elongation{ 0.0 }; // ecliptic longitude of the moon less that of the sun, degrees 0...360 + float m_phaseangle{ 0.f }; }; diff --git a/environment/sun.cpp b/environment/sun.cpp index c85fcacb..1e0ebe80 100644 --- a/environment/sun.cpp +++ b/environment/sun.cpp @@ -8,6 +8,29 @@ ////////////////////////////////////////////////////////////////////////////////////////// // cSun -- class responsible for dynamic calculation of position and intensity of the Sun, +namespace { + +// direction to a body on the sky of the observer (x to the west, y up, z to the north) from its hour angle and declination, +// at the given (refracted) elevation. degrees +glm::vec3 +horizontal_direction( double const Hourangle, double const Declination, double const Latitude, double const Elevation ) { + + auto const hourangle { glm::radians( Hourangle ) }; + auto const declination { glm::radians( Declination ) }; + auto const latitude { glm::radians( Latitude ) }; + auto const elevation { glm::radians( Elevation ) }; + auto const west { std::cos( declination ) * std::sin( hourangle ) }; + auto const north { std::cos( latitude ) * std::sin( declination ) - std::sin( latitude ) * std::cos( declination ) * std::cos( hourangle ) }; + auto const horizontal { std::hypot( west, north ) }; + if( horizontal < 1e-9 ) { + return { 0.f, ( Elevation >= 0.0 ? 1.f : -1.f ), 0.f }; + } + auto const scale { std::cos( elevation ) / horizontal }; + return glm::normalize( glm::vec3( static_cast( west * scale ), static_cast( std::sin( elevation ) ), static_cast( north * scale ) ) ); +} + +} // namespace + cSun::cSun() { setLocation( 19.00f, 52.00f ); // default location roughly in centre of Poland @@ -29,11 +52,8 @@ cSun::update() { m_observer.temp = Global.AirTemperature; move(); - glm::vec3 position( 0.f, 0.f, -1.f ); - position = glm::rotateX( position, glm::radians( static_cast( m_body.elevref ) ) ); - position = glm::rotateY( position, glm::radians( static_cast( -m_body.hrang ) ) ); - - m_position = glm::normalize( position ); + // the azimuth follows from the hour angle, the declination and the latitude; it is not the hour angle itself + m_position = horizontal_direction( m_body.hrang, m_body.declin, m_observer.latitude, m_body.elevref ); } void @@ -122,7 +142,7 @@ void cSun::move() { + 275 * localtime.wMonth / 9 + localtime.wDay - 730530 - + localut / 24.0; + + ( localut - m_observer.timezone ) / 24.0; // the orbital elements count days of universal time // Universal Coordinated (Greenwich standard) time m_observer.utime = localut - m_observer.timezone; @@ -171,7 +191,8 @@ void cSun::move() { m_body.rascen = clamp_circular( radtodeg * std::atan2( top, bottom ) ); // Greenwich mean sidereal time (hours) - m_observer.gmst = clamp_circular( 6.697375 + 0.0657098242 * daynumber + m_observer.utime, 24.0 ); + // (Schlyter: GMST0 = sun's mean longitude + 180 degrees, for the day number counted from 2000 Jan 0.0 UT) + m_observer.gmst = clamp_circular( 6.5991604 + 0.0657098242 * daynumber + m_observer.utime, 24.0 ); // local mean sidereal time (deg) m_observer.lmst = clamp_circular( m_observer.gmst * 15.0 + m_observer.longitude ); diff --git a/rendering/opengl46/opengl46renderer.cpp b/rendering/opengl46/opengl46renderer.cpp index 8d65d6bb..ad806c6f 100644 --- a/rendering/opengl46/opengl46renderer.cpp +++ b/rendering/opengl46/opengl46renderer.cpp @@ -886,6 +886,8 @@ void opengl46_renderer::draw_debug_ui() if (ImGui::Begin("Sky and exposure")) { m_sky.draw_debug_ui(); + ImGui::SliderFloat("sun disc, 1 = physical (gfx.sun.disc)", &Global.gfx_sun_disc, 0.f, 2.f); + ImGui::SliderFloat("moon disc (gfx.moon.disc)", &Global.gfx_moon_disc, 0.f, 8.f); ImGui::Separator(); m_exposure.draw_debug_ui(); ImGui::Separator(); @@ -2229,6 +2231,8 @@ bool opengl46_renderer::Render(world_environment *Environment) // celestial bodies m_celestial_shader->bind(); m_empty_vao->bind(); + // the bodies give their colour premultiplied, and as alpha how much of the sky behind them they hide (celestial.frag) + ::glBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA); // the fog and the clouds laid over the sky hide the sun and the moon; without the clouds the overcast fades them auto const overcastfactor{m_weather.clouds_drawn() ? 1.f : std::clamp(1.5f - Global.Overcast, 0.f, 1.f)}; @@ -2236,97 +2240,111 @@ bool opengl46_renderer::Render(world_environment *Environment) glm::vec3 suncolor = glm::mix(glm::vec3(1.0f, 242.0f / 255.0f, 231.0f / 255.0f), glm::vec3(235.0f / 255.0f, 140.0f / 255.0f, 36.0f / 255.0f), duskfactor); suncolor = m_sky.sun_disc_color(suncolor); + // the vertical of the sky on the quads of the bodies, which lie in the plane of the view + glm::vec2 zenith{glm::mat3{modelview} * glm::vec3{0.f, 1.f, 0.f}}; + zenith = (glm::length(zenith) > 1e-3f ? glm::normalize(zenith) : glm::vec2{0.f, 1.f}); + // a disc low over the horizon is flattened by the refraction: its lower edge is lifted more than the upper one (refraction + // after Bennett, for the apparent elevation). height over width of the disc + auto const flattening = [](float const Elevation, float const Radius) { + auto const refraction = [](float const Height) { + auto const height{std::max(Height, -2.f)}; + return 1.f / std::tan(glm::radians(height + 7.31f / (height + 4.4f))); // arc minutes + }; + auto const radius{glm::degrees(Radius)}; + return std::clamp(1.f - (refraction(Elevation - radius) - refraction(Elevation + radius)) / (2.f * radius * 60.f), 0.5f, 1.f); + }; + // the view of the frame gets the sun as bright as it is, which gives the glare of the bloom; the other views (the + // environment map, mirrors) keep the old dim disc, the reflections of the sun on surfaces come from their shading + auto const mainview{(m_renderpass.draw_mode == rendermode::color) && (m_current_viewport != nullptr) && m_current_viewport->main}; + std::array body{}; + // sun { - Bind_Texture(0, m_suntexture); - glm::vec4 color(suncolor.x, suncolor.y, suncolor.z, overcastfactor); - auto const sunvector = Environment->m_sun.getDirection(); + // the true size of the disc, 0.533 degrees across, the quad a little larger for a smooth edge + float const sunradius{glm::radians(0.2666f)}; + float const quadmargin{1.25f}; + // limb darkening per channel (linear law): the edge of the disc is darker and redder than its centre. the centre is + // this many times as bright as the disc on average + glm::vec3 const limbdarkening{0.50f, 0.60f, 0.75f}; + float const limbpeak{1.f / (1.f - 0.75f / 3.f)}; + // the brightest the centre of the disc gets: within the range of the half floats of the image, with the sky added + float const radiancemax{60000.f}; - /*float const size = std::lerp( // TODO: expose distance/scale factor from the moon object - 0.0325f, - 0.0275f, - std::clamp( Environment->m_sun.getAngle(), 0.f, 90.f ) / 90.f );*/ - model_ubs.param[0] = color; - model_ubs.param[1] = glm::vec4(glm::vec3(modelview * glm::vec4(sunvector, 1.0f)), 0.00463f /* size */); - model_ubs.param[2] = glm::vec4(0.0f, 1.0f, 1.0f, 0.0f); + auto const sunvector = Environment->m_sun.getDirection(); + glm::vec3 radiance{suncolor * overcastfactor}; + if (mainview) + { + // radiance of the disc: the illuminance of its sunlight over its solid angle. the shading lights a surface facing + // the sun with the diffuse colour times SUN_DIFFUSE_SCALE (1.5) and has no 1/pi in its lambert term + // (light_common.glsl), so its surfaces are pi times as bright as physical ones; the disc is put in the same units. + // the overcast without the clouds drawn hides the sun, where it only faded the old disc + auto const solidangle{2.f * glm::pi() * (1.f - std::cos(sunradius))}; + auto const hidden{m_weather.clouds_drawn() ? 0.f : glm::smoothstep(0.5f, 1.f, Global.Overcast)}; + auto const physical{glm::vec3{m_sunlight.diffuse} * (1.5f * glm::pi() / solidangle * Global.gfx_sun_disc)}; + radiance = glm::mix(physical, radiance, hidden); + auto const peak{std::max({radiance.r, radiance.g, radiance.b}) * limbpeak}; + if (peak > radiancemax) + radiance *= radiancemax / peak; + } + body[0] = glm::vec4{zenith, flattening(Environment->m_sun.getAngle(), sunradius), 1.f / quadmargin}; + body[1] = glm::vec4{0.f}; + body[2] = glm::vec4{limbdarkening, 0.f}; + glProgramUniform4fv(*m_celestial_shader, 0, static_cast(body.size()), glm::value_ptr(body[0])); + model_ubs.param[0] = glm::vec4{radiance, 1.f}; + model_ubs.param[1] = glm::vec4(glm::vec3(modelview * glm::vec4(sunvector, 1.0f)), std::tan(sunradius) * quadmargin); + model_ubs.param[2] = glm::vec4(0.f, 0.f, 0.f, 1.f); model_ubo->update(model_ubs); glDrawArrays(GL_TRIANGLE_STRIP, 0, 4); } // moon { + // the disc of the full moon in the picture of the phases (fx/moon, 3x3): centre and radius in texture units + glm::vec2 const discposition{0.5f, 0.5f}; + float const discradius{0.1635f}; + float const quadmargin{1.1f}; + // light of the earth on the dark part of the moon, against the light of the sun on the lit part, at the new moon + float const earthshine{0.005f}; + // extinction of the air per air mass (natural log), red, green and blue + glm::vec3 const extinction{0.10f, 0.17f, 0.30f}; + Bind_Texture(0, m_moontexture); glm::vec3 mooncolor(1.0f, 0.9490196f, 0.9058824f); - glm::vec4 color(mooncolor.r, mooncolor.g, mooncolor.b, - // fade the moon if it's near the sun in the sky, especially during the day - std::max(0.f, 1.0 - 0.5 * Global.fLuminance - 0.65 * std::max(0.f, glm::dot(Environment->m_sun.getDirection(), Environment->m_moon.getDirection()))) * overcastfactor); + // fade the moon if it's near the sun in the sky, especially during the day + float const visibility = std::max(0.f, 1.0 - 0.5 * Global.fLuminance - 0.65 * std::max(0.f, glm::dot(Environment->m_sun.getDirection(), Environment->m_moon.getDirection()))) * overcastfactor; auto const moonvector = Environment->m_moon.getDirection(); + auto const moonradius{Environment->m_moon.getAngularRadius()}; + auto const moonelevation{Environment->m_moon.getAngle()}; + auto const moonview{glm::vec3(modelview * glm::vec4(moonvector, 1.0f))}; - // choose the moon appearance variant, based on current moon phase - // NOTE: implementation specific, 8 variants are laid out in 3x3 arrangement - // from new moon onwards, top left to right bottom (last spot is left for future use, if any) - auto const moonphase = Environment->m_moon.getPhase(); - float moonu, moonv; - if (moonphase < 1.84566f) + if (Environment->m_moon.getPhase() == 50) { - moonv = 1.0f - 0.0f; - moonu = 0.0f; - } - else if (moonphase < 5.53699f) - { - moonv = 1.0f - 0.0f; - moonu = 0.333f; - } - else if (moonphase < 9.22831f) - { - moonv = 1.0f - 0.0f; - moonu = 0.667f; - } - else if (moonphase < 12.91963f) - { - moonv = 1.0f - 0.333f; - moonu = 0.0f; - } - else if (moonphase < 16.61096f) - { - moonv = 1.0f - 0.333f; - moonu = 0.333f; - } - else if (moonphase < 20.30228f) - { - moonv = 1.0f - 0.333f; - moonu = 0.667f; - } - else if (moonphase < 23.99361f) - { - moonv = 1.0f - 0.667f; - moonu = 0.0f; - } - else if (moonphase < 27.68493f) - { - moonv = 1.0f - 0.667f; - moonu = 0.333f; - } - else if (moonphase == 50) //9th slot used for easter egg - { - moonv = 1.0f - 0.667f; - moonu = 0.66f; + // 9th slot of the picture used for easter egg, drawn as it is + model_ubs.param[0] = glm::vec4(mooncolor, visibility); + model_ubs.param[1] = glm::vec4(moonview, std::tan(moonradius) * (0.1665f / discradius)); + model_ubs.param[2] = glm::vec4(0.66f, 1.0f - 0.667f, 0.333f, 0.f); } else { - moonv = 1.0f - 0.0f; - moonu = 0.0f; + // the air dims and reddens the moon low over the horizon (air mass after Kasten and Young 1989), against its colour + // high in the sky + auto const zenithangle{90.f - std::max(moonelevation, -1.f)}; + auto const airmass{1.f / (std::cos(glm::radians(zenithangle)) + 0.50572f * std::pow(std::max(96.07995f - zenithangle, 0.1f), -1.6364f))}; + auto const transmittance{glm::exp(-extinction * (airmass - 1.f))}; + // the dark part of the moon hides the stars behind it. the sky behind it isn't hidden in the twilight (its light is + // scattered in front of the moon), so this waits until the sky is dark + auto const night{1.f - glm::smoothstep(-15.f, -9.f, Environment->m_sun.getAngle())}; + // the moon is lit from the direction of the sun, which shapes its phase; the earth, full at the new moon, lights its + // dark part a little + auto const earthphase{0.5f * (1.f - std::cos(glm::radians(Environment->m_moon.getPhaseAngle())))}; + body[0] = glm::vec4{zenith, flattening(moonelevation, moonradius), 1.f / quadmargin}; + body[1] = glm::vec4{glm::normalize(glm::mat3{modelview} * Environment->m_sun.getDirection()), earthshine * earthphase}; + body[2] = glm::vec4{0.f}; + glProgramUniform4fv(*m_celestial_shader, 0, static_cast(body.size()), glm::value_ptr(body[0])); + model_ubs.param[0] = glm::vec4(mooncolor * transmittance * (visibility * Global.gfx_moon_disc), night * overcastfactor); + model_ubs.param[1] = glm::vec4(moonview, std::tan(moonradius) * quadmargin); + model_ubs.param[2] = glm::vec4(discposition, discradius, 2.f); } - - /* - float const size = std::lerp( // TODO: expose distance/scale factor from the moon object - 0.0160f, - 0.0135f, - std::clamp( Environment->m_moon.getAngle(), 0.f, 90.f ) / 90.f );*/ - - model_ubs.param[0] = color; - model_ubs.param[1] = glm::vec4(glm::vec3(modelview * glm::vec4(moonvector, 1.0f)), 0.00451f /* size */); - model_ubs.param[2] = glm::vec4(moonu, moonv, 0.333f, 0.0f); model_ubo->update(model_ubs); glDrawArrays(GL_TRIANGLE_STRIP, 0, 4); } @@ -2354,8 +2372,11 @@ bool opengl46_renderer::Render(world_environment *Environment) m_shadowcolor = glm::mix(glm::vec4{colors::shadow}, glm::vec4{colors::white}, std::clamp(-Environment->m_sun.getAngle(), 0.f, 6.f) / 6.f); if ((Environment->m_sun.getAngle() < -18.f) && (Environment->m_moon.getAngle() > 0.f)) { - // turn on moon shadows after nautical twilight, if the moon is actually up - m_shadowcolor = colors::shadow; + // turn on moon shadows after astronomical twilight, if the moon is actually up. they come in with the moonlight + // (simulation/simulationenvironment.cpp) and are as deep as it is strong: full from the gibbous moon, faint from a crescent + auto const moonshadow{std::clamp(4.f * Environment->m_moon.getBrightness(), 0.f, 1.f) * std::clamp((-18.f - Environment->m_sun.getAngle()) / 2.f, 0.f, 1.f) * + std::clamp(Environment->m_moon.getAngle() / 5.f, 0.f, 1.f)}; + m_shadowcolor = glm::mix(glm::vec4{colors::white}, glm::vec4{colors::shadow}, moonshadow); } // soften shadows depending on sky overcast factor, unless the shadows of the clouds hold back the sun where they fall if (false == m_weather.cloud_shadows()) diff --git a/shaders/celestial.frag b/shaders/celestial.frag index 7ef0756f..eee97554 100644 --- a/shaders/celestial.frag +++ b/shaders/celestial.frag @@ -1,4 +1,11 @@ +// the sun, the moon and the textured sprite of the easter egg, drawn over the sky on a quad around the direction of the body +// (celestial.vert). param[0]: colour (rgb) and alpha, param[1]: position in view space and half the size of the quad, +// param[2].xy: corner (sprite) or centre (moon) of the picture in the texture, param[2].z: its size (sprite) or the radius of +// the disc in it (moon), in texture units, param[2].w: kind of body, 0: textured sprite, 1: sun, 2: moon. +// the colour goes out premultiplied, the alpha is how much of the sky behind the body it hides (blending ONE, ONE_MINUS_SRC_ALPHA) + in vec2 f_coord; +in vec2 f_local; #texture (tex1, 0, sRGB_A) uniform sampler2D tex1; @@ -6,6 +13,14 @@ uniform sampler2D tex1; #include #include +// shape and light of the sun and the moon, set by opengl46_renderer::Render(world_environment *): +// [0].xy: direction to the zenith on the quad, [0].z: height over width of the disc (flattened by the refraction near the +// horizon), [0].w: radius of the disc in units of half the quad +// [1].xyz: direction to the sun in view space, [1].w: light of the earth on the dark part of the moon, against the light of +// the sun on the lit part +// [2].rgb: limb darkening of the sun per channel (linear law) +layout(location = 0) uniform vec4 body[3]; + layout(location = 0) out vec4 out_color; #if MOTIONBLUR_ENABLED layout(location = 1) out vec4 out_motion; @@ -13,11 +28,54 @@ layout(location = 1) out vec4 out_motion; void main() { - vec4 tex_color = texture(tex1, f_coord); + int kind = int(param[2].w + 0.5); + vec4 color = vec4(0.0); + if (kind == 0) + { + vec4 tex_color = texture(tex1, f_coord) * param[0]; + color = vec4(tex_color.rgb * tex_color.a, 0.0); + } + else + { + // the point on the disc in its radii, x along the horizon, y to the zenith, with the flattening undone + vec2 zenith = body[0].xy; + vec2 across = vec2(zenith.y, -zenith.x); + vec2 p = vec2(dot(f_local, across), dot(f_local, zenith)) / body[0].w; + p.y /= body[0].z; + float r = length(p); + float coverage = clamp((1.0 - r) / max(fwidth(r), 1e-4) + 0.5, 0.0, 1.0); + // the horizon hides what is under it + vec3 direction = normalize(param[1].xyz + vec3(f_local * param[1].w, 0.0)); + float height = (mat3(inv_view) * direction).y; + coverage *= clamp(height / max(fwidth(height), 1e-7) + 0.5, 0.0, 1.0); + // cosine of the angle between the line of sight and the surface of the sphere + float mu = sqrt(max(1.0 - r * r, 0.0)); + if (kind == 1) + { + // limb darkening, kept to the average brightness of the disc + vec3 limb = (vec3(1.0) - body[2].rgb * (1.0 - mu)) / (vec3(1.0) - body[2].rgb / 3.0); + color = vec4(param[0].rgb * limb * coverage, 0.0); + } + else + { + // the moon as a sphere lit from the direction of the sun. the frame of the disc in view space, z to the viewer + vec3 z = -normalize(param[1].xyz); + vec3 x = vec3(across, 0.0); + x = normalize(x - z * dot(x, z)); + vec3 y = cross(z, x); + vec2 q = p / max(r, 1.0); + vec3 normal = q.x * x + q.y * y + mu * z; + float mu0 = max(dot(normal, normalize(body[1].xyz)), 0.0); + // Lommel-Seeliger law: an even disc at the full moon, as the real one, no darkening towards the limb + float lit = 2.0 * mu0 / (mu0 + max(mu, 1e-3)); + vec3 albedo = texture(tex1, param[2].xy + p * param[2].z).rgb; + color = vec4(albedo * param[0].rgb * ((lit + body[1].w) * coverage), coverage * param[0].a); + } + } #if POSTFX_ENABLED - out_color = tex_color * param[0]; + out_color = color; #else - out_color = tonemap(tex_color * param[0]); + out_color = vec4(tonemap(vec4(color.rgb, 1.0)).rgb, color.a); #endif #if MOTIONBLUR_ENABLED out_motion = vec4(0.0f); diff --git a/shaders/celestial.vert b/shaders/celestial.vert index 168ba49f..f11cdb79 100644 --- a/shaders/celestial.vert +++ b/shaders/celestial.vert @@ -7,6 +7,8 @@ const vec2 vert[4] = vec2[] ); out vec2 f_coord; +// the corner of the quad, -1...1 across it +out vec2 f_local; #include @@ -14,6 +16,8 @@ void main() { gl_Position = projection * (vec4(param[1].xyz, 1.0) + vec4(vert[gl_VertexID] * param[1].w, 0.0, 0.0)); + f_local = vert[gl_VertexID]; + vec2 coord = param[2].xy; float size = param[2].z; diff --git a/simulation/simulationenvironment.cpp b/simulation/simulationenvironment.cpp index 5ac804b0..db90bffb 100644 --- a/simulation/simulationenvironment.cpp +++ b/simulation/simulationenvironment.cpp @@ -137,8 +137,17 @@ world_environment::update() { float keylightintensity; glm::vec3 keylightcolor; + // light of the night sky the key light gets on top of its own; the moon gives less of it in its thinner phases + auto lightfloor { 0.05f }; + // slightly cool moonlight, the way the eye sees a moonlit night + glm::vec3 const moonlightcolor { 0.84f, 0.91f, 1.0f }; + auto const moonlightshare { m_moon.getBrightness() }; + // the light of the sun is gone at the end of the twilight (-18 degrees), the moon's comes in over the next two degrees + // instead of all at once + auto const moonfadein { std::clamp( ( -18.f - m_sun.getAngle() ) / 2.f, 0.f, 1.f ) }; + auto moonkeylight { false }; Global.SunAngle = m_sun.getAngle(); - if (moonlightlevel > sunlightlevel && Global.SunAngle < -20) + if (moonlightlevel > sunlightlevel && Global.SunAngle < -18) { // rare situations when the moon is brighter than the sun, typically at night Global.SunAngle = m_moon.getAngle(); @@ -148,7 +157,9 @@ world_environment::update() { m_lightintensity = moonlightlevel; // if the moon is up, it overrides the twilight twilightfactor = 0.0f; - keylightcolor = glm::vec3( 1.0f, 0.9490196f, 0.7921569f ); + lightfloor *= moonlightshare; + keylightcolor = moonlightcolor * moonfadein; + moonkeylight = true; } else { // regular situation with sun as the key light @@ -173,7 +184,7 @@ world_environment::update() { // intensity combines intensity of the sun and the light reflected by the sky dome // it'd be more technically correct to have just the intensity of the sun here, // but whether it'd _look_ better is something to be tested - auto const intensity = std::min( 1.15f * ( 0.05f + keylightintensity + skydomehsv.z ), 1.25f ); + auto const intensity = std::min( 1.15f * ( lightfloor + keylightintensity + skydomehsv.z ), 1.25f ); // the impact of sun component is reduced proportionally to overcast level, as overcast increases role of ambient light auto const diffuselevel = std::lerp( keylightintensity, intensity * ( 1.0f - twilightfactor ), 1.0f - std::min( 1.f, Global.Overcast ) * 0.75f ); // ...update light colours and intensity. @@ -188,6 +199,13 @@ world_environment::update() { Global.DayLight.ambient[ 0 ] = std::lerp( skydomehsv.z, skydomecolour.r, ambienttone ) * ambientintensitynightfactor; Global.DayLight.ambient[ 1 ] = std::lerp( skydomehsv.z, skydomecolour.g, ambienttone ) * ambientintensitynightfactor; Global.DayLight.ambient[ 2 ] = std::lerp( skydomehsv.z, skydomecolour.b, ambienttone ) * ambientintensitynightfactor; + if( moonkeylight ) { + // the sky lit by the moon adds its own light, as the sky of the day does to the sunlight + auto const moonlit { moonlightcolor * ( 0.25f * moonlightlevel * moonfadein ) }; + Global.DayLight.ambient[ 0 ] += moonlit.r; + Global.DayLight.ambient[ 1 ] += moonlit.g; + Global.DayLight.ambient[ 2 ] += moonlit.b; + } Global.fLuminance = intensity; diff --git a/utilities/Globals.cpp b/utilities/Globals.cpp index ae3c1d73..06722592 100644 --- a/utilities/Globals.cpp +++ b/utilities/Globals.cpp @@ -1708,6 +1708,8 @@ global_settings::ConfigParse_gfx_sky( cParser &Parser, std::string_view const To else if( Token == "gfx.clouds.shadows.range" ) { readvalue( gfx_clouds_shadows_range, 1000.f, 100000.f ); } else if( Token == "gfx.stars" ) { readswitch( gfx_stars ); } else if( Token == "gfx.stars.brightness" ) { readvalue( gfx_stars_brightness, 0.f, 16.f ); } + else if( Token == "gfx.sun.disc" ) { readvalue( gfx_sun_disc, 0.f, 4.f ); } + else if( Token == "gfx.moon.disc" ) { readvalue( gfx_moon_disc, 0.f, 16.f ); } else { return false; } return true; @@ -2013,6 +2015,8 @@ global_settings::export_as_text( std::ostream &Output ) const { export_as_text( Output, "gfx.clouds.shadows.range", gfx_clouds_shadows_range ); export_as_text( Output, "gfx.stars", gfx_stars ); export_as_text( Output, "gfx.stars.brightness", gfx_stars_brightness ); + export_as_text( Output, "gfx.sun.disc", gfx_sun_disc ); + export_as_text( Output, "gfx.moon.disc", gfx_moon_disc ); export_as_text( Output, "python.enabled", python_enabled ); export_as_text( Output, "python.threadedupload", python_threadedupload ); export_as_text( Output, "python.uploadmain", python_uploadmain ); diff --git a/utilities/Globals.h b/utilities/Globals.h index 1a506ffd..ad03c982 100644 --- a/utilities/Globals.h +++ b/utilities/Globals.h @@ -400,6 +400,10 @@ struct global_settings { // stars from the catalogue, and their brightness bool gfx_stars { true }; float gfx_stars_brightness { 1.f }; + // brightness of the sun disc against the radiance its sunlight gives it (1: physical, held within the range of the image), + // and of the moon disc + float gfx_sun_disc { 1.f }; + float gfx_moon_disc { 1.f }; float ui_fontsize = 13.0f; float ui_scale = 1.0f;