From 272c871d15a3c4e91a2d27862d5124df662e8133 Mon Sep 17 00:00:00 2001 From: "C. J. Howard" Date: Wed, 26 May 2021 13:46:05 +0800 Subject: [PATCH] Integrate blackbody lighting with atmospheric scattering --- src/ecs/components/atmosphere-component.hpp | 12 +- src/ecs/components/blackbody-component.hpp | 15 ++ src/ecs/systems/astronomy-system.cpp | 279 +++++++++++--------- src/ecs/systems/astronomy-system.hpp | 18 +- src/game/states/play-state.cpp | 7 +- src/physics/atmosphere.hpp | 50 ++++ src/physics/constants.hpp | 24 +- src/physics/light/photometry.hpp | 4 +- src/renderer/passes/sky-pass.cpp | 32 ++- src/renderer/passes/sky-pass.hpp | 3 + 10 files changed, 290 insertions(+), 154 deletions(-) create mode 100644 src/physics/atmosphere.hpp diff --git a/src/ecs/components/atmosphere-component.hpp b/src/ecs/components/atmosphere-component.hpp index 16cc1cf..af64a30 100644 --- a/src/ecs/components/atmosphere-component.hpp +++ b/src/ecs/components/atmosphere-component.hpp @@ -27,19 +27,19 @@ namespace ecs { /// Atmosphere struct atmosphere_component { - /// Radius of the outer atmosphere, in meters. - double exosphere_radius; + /// Altitude of the outer atmosphere, in meters. + double exosphere_altitude; - /// Rayleigh scale height + /// Rayleigh scale height, in meters. double rayleigh_scale_height; - /// Mie scale height + /// Mie scale height, in meters. double mie_scale_height; - /// Rayleigh scattering coefficients + /// (Dependent) Rayleigh scattering coefficients double3 rayleigh_scattering_coefficients; - /// Mie scattering coefficients + /// (Dependent) Mie scattering coefficients double3 mie_scattering_coefficients; }; diff --git a/src/ecs/components/blackbody-component.hpp b/src/ecs/components/blackbody-component.hpp index 193bb8f..1709a2a 100644 --- a/src/ecs/components/blackbody-component.hpp +++ b/src/ecs/components/blackbody-component.hpp @@ -27,6 +27,21 @@ struct blackbody_component { /// Effective temperature, in Kelvin. double temperature; + + /// Blackbody radius, in meters. + double radius; + + /// (Dependent) Blackbody radiant flux, in watts. + double radiant_flux; + + /// (Dependent) Blackbody luminous flux, in lumens. + double luminous_flux; + + /// (Dependent) Blackbody luminous intensity, in lumens per steradian. + double luminous_intensity; + + /// (Dependent) ACEScg color + double3 color; }; } // namespace ecs diff --git a/src/ecs/systems/astronomy-system.cpp b/src/ecs/systems/astronomy-system.cpp index 687eb55..297dba8 100644 --- a/src/ecs/systems/astronomy-system.cpp +++ b/src/ecs/systems/astronomy-system.cpp @@ -36,47 +36,64 @@ namespace ecs { /** - * Calculates the optical depth between two points. + * Approximates the density of exponentially-distributed atmospheric particles between two points using the trapezoidal rule. * - * @param start Start point. - * @param end End point. - * @param sample_count Number of samples to take between the start and end points. - * @param scale_height Scale height of the atmospheric particles to measure. + * @param a Start point. + * @param b End point. + * @param r Radius of the planet. + * @param sh Scale height of the atmospheric particles. + * @param n Number of samples. */ template -T transmittance(math::vector3 start, math::vector3 end, std::size_t sample_count, T scale_height) +T optical_depth(const math::vector3& a, const math::vector3& b, T r, T sh, std::size_t n) { - const T inverse_scale_height = T(1) / -scale_height; + T inverse_sh = T(-1) / sh; - math::vector3 direction = end - start; - T distance = length(direction); - direction /= distance; + T h = math::length(b - a) / T(n); - // Calculate the distance between each sample point - T step_distance = distance / T(sample_count); + math::vector3 dy = (b - a) / T(n); + math::vector3 y = a + dy; - // Sum the atmospheric particle densities at each sample point - T total_density = 0.0; - math::vector3 sample_point = start; - for (std::size_t i = 0; i < sample_count; ++i) + T f_x = std::exp((length(a) - r) * inverse_sh); + T f_y = std::exp((length(y) - r) * inverse_sh); + T sum = (f_x + f_y); + + for (std::size_t i = 1; i < n; ++i) { - // Determine altitude of sample point - T altitude = length(sample_point); - - // Calculate atmospheric particle density at sample altitude - T density = exp(altitude * inverse_scale_height); - - // Add density to the total density - total_density += density; - - // Advance to next sample point - sample_point += direction * step_distance; + f_x = f_y; + y += dy; + f_y = std::exp((length(y) - r) * inverse_sh); + sum += (f_x + f_y); } - // Scale the total density by the step distance - return total_density * step_distance; + return sum / T(2) * h; +} + +template +math::vector3 transmittance(T depth_r, T depth_m, T depth_o, const math::vector3& beta_r, const math::vector3& beta_m) +{ + math::vector3 transmittance_r = beta_r * depth_r; + math::vector3 transmittance_m = beta_m * depth_m; + math::vector3 transmittance_o = {0, 0, 0}; + + math::vector3 t = transmittance_r + transmittance_m + transmittance_o; + t.x = std::exp(-t.x); + t.y = std::exp(-t.y); + t.z = std::exp(-t.z); + + return t; +} + +double calc_beta_r(double wavelength, double ior, double density) +{ + double wavelength2 = wavelength * wavelength; + double ior2m1 = ior * ior - 1.0; + double num = 8.0 * (math::pi * math::pi * math::pi) * ior2m1 * ior2m1; + double den = 3.0 * density * (wavelength2 * wavelength2); + return num / den; } + astronomy_system::astronomy_system(ecs::registry& registry): entity_system(registry), universal_time(0.0), @@ -86,7 +103,13 @@ astronomy_system::astronomy_system(ecs::registry& registry): reference_body_axial_rotation(0.0), sun_light(nullptr), sky_pass(nullptr) -{} +{ + registry.on_construct().connect<&astronomy_system::on_blackbody_construct>(this); + registry.on_replace().connect<&astronomy_system::on_blackbody_replace>(this); + + registry.on_construct().connect<&astronomy_system::on_atmosphere_construct>(this); + registry.on_replace().connect<&astronomy_system::on_atmosphere_replace>(this); +} void astronomy_system::update(double t, double dt) { @@ -130,59 +153,45 @@ void astronomy_system::update(double t, double dt) transform.local.translation = math::type_cast(r_topocentric); }); - // Get atmosphere component of reference body, if any - if (registry.has(reference_body)) - { - const ecs::atmosphere_component& atmosphere = registry.get(reference_body); - } - - if (sun_light != nullptr) - { - const math::vector3 sun_position_inertial = {0, 0, 0}; - const math::vector3 sun_forward_inertial = math::normalize(reference_orbit.state.r - sun_position_inertial); - const math::vector3 sun_up_inertial = {0, 0, 1}; + // Update blackbody lighting + registry.view().each( + [&](ecs::entity entity, auto& blackbody, auto& orbit) + { + // Calculate blackbody inertial basis + double3 blackbody_forward_inertial = math::normalize(reference_orbit.state.r - orbit.state.r); + double3 blackbody_up_inertial = {0, 0, 1}; - // Transform sun position, forward, and up vectors into topocentric space - const math::vector3 sun_position_topocentric = inertial_to_topocentric * sun_position_inertial; - const math::vector3 sun_forward_topocentric = inertial_to_topocentric.rotation * sun_forward_inertial; - const math::vector3 sun_up_topocentric = inertial_to_topocentric.rotation * sun_up_inertial; + // Transform blackbody inertial position and basis into topocentric space + double3 blackbody_position_topocentric = inertial_to_topocentric * orbit.state.r; + double3 blackbody_forward_topocentric = inertial_to_topocentric.rotation * blackbody_forward_inertial; + double3 blackbody_up_topocentric = inertial_to_topocentric.rotation * blackbody_up_inertial; - // Update sun light transform - sun_light->set_translation(math::type_cast(sun_position_topocentric)); - sun_light->set_rotation - ( - math::look_rotation - ( - math::type_cast(sun_forward_topocentric), - math::type_cast(sun_up_topocentric) - ) - ); + // Calculate distance from observer to blackbody + const double meters_per_au = 1.496e+11; + double blackbody_distance = math::length(blackbody_position_topocentric) * meters_per_au; - // Convert sun topocentric Cartesian coordinates to spherical coordinates - math::vector3 sun_az_el = geom::cartesian::to_spherical(ezs_to_sez * sun_position_topocentric); - sun_az_el.z = math::pi - sun_az_el.z; + // Calculate blackbody illuminance according to distance + double blackbody_illuminance = blackbody.luminous_intensity / (blackbody_distance * blackbody_distance); - //std::cout << "el: " << math::degrees(sun_az_el.y) << "; az: " << math::degrees(sun_az_el.z) << std::endl; + // Get blackbody color + double3 blackbody_color = blackbody.color; - // If the reference body has an atmosphere - if (registry.has(reference_body)) + // Get atmosphere component of reference body, if any + if (this->registry.has(reference_body)) { - // Get the atmosphere component of the reference body - const auto& atmosphere = registry.get(reference_body); + const ecs::atmosphere_component& atmosphere = this->registry.get(reference_body); - const double meters_per_au = 1.496e+11; const double earth_radius_au = 4.26352e-5; const double earth_radius_m = earth_radius_au * meters_per_au; - const double observer_altitude_m = (observer_location[0] - earth_radius_au) * meters_per_au; - // Altitude of observer in meters + // Altitude of observer in meters geom::ray sample_ray; - sample_ray.origin = {0, observer_altitude_m, 0}; - sample_ray.direction = math::normalize(sun_position_topocentric); + sample_ray.origin = {0, observer_location[0] * meters_per_au, 0}; + sample_ray.direction = math::normalize(blackbody_position_topocentric); geom::sphere exosphere; - exosphere.center = {0, -earth_radius_m, 0}; - exosphere.radius = atmosphere.exosphere_radius; + exosphere.center = {0, 0, 0}; + exosphere.radius = earth_radius_m + atmosphere.exosphere_altitude; auto intersection_result = geom::ray_sphere_intersection(sample_ray, exosphere); @@ -191,50 +200,44 @@ void astronomy_system::update(double t, double dt) double3 sample_start = sample_ray.origin; double3 sample_end = sample_ray.extrapolate(std::get<2>(intersection_result)); - double transmittance_rayleigh = transmittance(sample_start, sample_end, 16, atmosphere.rayleigh_scale_height); - double transmittance_mie = transmittance(sample_start, sample_end, 16, atmosphere.mie_scale_height); - - // Calculate attenuation due to atmospheric scattering - double3 scattering_attenuation = - atmosphere.rayleigh_scattering_coefficients * transmittance_rayleigh + - atmosphere.mie_scattering_coefficients * transmittance_mie; - scattering_attenuation.x = std::exp(-scattering_attenuation.x); - scattering_attenuation.y = std::exp(-scattering_attenuation.y); - scattering_attenuation.z = std::exp(-scattering_attenuation.z); - - double scattering_mean = (scattering_attenuation.x + scattering_attenuation.y + scattering_attenuation.z) / 3.0; - - const double sun_temperature = 5777.0; - const double sun_radius = 6.957e+8; - const double sun_surface_area = 4.0 * math::pi * sun_radius * sun_radius; - - // Calculate distance attenuation - double sun_distance_m = math::length(sun_position_topocentric) * meters_per_au; - double distance_attenuation = 1.0 / (sun_distance_m * sun_distance_m); - - double sun_luminous_flux = blackbody_luminous_flux(sun_temperature, sun_radius); - double sun_luminous_intensity = sun_luminous_flux / (4.0 * math::pi); - double sun_illuminance = sun_luminous_intensity / (sun_distance_m * sun_distance_m); - - std::cout << "distance atten: " << distance_attenuation << std::endl; - std::cout << "scatter atten: " << scattering_attenuation << std::endl; - - std::cout << "luminous flux: " << sun_luminous_flux << std::endl; - std::cout << "luminous intensity: " << sun_luminous_intensity << std::endl; - std::cout << "illuminance: " << sun_illuminance * scattering_mean << std::endl; - + double optical_depth_r = optical_depth(sample_start, sample_end, earth_radius_m, atmosphere.rayleigh_scale_height, 32); + double optical_depth_k = optical_depth(sample_start, sample_end, earth_radius_m, atmosphere.mie_scale_height, 32); + double optical_depth_o = 0.0; - // Calculate sun color - double3 color_xyz = color::cct::to_xyz(sun_temperature); - double3 color_acescg = color::xyz::to_acescg(color_xyz); + double3 attenuation = transmittance(optical_depth_r, optical_depth_k, optical_depth_o, atmosphere.rayleigh_scattering_coefficients, atmosphere.mie_scattering_coefficients); - sun_light->set_color(math::type_cast(color_acescg * scattering_attenuation)); - - sun_light->set_intensity(sun_illuminance); + // Attenuate blackbody color + blackbody_color *= attenuation; } } - } + + if (sun_light != nullptr) + { + // Update blackbody light transform + sun_light->set_translation(math::type_cast(blackbody_position_topocentric)); + sun_light->set_rotation + ( + math::look_rotation + ( + math::type_cast(blackbody_forward_topocentric), + math::type_cast(blackbody_up_topocentric) + ) + ); + + // Update blackbody light color and intensity + sun_light->set_color(math::type_cast(blackbody_color)); + sun_light->set_intensity(static_cast(blackbody_illuminance)); + + // Pass blackbody params to sky pas + if (this->sky_pass) + { + this->sky_pass->set_sun_object(sun_light); + this->sky_pass->set_sun_color(math::type_cast(blackbody.color * blackbody_illuminance)); + } + } + }); + // Update sky pass topocentric frame if (sky_pass != nullptr) { sky_pass->set_topocentric_frame @@ -245,8 +248,6 @@ void astronomy_system::update(double t, double dt) math::type_cast(inertial_to_topocentric.rotation) } ); - - sky_pass->set_sun_object(sun_light); } } @@ -307,15 +308,26 @@ void astronomy_system::set_sky_pass(::sky_pass* pass) this->sky_pass = pass; } -double astronomy_system::blackbody_luminous_flux(double t, double r) +void astronomy_system::on_blackbody_construct(ecs::registry& registry, ecs::entity entity, ecs::blackbody_component& blackbody) { + on_blackbody_replace(registry, entity, blackbody); +} + +void astronomy_system::on_blackbody_replace(ecs::registry& registry, ecs::entity entity, ecs::blackbody_component& blackbody) +{ + // Calculate surface area of spherical blackbody + const double surface_area = double(4) * math::pi * blackbody.radius * blackbody.radius; + + // Calculate radiant flux + blackbody.radiant_flux = physics::light::blackbody::radiant_flux(blackbody.temperature, surface_area); + // Blackbody spectral power distribution function - auto spd = [t](double x) -> double + auto spd = [blackbody](double x) -> double { // Convert nanometers to meters x *= double(1e-9); - return physics::light::blackbody::spectral_radiance(t, x, physics::constants::speed_of_light); + return physics::light::blackbody::spectral_radiance(blackbody.temperature, x, physics::constants::speed_of_light); }; // Luminous efficiency function (photopic) @@ -331,16 +343,41 @@ double astronomy_system::blackbody_luminous_flux(double t, double r) // Calculate luminous efficiency const double efficiency = physics::light::luminous_efficiency(spd, lef, samples.begin(), samples.end()); - // Calculate surface area of spherical blackbody - const double a = double(4) * math::pi * r * r; + // Convert radiant flux to luminous flux + blackbody.luminous_flux = physics::light::watts_to_lumens(blackbody.radiant_flux, efficiency); - // Calculate radiant flux - const double radiant_flux = physics::light::blackbody::radiant_flux(t, a); + // Calculate luminous intensity from luminous flux + blackbody.luminous_intensity = blackbody.luminous_flux / (4.0 * math::pi); - // Convert radiant flux to luminous flux - const double luminous_flux = physics::light::watts_to_lumens(radiant_flux, efficiency); + // Calculate blackbody color from temperature + double3 color_xyz = color::cct::to_xyz(blackbody.temperature); + blackbody.color = color::xyz::to_acescg(color_xyz); +} + +void astronomy_system::on_atmosphere_construct(ecs::registry& registry, ecs::entity entity, ecs::atmosphere_component& atmosphere) +{ + on_atmosphere_replace(registry, entity, atmosphere); +} + +void astronomy_system::on_atmosphere_replace(ecs::registry& registry, ecs::entity entity, ecs::atmosphere_component& atmosphere) +{ + // ACEScg wavelengths determined by matching wavelengths to XYZ, transforming XYZ to ACEScg, then selecting the max wavelengths for R, G, and B. + const double3 acescg_wavelengths_nm = {600.0, 540.0, 450.0}; + const double3 acescg_wavelengths_m = acescg_wavelengths_nm * 1.0e-9; + + // Calculate Rayleigh scattering coefficients + const double air_ior = 1.0003; + const double molecular_density = 2.545e25; + double3 beta_r; + atmosphere.rayleigh_scattering_coefficients = + { + calc_beta_r(acescg_wavelengths_m.x, air_ior, molecular_density), + calc_beta_r(acescg_wavelengths_m.y, air_ior, molecular_density), + calc_beta_r(acescg_wavelengths_m.z, air_ior, molecular_density) + }; - return luminous_flux; + // Calculate Mie scattering coefficients + atmosphere.mie_scattering_coefficients = {2.0e-6, 2.0e-6, 2.0e-6}; } } // namespace ecs diff --git a/src/ecs/systems/astronomy-system.hpp b/src/ecs/systems/astronomy-system.hpp index a7f104e..018ec85 100644 --- a/src/ecs/systems/astronomy-system.hpp +++ b/src/ecs/systems/astronomy-system.hpp @@ -26,11 +26,13 @@ #include "utility/fundamental-types.hpp" #include "physics/frame.hpp" #include "renderer/passes/sky-pass.hpp" +#include "ecs/components/blackbody-component.hpp" +#include "ecs/components/atmosphere-component.hpp" namespace ecs { /** - * Calculates apparent properties of celestial bodies relative to an observer (magnitude, angular radius, horizontal coordinates) and modifies their model component and/or light component to render them accordingly. + * Calculates apparent properties of celestial bodies relative to an observer. */ class astronomy_system: public entity_system @@ -86,14 +88,12 @@ public: void set_sky_pass(sky_pass* pass); private: - /** - * Calculates the luminous flux of a spherical blackbody in vacuum using the CIE 1931 standard observer photopic luminosity function. - * - * @param t Temperature of the blackbody, in kelvin. - * @param r Radius of the blackbody, in meters. - * @return Luminous flux, in lumens. - */ - static double blackbody_luminous_flux(double t, double r); + void on_blackbody_construct(ecs::registry& registry, ecs::entity entity, ecs::blackbody_component& blackbody); + void on_blackbody_replace(ecs::registry& registry, ecs::entity entity, ecs::blackbody_component& blackbody); + + void on_atmosphere_construct(ecs::registry& registry, ecs::entity entity, ecs::atmosphere_component& atmosphere); + void on_atmosphere_replace(ecs::registry& registry, ecs::entity entity, ecs::atmosphere_component& atmosphere); + double universal_time; double time_scale; diff --git a/src/game/states/play-state.cpp b/src/game/states/play-state.cpp index fd8cee1..6627871 100644 --- a/src/game/states/play-state.cpp +++ b/src/game/states/play-state.cpp @@ -114,7 +114,8 @@ void play_state_enter(game_context* ctx) orbit.elements.ta = math::radians(0.0); ecs::blackbody_component blackbody; - blackbody.temperature = 5772.0; + blackbody.temperature = 5777.0; + blackbody.radius = 6.957e+8; ecs::transform_component transform; transform.local = math::identity_transform; @@ -139,11 +140,9 @@ void play_state_enter(game_context* ctx) const double earth_radius_m = 6378e3; ecs::atmosphere_component atmosphere; - atmosphere.exosphere_radius = earth_radius_m + 100e3; + atmosphere.exosphere_altitude = 80e3; atmosphere.rayleigh_scale_height = 8000.0; atmosphere.mie_scale_height = 1200.0; - atmosphere.rayleigh_scattering_coefficients = {5.8e-6, 1.35e-5, 3.31e-5}; - atmosphere.mie_scattering_coefficients = {2e-6, 2e-6, 2e-6}; ecs::transform_component transform; transform.local = math::identity_transform; diff --git a/src/physics/atmosphere.hpp b/src/physics/atmosphere.hpp new file mode 100644 index 0000000..37f1b7a --- /dev/null +++ b/src/physics/atmosphere.hpp @@ -0,0 +1,50 @@ +/* + * Copyright (C) 2021 Christopher J. Howard + * + * This file is part of Antkeeper source code. + * + * Antkeeper source code is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * Antkeeper source code is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with Antkeeper source code. If not, see . + */ + +#ifndef ANTKEEPER_PHYSICS_ATMOSPHERE_HPP +#define ANTKEEPER_PHYSICS_ATMOSPHERE_HPP + +#include "physics/constants.hpp" +#include + +namespace physics { + +/// Atmosphere-related functions. +namespace atmosphere { + +/** + * Calculates the density of exponentially-distributed atmospheric particles at a given altitude. + * + * @param d0 Density at sea level. + * @param z Height above sea level. + * @param sh Scale height of the particle type. + * @return Particle density at altitude. + * + * @see https://en.wikipedia.org/wiki/Scale_height + */ +T density(T d0, T z, T sh) +{ + return d0 * std::exp(-z / sh); +} + +} // namespace atmosphere + +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_ATMOSPHERE_HPP diff --git a/src/physics/constants.hpp b/src/physics/constants.hpp index a13098f..ccf68d8 100644 --- a/src/physics/constants.hpp +++ b/src/physics/constants.hpp @@ -25,18 +25,34 @@ namespace physics { /// Physical constants namespace constants { +/** + * Avogadro number (N). + * + * @see https://physics.nist.gov/cgi-bin/cuu/Value?na + */ +template +constexpr T avogadro = T(6.02214076e+23); + /** * Boltzmann constant (k), in joule per kelvin. * - * @see https://physics.nist.gov/cgi-bin/cuu/Value?k|search_for=universal_in! + * @see https://physics.nist.gov/cgi-bin/cuu/Value?k */ template constexpr T boltzmann = T(1.380649e-23); +/** + * Molar gas constant (R), in joule per kelvin per mole. + * + * @see https://en.wikipedia.org/wiki/Gas_constant + */ +template +constexpr T gas = T(8.31446261815324); + /** * Gravitational constant (G), in cubic meter per second squared per kilogram. * - * @see https://physics.nist.gov/cgi-bin/cuu/Value?G|search_for=universal_in! + * @see https://physics.nist.gov/cgi-bin/cuu/Value?G */ template constexpr T gravitational = T(6.67430e-11); @@ -44,7 +60,7 @@ constexpr T gravitational = T(6.67430e-11); /** * Planck constant (h), in joule per hertz. * - * @see https://physics.nist.gov/cgi-bin/cuu/Value?h|search_for=universal_in! + * @see https://physics.nist.gov/cgi-bin/cuu/Value?h */ template constexpr T planck = T(6.62607015e-34); @@ -60,7 +76,7 @@ constexpr T stefan_boltzmann = T(5.670374419184429453970996731889230875840122970 /** * Speed of light in vacuum (c), in meters per second. * - * @see https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=universal_in! + * @see https://physics.nist.gov/cgi-bin/cuu/Value?c */ template constexpr T speed_of_light = T(299792458); diff --git a/src/physics/light/photometry.hpp b/src/physics/light/photometry.hpp index 3ea61af..bd86ec0 100644 --- a/src/physics/light/photometry.hpp +++ b/src/physics/light/photometry.hpp @@ -45,8 +45,8 @@ T luminous_efficiency(UnaryOp1 spd, UnaryOp2 lef, InputIt first, InputIt last) return spd(x) * lef(x); }; - const T num = math::quadrature::trapezoid(spd_lef, first, last); - const T den = math::quadrature::trapezoid(spd, first, last); + const T num = math::quadrature::simpson(spd_lef, first, last); + const T den = math::quadrature::simpson(spd, first, last); return num / den; } diff --git a/src/renderer/passes/sky-pass.cpp b/src/renderer/passes/sky-pass.cpp index 9918161..9dcef53 100644 --- a/src/renderer/passes/sky-pass.cpp +++ b/src/renderer/passes/sky-pass.cpp @@ -43,6 +43,7 @@ #include "geom/cartesian.hpp" #include "geom/spherical.hpp" #include "physics/orbit/orbit.hpp" +#include "physics/light/photometry.hpp" #include #include #include @@ -68,6 +69,7 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe julian_day_tween(0.0, math::lerp), horizon_color_tween(float3{0.0f, 0.0f, 0.0f}, math::lerp), zenith_color_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp), + sun_color_tween(float3{1.0f, 1.0f, 1.0f}, math::lerp), topocentric_frame_translation({0, 0, 0}, math::lerp), topocentric_frame_rotation(math::quaternion::identity(), math::nlerp), sun_object(nullptr) @@ -122,17 +124,17 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe // Calculate XYZ color from color temperature double3 color_xyz = color::cct::to_xyz(cct); - // Transform XYZ from (assumed) D65 illuminant to ACES illuminant. - //color_xyz = color::xyz::cat::d65_to_aces(color_xyz); - // Transform XYZ color to ACEScg colorspace double3 color_acescg = color::xyz::to_acescg(color_xyz); - // Convert apparent magnitude to lux + // Convert apparent magnitude to irradiance W/m2 double vmag_lux = astro::vmag_to_lux(vmag); - // Normalized color luminance and scale by apparent magnitude - double3 scaled_color = color_acescg * vmag_lux; + // Convert irradiance to illuminance (using luminous efficiency of sun) + double illuminance = physics::light::watts_to_lumens(vmag_lux, 0.13); + + // Scale color by illuminance + double3 scaled_color = color_acescg * illuminance; // Build vertex *(star_vertex++) = static_cast(position_inertial.x); @@ -204,6 +206,8 @@ void sky_pass::render(render_context* context) const float julian_day = julian_day_tween.interpolate(context->alpha); float3 horizon_color = horizon_color_tween.interpolate(context->alpha); float3 zenith_color = zenith_color_tween.interpolate(context->alpha); + float3 sun_color = sun_color_tween.interpolate(context->alpha); + // Construct tweened inertial to topocentric frame physics::frame topocentric_frame = @@ -233,6 +237,8 @@ void sky_pass::render(render_context* context) const horizon_color_input->upload(horizon_color); if (zenith_color_input) zenith_color_input->upload(zenith_color); + if (sun_color_input) + sun_color_input->upload(sun_color); if (mouse_input) mouse_input->upload(mouse_position); if (resolution_input) @@ -267,11 +273,14 @@ void sky_pass::render(render_context* context) const rasterizer->draw_arrays(*sky_model_vao, sky_model_drawing_mode, sky_model_start_index, sky_model_index_count); } + glEnable(GL_BLEND); + //glBlendFunc(GL_SRC_ALPHA, GL_ONE); + glBlendFunc(GL_ONE, GL_ONE); + // Draw moon model if (moon_position.y >= -moon_angular_radius) { - glEnable(GL_BLEND); - glBlendFunc(GL_SRC_ALPHA, GL_ONE); + float moon_distance = (clip_near + clip_far) * 0.5f; float moon_radius = moon_angular_radius * moon_distance; @@ -370,6 +379,7 @@ void sky_pass::set_sky_model(const model* model) model_view_projection_input = sky_shader_program->get_input("model_view_projection"); horizon_color_input = sky_shader_program->get_input("horizon_color"); zenith_color_input = sky_shader_program->get_input("zenith_color"); + sun_color_input = sky_shader_program->get_input("sun_color"); mouse_input = sky_shader_program->get_input("mouse"); resolution_input = sky_shader_program->get_input("resolution"); time_input = sky_shader_program->get_input("time"); @@ -437,6 +447,7 @@ void sky_pass::update_tweens() zenith_color_tween.update(); topocentric_frame_translation.update(); topocentric_frame_rotation.update(); + sun_color_tween.update(); } void sky_pass::set_time_of_day(float time) @@ -503,6 +514,11 @@ void sky_pass::set_zenith_color(const float3& color) zenith_color_tween[1] = color; } +void sky_pass::set_sun_color(const float3& color) +{ + sun_color_tween[1] = color; +} + void sky_pass::handle_event(const mouse_moved_event& event) { mouse_position = {static_cast(event.x), static_cast(event.y)}; diff --git a/src/renderer/passes/sky-pass.hpp b/src/renderer/passes/sky-pass.hpp index 5b0819d..cd7b724 100644 --- a/src/renderer/passes/sky-pass.hpp +++ b/src/renderer/passes/sky-pass.hpp @@ -55,6 +55,7 @@ public: void set_sky_model(const model* model); void set_horizon_color(const float3& color); void set_zenith_color(const float3& color); + void set_sun_color(const float3& color); void set_time_of_day(float time); void set_blue_noise_map(const gl::texture_2d* texture); void set_sky_gradient(const gl::texture_2d* texture, const gl::texture_2d* texture2); @@ -91,6 +92,7 @@ private: const gl::shader_input* sky_gradient_input; const gl::shader_input* sky_gradient2_input; const gl::shader_input* exposure_input; + const gl::shader_input* sun_color_input; gl::shader_program* moon_shader_program; const gl::shader_input* moon_model_view_projection_input; @@ -133,6 +135,7 @@ private: tween julian_day_tween; tween horizon_color_tween; tween zenith_color_tween; + tween sun_color_tween; tween topocentric_frame_translation; tween> topocentric_frame_rotation;