diff --git a/src/entity/components/atmosphere.hpp b/src/entity/components/atmosphere.hpp index 7e9330d..7e418c3 100644 --- a/src/entity/components/atmosphere.hpp +++ b/src/entity/components/atmosphere.hpp @@ -49,14 +49,32 @@ struct atmosphere /// Mie phase function anisotropy factor. double mie_anisotropy; - /// (Dependent) Rayleigh scattering coefficients at sea level. + /// Airglow, in lux. + double airglow; + + /// Molar concentration of air at sea level, in mol/m-3. + double air_concentration; + + /// Concentration of ozone in the atmosphere, unitless. + double ozone_concentration; + + /// Elevation of the lower limit of the ozone triangular distribution, in meters. + double ozone_lower_limit; + + /// Elevation of the upper limit of the ozone triangular distribution, in meters. + double ozone_upper_limit; + + /// Elevation of the mode of the ozone triangular distribution, in meters. + double ozone_mode; + + /// (Dependent) Rayleigh scattering coefficients. double3 rayleigh_scattering; - /// (Dependent) Mie scattering coefficients at sea level. + /// (Dependent) Mie scattering coefficients. double3 mie_scattering; - /// Airglow, in lux. - double airglow; + /// (Dependent) Ozone absorption coefficients. + double3 ozone_absorption; }; } // namespace component diff --git a/src/entity/systems/astronomy.cpp b/src/entity/systems/astronomy.cpp index 79ceff4..50c472f 100644 --- a/src/entity/systems/astronomy.cpp +++ b/src/entity/systems/astronomy.cpp @@ -30,7 +30,7 @@ #include "physics/light/photometry.hpp" #include "physics/light/luminosity.hpp" #include "physics/light/refraction.hpp" -#include "physics/atmosphere.hpp" +#include "physics/gas/atmosphere.hpp" #include "geom/cartesian.hpp" #include "astro/apparent-size.hpp" #include "geom/solid-angle.hpp" @@ -41,11 +41,11 @@ namespace entity { namespace system { 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(T depth_r, T depth_m, T depth_o, const math::vector3& beta_r, const math::vector3& beta_m, const math::vector3& beta_o) { math::vector3 transmittance_r = beta_r * depth_r; - math::vector3 transmittance_m = beta_m * 1.1 * depth_m; - math::vector3 transmittance_o = {0, 0, 0}; + math::vector3 transmittance_m = beta_m * (T(10)/T(9)) * depth_m; + math::vector3 transmittance_o = beta_o * depth_o; math::vector3 t = transmittance_r + transmittance_m + transmittance_o; t.x = std::exp(-t.x); @@ -218,7 +218,7 @@ void astronomy::update(double t, double dt) { const entity::component::atmosphere& reference_atmosphere = registry.get(reference_entity); - // Altitude of observer in meters + // Altitude of observer in meters geom::ray sample_ray; sample_ray.origin = {0, reference_body.radius + observer_location[0], 0}; sample_ray.direction = math::normalize(blackbody_position_eus); @@ -234,11 +234,11 @@ void astronomy::update(double t, double dt) double3 sample_start = sample_ray.origin; double3 sample_end = sample_ray.extrapolate(std::get<2>(intersection_result)); - double optical_depth_r = physics::atmosphere::optical_depth(sample_start, sample_end, reference_body.radius, reference_atmosphere.rayleigh_scale_height, 32); - double optical_depth_k = physics::atmosphere::optical_depth(sample_start, sample_end, reference_body.radius, reference_atmosphere.mie_scale_height, 32); - double optical_depth_o = 0.0; + double optical_depth_r = physics::gas::atmosphere::optical_depth_exp(sample_start, sample_end, reference_body.radius, reference_atmosphere.rayleigh_scale_height, 16); + double optical_depth_m = physics::gas::atmosphere::optical_depth_exp(sample_start, sample_end, reference_body.radius, reference_atmosphere.mie_scale_height, 16); + double optical_depth_o = physics::gas::atmosphere::optical_depth_tri(sample_start, sample_end, reference_body.radius, reference_atmosphere.ozone_lower_limit, reference_atmosphere.ozone_upper_limit, reference_atmosphere.ozone_mode, 16); - atmospheric_transmittance = transmittance(optical_depth_r, optical_depth_k, optical_depth_o, reference_atmosphere.rayleigh_scattering, reference_atmosphere.mie_scattering); + atmospheric_transmittance = transmittance(optical_depth_r, optical_depth_m, optical_depth_o, reference_atmosphere.rayleigh_scattering, reference_atmosphere.mie_scattering, reference_atmosphere.ozone_absorption); } // Add airglow to sky light illuminance @@ -419,9 +419,17 @@ void astronomy::update(double t, double dt) { const entity::component::atmosphere& reference_atmosphere = registry.get(reference_entity); - sky_pass->set_scale_heights(reference_atmosphere.rayleigh_scale_height, reference_atmosphere.mie_scale_height); + sky_pass->set_particle_distributions + ( + static_cast(reference_atmosphere.rayleigh_scale_height), + static_cast(reference_atmosphere.mie_scale_height), + static_cast(reference_atmosphere.ozone_lower_limit), + static_cast(reference_atmosphere.ozone_upper_limit), + static_cast(reference_atmosphere.ozone_mode) + ); sky_pass->set_scattering_coefficients(math::type_cast(reference_atmosphere.rayleigh_scattering), math::type_cast(reference_atmosphere.mie_scattering)); sky_pass->set_mie_anisotropy(reference_atmosphere.mie_anisotropy); + sky_pass->set_absorption_coefficients(math::type_cast(reference_atmosphere.ozone_absorption)); sky_pass->set_atmosphere_radii(reference_body.radius, reference_body.radius + reference_atmosphere.exosphere_altitude); } } diff --git a/src/entity/systems/atmosphere.cpp b/src/entity/systems/atmosphere.cpp index 8d3b15b..f917386 100644 --- a/src/entity/systems/atmosphere.cpp +++ b/src/entity/systems/atmosphere.cpp @@ -18,15 +18,18 @@ */ #include "entity/systems/atmosphere.hpp" -#include "physics/atmosphere.hpp" +#include "physics/gas/ozone.hpp" +#include "physics/gas/atmosphere.hpp" +#include "physics/number-density.hpp" +#include "color/srgb.hpp" namespace entity { namespace system { atmosphere::atmosphere(entity::registry& registry): updatable(registry), - rgb_wavelengths_nm{0, 0, 0}, - rgb_wavelengths_m{0, 0, 0} + rgb_wavelengths{0, 0, 0}, + rgb_ozone_cross_sections{0, 0, 0} { registry.on_construct().connect<&atmosphere::on_atmosphere_construct>(this); registry.on_replace().connect<&atmosphere::on_atmosphere_replace>(this); @@ -37,8 +40,12 @@ void atmosphere::update(double t, double dt) void atmosphere::set_rgb_wavelengths(const double3& wavelengths) { - rgb_wavelengths_nm = wavelengths; - rgb_wavelengths_m = wavelengths * 1e-9; + rgb_wavelengths = wavelengths; +} + +void atmosphere::set_rgb_ozone_cross_sections(const double3& cross_sections) +{ + rgb_ozone_cross_sections = cross_sections; } void atmosphere::update_coefficients(entity::id entity_id) @@ -51,25 +58,40 @@ void atmosphere::update_coefficients(entity::id entity_id) component::atmosphere& atmosphere = registry.get(entity_id); // Calculate polarization factors - const double rayleigh_polarization = physics::atmosphere::polarization(atmosphere.index_of_refraction, atmosphere.rayleigh_density); - const double mie_polarization = physics::atmosphere::polarization(atmosphere.index_of_refraction, atmosphere.mie_density); + const double rayleigh_polarization = physics::gas::atmosphere::polarization(atmosphere.index_of_refraction, atmosphere.rayleigh_density); + const double mie_polarization = physics::gas::atmosphere::polarization(atmosphere.index_of_refraction, atmosphere.mie_density); - // Calculate Rayleigh scattering coefficients - atmosphere.rayleigh_scattering = + // Calculate Rayleigh scattering coefficients for sRGB wavelengths + const double3 rayleigh_scattering_srgb = { - physics::atmosphere::scattering_rayleigh(rgb_wavelengths_m.x, atmosphere.rayleigh_density, rayleigh_polarization), - physics::atmosphere::scattering_rayleigh(rgb_wavelengths_m.y, atmosphere.rayleigh_density, rayleigh_polarization), - physics::atmosphere::scattering_rayleigh(rgb_wavelengths_m.z, atmosphere.rayleigh_density, rayleigh_polarization) + physics::gas::atmosphere::scattering_rayleigh(rgb_wavelengths.x, atmosphere.rayleigh_density, rayleigh_polarization), + physics::gas::atmosphere::scattering_rayleigh(rgb_wavelengths.y, atmosphere.rayleigh_density, rayleigh_polarization), + physics::gas::atmosphere::scattering_rayleigh(rgb_wavelengths.z, atmosphere.rayleigh_density, rayleigh_polarization) }; + // Transform Rayleigh scattering coefficients from sRGB to ACEScg + atmosphere.rayleigh_scattering = color::srgb::to_acescg(rayleigh_scattering_srgb); + // Calculate Mie scattering coefficients - const double mie_scattering = physics::atmosphere::scattering_mie(atmosphere.mie_density, mie_polarization); + const double mie_scattering = physics::gas::atmosphere::scattering_mie(atmosphere.mie_density, mie_polarization); atmosphere.mie_scattering = { mie_scattering, mie_scattering, mie_scattering }; + + // Calculate ozone absorption coefficients for sRGB wavelengths + const double n_air = physics::number_density(atmosphere.air_concentration); + const double3 ozone_absorption_srgb = + { + physics::gas::ozone::absorption(rgb_ozone_cross_sections.x, n_air, atmosphere.ozone_concentration), + physics::gas::ozone::absorption(rgb_ozone_cross_sections.y, n_air, atmosphere.ozone_concentration), + physics::gas::ozone::absorption(rgb_ozone_cross_sections.z, n_air, atmosphere.ozone_concentration) + }; + + // Transform ozone absorption coefficients from sRGB to ACEScg + atmosphere.ozone_absorption = color::srgb::to_acescg(ozone_absorption_srgb); } void atmosphere::on_atmosphere_construct(entity::registry& registry, entity::id entity_id, entity::component::atmosphere& atmosphere) diff --git a/src/entity/systems/atmosphere.hpp b/src/entity/systems/atmosphere.hpp index 03e600b..68aebe3 100644 --- a/src/entity/systems/atmosphere.hpp +++ b/src/entity/systems/atmosphere.hpp @@ -42,18 +42,25 @@ public: /** * Sets the wavelengths of red, green, and blue light. * - * @param wavelengths Vector containing the wavelengths of red (x), green (y), and blue (z) light, in nanometers. + * @param wavelengths Vector containing the wavelengths of red (x), green (y), and blue (z) light, in meters. */ void set_rgb_wavelengths(const double3& wavelengths); + /** + * Sets ozone cross sections for red, green, and blue wavelengths. + * + * @param cross_sections Vector containing the ozone cross sections of red (x), green (y), and blue (z) light, in m-2/molecule. + */ + void set_rgb_ozone_cross_sections(const double3& cross_sections); + private: void update_coefficients(entity::id entity_id); void on_atmosphere_construct(entity::registry& registry, entity::id entity_id, entity::component::atmosphere& atmosphere); void on_atmosphere_replace(entity::registry& registry, entity::id entity_id, entity::component::atmosphere& atmosphere); - double3 rgb_wavelengths_nm; - double3 rgb_wavelengths_m; + double3 rgb_wavelengths; + double3 rgb_ozone_cross_sections; }; } // namespace system diff --git a/src/entity/systems/blackbody.cpp b/src/entity/systems/blackbody.cpp index 8cc01c0..5eb4cd4 100644 --- a/src/entity/systems/blackbody.cpp +++ b/src/entity/systems/blackbody.cpp @@ -27,9 +27,7 @@ namespace entity { namespace system { blackbody::blackbody(entity::registry& registry): - updatable(registry), - rgb_wavelengths_nm{0, 0, 0}, - rgb_wavelengths_m{0, 0, 0} + updatable(registry) { // Construct a range of sample wavelengths in the visible spectrum visible_wavelengths_nm.resize(780 - 280); @@ -45,12 +43,6 @@ blackbody::blackbody(entity::registry& registry): void blackbody::update(double t, double dt) {} -void blackbody::set_rgb_wavelengths(const double3& wavelengths) -{ - rgb_wavelengths_nm = wavelengths; - rgb_wavelengths_m = wavelengths * 1e-9; -} - void blackbody::update_luminance(entity::id entity_id) { // Abort if entity has no blackbody component diff --git a/src/entity/systems/blackbody.hpp b/src/entity/systems/blackbody.hpp index 6c3fa3b..7512398 100644 --- a/src/entity/systems/blackbody.hpp +++ b/src/entity/systems/blackbody.hpp @@ -41,13 +41,6 @@ public: virtual void update(double t, double dt); - /** - * Sets the wavelengths of red, green, and blue light. - * - * @param wavelengths Vector containing the wavelengths of red (x), green (y), and blue (z) light, in nanometers. - */ - void set_rgb_wavelengths(const double3& wavelengths); - private: void update_luminance(entity::id entity_id); @@ -57,8 +50,6 @@ private: void on_celestial_body_construct(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body); void on_celestial_body_replace(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& celestial_body); - double3 rgb_wavelengths_nm; - double3 rgb_wavelengths_m; std::vector visible_wavelengths_nm; }; diff --git a/src/game/state/boot.cpp b/src/game/state/boot.cpp index 5254e58..b3d9a8c 100644 --- a/src/game/state/boot.cpp +++ b/src/game/state/boot.cpp @@ -92,6 +92,8 @@ #include "game/menu.hpp" #include "game/graphics.hpp" #include "utility/timestamp.hpp" +#include "color/color.hpp" +#include "physics/gas/ozone.hpp" #include #include #include @@ -832,8 +834,19 @@ void boot::setup_systems() const auto& viewport_dimensions = ctx.app->get_viewport_dimensions(); float4 viewport = {0.0f, 0.0f, static_cast(viewport_dimensions[0]), static_cast(viewport_dimensions[1])}; - // RGB wavelengths determined by matching wavelengths to XYZ, transforming XYZ to ACEScg, then selecting the max wavelengths for R, G, and B. - const double3 rgb_wavelengths_nm = {602.224, 541.069, 448.143}; + // RGB wavelengths + const double3 rgb_wavelengths_nm = {680, 550, 440}; + const double3 rgb_wavelengths_m = rgb_wavelengths_nm * 1e-9; + + // Ozone cross sections + double3 rgb_ozone_cross_sections_293k = + { + physics::gas::ozone::cross_section_680nm_293k, + physics::gas::ozone::cross_section_550nm_293k, + physics::gas::ozone::cross_section_440nm_293k + }; + + rgb_ozone_cross_sections_293k = rgb_ozone_cross_sections_293k; // Setup terrain system ctx.terrain_system = new entity::system::terrain(*ctx.entity_registry); @@ -891,11 +904,11 @@ void boot::setup_systems() // Setup blackbody system ctx.blackbody_system = new entity::system::blackbody(*ctx.entity_registry); - ctx.blackbody_system->set_rgb_wavelengths(rgb_wavelengths_nm); // Setup atmosphere system ctx.atmosphere_system = new entity::system::atmosphere(*ctx.entity_registry); - ctx.atmosphere_system->set_rgb_wavelengths(rgb_wavelengths_nm); + ctx.atmosphere_system->set_rgb_wavelengths(rgb_wavelengths_m); + ctx.atmosphere_system->set_rgb_ozone_cross_sections(rgb_ozone_cross_sections_293k); // Setup astronomy system ctx.astronomy_system = new entity::system::astronomy(*ctx.entity_registry); diff --git a/src/game/state/nuptial-flight.cpp b/src/game/state/nuptial-flight.cpp index a4589e1..57713f4 100644 --- a/src/game/state/nuptial-flight.cpp +++ b/src/game/state/nuptial-flight.cpp @@ -98,7 +98,7 @@ nuptial_flight::nuptial_flight(game::context& ctx): game::world::cosmogenesis(ctx); // Create boids - for (int i = 0; i < 50; ++i) + for (int i = 0; i < 100; ++i) { entity::id boid_eid = ctx.entity_registry->create(); @@ -143,7 +143,7 @@ nuptial_flight::nuptial_flight(game::context& ctx): } // Load biome - game::load::biome(ctx, "debug.bio"); + game::load::biome(ctx, "desert-scrub.bio"); // Set world time game::world::set_time(ctx, 2022, 6, 21, 12, 0, 0.0); @@ -353,7 +353,7 @@ void nuptial_flight::enable_keeper_controls() bool mouse_look_toggle = false; ctx.mouse_look = false; const double time_scale = 60.0; - const double ff_time_scale = 50000.0; + const double ff_time_scale = time_scale * 200; if (ctx.config->contains("mouse_tilt_sensitivity")) mouse_tilt_sensitivity = math::radians((*ctx.config)["mouse_tilt_sensitivity"].get()); @@ -668,22 +668,22 @@ void nuptial_flight::enable_keeper_controls() } ); - ctx.controls["increase_exposure"]->set_activated_callback + ctx.controls["increase_exposure"]->set_active_callback ( - [&ctx = this->ctx]() + [&ctx = this->ctx](float) { //ctx.astronomy_system->set_exposure_offset(ctx.astronomy_system->get_exposure_offset() - 1.0f); - ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() + 1.0f); + ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() + 3.0f * (1.0f / 60.0f)); ctx.logger->log("EV100: " + std::to_string(ctx.surface_camera->get_exposure())); } ); - ctx.controls["decrease_exposure"]->set_activated_callback + ctx.controls["decrease_exposure"]->set_active_callback ( - [&ctx = this->ctx]() + [&ctx = this->ctx](float) { //ctx.astronomy_system->set_exposure_offset(ctx.astronomy_system->get_exposure_offset() + 1.0f); - ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() - 1.0f); + ctx.surface_camera->set_exposure(ctx.surface_camera->get_exposure() - 3.0f * (1.0f / 60.0f)); ctx.logger->log("EV100: " + std::to_string(ctx.surface_camera->get_exposure())); } ); @@ -735,7 +735,7 @@ void nuptial_flight::enable_ant_controls() bool gamepad_invert_tilt = false; bool gamepad_invert_pan = false; const double time_scale = 60.0; - const double ff_time_scale = 50000.0; + const double ff_time_scale = time_scale * 200; if (ctx.config->contains("mouse_tilt_sensitivity")) mouse_tilt_sensitivity = math::radians((*ctx.config)["mouse_tilt_sensitivity"].get()); diff --git a/src/physics/atmosphere.hpp b/src/physics/gas/atmosphere.hpp similarity index 53% rename from src/physics/atmosphere.hpp rename to src/physics/gas/atmosphere.hpp index 2e26cd9..9420fc1 100644 --- a/src/physics/atmosphere.hpp +++ b/src/physics/gas/atmosphere.hpp @@ -17,35 +17,20 @@ * along with Antkeeper source code. If not, see . */ -#ifndef ANTKEEPER_PHYSICS_ATMOSPHERE_HPP -#define ANTKEEPER_PHYSICS_ATMOSPHERE_HPP +#ifndef ANTKEEPER_PHYSICS_GAS_ATMOSPHERE_HPP +#define ANTKEEPER_PHYSICS_GAS_ATMOSPHERE_HPP #include "physics/constants.hpp" #include "math/constants.hpp" +#include #include namespace physics { +namespace gas { /// 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 - * @see https://en.wikipedia.org/wiki/Barometric_formula - */ -template -T density(T d0, T z, T sh) -{ - return d0 * std::exp(-z / sh); -} - /** * Calculates a particle polarizability factor used in computing scattering coefficients. * @@ -66,10 +51,10 @@ T polarization(T ior, T density) } /** - * Calculates a Rayleigh scattering coefficient at sea level (wavelength-dependent). + * Calculates a Rayleigh scattering coefficient (wavelength-dependent). * * @param wavelength Wavelength of light, in meters. - * @param density Molecular density of Rayleigh particles at sea level. + * @param density Molecular density of Rayleigh particles. * @param polarization Rayleigh particle polarization factor. * * @see atmosphere::polarization @@ -81,15 +66,17 @@ template T scattering_rayleigh(T wavelength, T density, T polarization) { const T wavelength2 = wavelength * wavelength; - return T(4) * math::pi * density / (wavelength2 * wavelength2) * polarization; + return math::four_pi * density / (wavelength2 * wavelength2) * polarization; } /** - * Calculates a Mie scattering coefficient at sea level (wavelength-independent). + * Calculates a Mie scattering coefficient (wavelength-independent). * - * @param density Molecular density of Mie particles at sea level. + * @param density Molecular density of Mie particles. * @param polarization Mie particle polarization factor. * + * @return Mie scattering coefficient. + * * @see atmosphere::polarization * * @see Elek, Oskar. (2009). Rendering Parametrizable Planetary Atmospheres with Multiple Scattering in Real-Time. @@ -98,7 +85,22 @@ T scattering_rayleigh(T wavelength, T density, T polarization) template T scattering_mie(T density, T polarization) { - return T(4) * math::pi * density * polarization; + return math::four_pi * density * polarization; +} + +/** + * Calculates a Mie absorption coefficient (wavelength-independent). + * + * @param scattering Mie scattering coefficient. + * + * @return Mie absorption coefficient. + * + * @see Bruneton, E. and Neyret, F. (2008), Precomputed Atmospheric Scattering. Computer Graphics Forum, 27: 1079-1086. https://doi.org/10.1111/j.1467-8659.2008.01245.x + */ +template +T absorption_mie(T scattering) +{ + return scattering / T(9); } /** @@ -138,7 +140,7 @@ T albedo(T s, T e) * @return Optical depth between @p a and @p b. */ template -T optical_depth(const math::vector3& a, const math::vector3& b, T r, T sh, std::size_t n) +T optical_depth_exp(const math::vector3& a, const math::vector3& b, T r, T sh, std::size_t n) { sh = T(-1) / sh; @@ -162,8 +164,95 @@ T optical_depth(const math::vector3& a, const math::vector3& b, T r, T sh, return sum / T(2) * h; } +/** + * Approximates the optical depth of triangularly-distributed atmospheric particles between two points using the trapezoidal rule. + * + * @param p0 Start point. + * @param p1 End point. + * @param r Radius of the planet. + * @param a Distribution lower limit. + * @param b Distribution upper limit. + * @param c Distribution upper mode. + * @param n Number of samples. + * @return Optical depth between @p a and @p b. + */ +template +T optical_depth_tri(const math::vector3& p0, const math::vector3& p1, T r, T a, T b, T c, std::size_t n) +{ + a = T(1) / (a - c); + b = T(1) / (b - c); + + const T h = math::length(p1 - p0) / T(n); + + math::vector3 dy = (p1 - p0) / T(n); + math::vector3 y = p0 + dy; + + T z = math::length(p0) - r; + T f_x = std::max(T(0), std::max(T(0), c - z) * a - std::max(T(0), z - c) * b + T(1)); + + z = math::length(y) - r; + T f_y = std::max(T(0), std::max(T(0), c - z) * a - std::max(T(0), z - c) * b + T(1)); + T sum = (f_x + f_y); + + for (std::size_t i = 1; i < n; ++i) + { + f_x = f_y; + y += dy; + + z = math::length(y) - r; + f_y = std::max(T(0), std::max(T(0), c - z) * a - std::max(T(0), z - c) * b + T(1)); + + sum += (f_x + f_y); + } + + return sum / T(2) * h; +} + +/// Atmospheric density functions. +namespace density { + + /** + * Calculates the density of exponentially-distributed atmospheric particles at a given elevation. + * + * @param d0 Density at sea level. + * @param z Height above sea level. + * @param sh Scale height of the particle type. + * + * @return Particle density at elevation @p z. + * + * @see https://en.wikipedia.org/wiki/Barometric_formula + * @see https://en.wikipedia.org/wiki/Scale_height + */ + template + T exponential(T d0, T z, T sh) + { + return d0 * std::exp(-z / sh); + } + + /** + * Calculates the density of triangularly-distributed atmospheric particles at a given elevation. + * + * @param d0 Density at sea level. + * @param z Height above sea level. + * @param a Distribution lower limit. + * @param b Distribution upper limit. + * @param c Distribution mode. + * + * @return Particle density at elevation @p z. + * + * @see https://en.wikipedia.org/wiki/Triangular_distribution + */ + template + T triangular(T d0, T z, T a, T b, T c) + { + return d0 * max(T(0), max(T(0), c - z) / (a - c) - max(T(0), z - c) / (b - c) + T(1)); + } + +} // namespace density + } // namespace atmosphere +} // namespace gas } // namespace physics -#endif // ANTKEEPER_PHYSICS_ATMOSPHERE_HPP +#endif // ANTKEEPER_PHYSICS_GAS_ATMOSPHERE_HPP diff --git a/src/physics/gas/gas.hpp b/src/physics/gas/gas.hpp new file mode 100644 index 0000000..44e8dd4 --- /dev/null +++ b/src/physics/gas/gas.hpp @@ -0,0 +1,33 @@ +/* + * 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_GAS_HPP +#define ANTKEEPER_PHYSICS_GAS_HPP + +namespace physics { + +/// Gas-related functions. +namespace light {} + +} // namespace physics + +#include "physics/gas/atmosphere.hpp" +#include "physics/gas/ozone.hpp" + +#endif // ANTKEEPER_PHYSICS_GAS_HPP diff --git a/src/physics/gas/ozone.hpp b/src/physics/gas/ozone.hpp new file mode 100644 index 0000000..c8d17f4 --- /dev/null +++ b/src/physics/gas/ozone.hpp @@ -0,0 +1,75 @@ +/* + * 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_GAS_OZONE_HPP +#define ANTKEEPER_PHYSICS_GAS_OZONE_HPP + +namespace physics { +namespace gas { + +/// Ozone-related constants and functions. +namespace ozone { + +/** + * Cross section of ozone at wavelength 680nm and temperature 293K, in m-2/molecule. + * + * @see https://www.iup.uni-bremen.de/gruppen/molspec/databases/referencespectra/o3spectra2011/index.html + */ +template +constexpr T cross_section_680nm_293k = T(1.36820899679147e-25); + +/** + * Cross section of ozone at wavelength 550nm and temperature 293K, in m-2/molecule. + * + * @see https://www.iup.uni-bremen.de/gruppen/molspec/databases/referencespectra/o3spectra2011/index.html + */ +template +constexpr T cross_section_550nm_293k = T(3.31405330400124e-25); + +/** + * Cross section of ozone at wavelength 550nm and temperature 293K, in m-2/molecule. + * + * @see https://www.iup.uni-bremen.de/gruppen/molspec/databases/referencespectra/o3spectra2011/index.html + */ +template +constexpr T cross_section_440nm_293k = T(1.36017282525383e-26); + +/** + * Calculates an ozone absorption coefficient (wavelength-dependent). + * + * @param cross_section Ozone cross section of a wavelength, in m-2/molecule. + * @param n Molecular number density of standard air, in molecules/m-3. + * @param c Concentration of ozone in the atmosphere, unitless. + * + * @return Ozone absorption coefficient. + * + * @see https://skyrenderer.blogspot.com/2012/10/ozone-absorption.html + */ +template +T absorption(T cross_section, T n, T c) +{ + return cross_section * n * c; +} + +} // namespace ozone + +} // namespace gas +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_GAS_OZONE_HPP diff --git a/src/physics/number-density.hpp b/src/physics/number-density.hpp new file mode 100644 index 0000000..3eea24b --- /dev/null +++ b/src/physics/number-density.hpp @@ -0,0 +1,43 @@ +/* + * 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_NUMBER_DENSITY_HPP +#define ANTKEEPER_PHYSICS_NUMBER_DENSITY_HPP + +#include "physics/constants.hpp" + +namespace physics { + +/** + * Calculates the number density of a substance. + * + * @param c Molar concentration, in mol/m-3. + * @return Number density, in m-3. + * + * @see https://en.wikipedia.org/wiki/Number_density + */ +template +T number_density(T c) +{ + return physics::constants::avogadro * c; +} + +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_NUMBER_DENSITY_HPP diff --git a/src/physics/physics.hpp b/src/physics/physics.hpp index 1aa33d9..5d31e3c 100644 --- a/src/physics/physics.hpp +++ b/src/physics/physics.hpp @@ -25,8 +25,9 @@ namespace physics {} #include "constants.hpp" #include "frame.hpp" +#include "number-density.hpp" #include "planck.hpp" - +#include "gas/gas.hpp" #include "light/light.hpp" #include "orbit/orbit.hpp" #include "time/time.hpp" diff --git a/src/render/passes/sky-pass.cpp b/src/render/passes/sky-pass.cpp index 80873fb..254b347 100644 --- a/src/render/passes/sky-pass.cpp +++ b/src/render/passes/sky-pass.cpp @@ -160,14 +160,18 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const if (sun_illuminance_input) sun_illuminance_input->upload(sun_illuminance); - if (scale_height_rm_input) - scale_height_rm_input->upload(scale_height_rm); + if (distribution_rm_input) + distribution_rm_input->upload(distribution_rm); + if (distribution_o_input) + distribution_o_input->upload(distribution_o); if (rayleigh_scattering_input) rayleigh_scattering_input->upload(rayleigh_scattering); if (mie_scattering_input) mie_scattering_input->upload(mie_scattering); if (mie_anisotropy_input) mie_anisotropy_input->upload(mie_anisotropy); + if (ozone_absorption_input) + ozone_absorption_input->upload(ozone_absorption); if (atmosphere_radii_input) atmosphere_radii_input->upload(atmosphere_radii); @@ -299,10 +303,12 @@ void sky_pass::set_sky_model(const model* model) sun_luminance_input = sky_shader_program->get_input("sun_luminance"); sun_illuminance_input = sky_shader_program->get_input("sun_illuminance"); sun_angular_radius_input = sky_shader_program->get_input("sun_angular_radius"); - scale_height_rm_input = sky_shader_program->get_input("scale_height_rm"); + distribution_rm_input = sky_shader_program->get_input("distribution_rm"); + distribution_o_input = sky_shader_program->get_input("distribution_o"); rayleigh_scattering_input = sky_shader_program->get_input("rayleigh_scattering"); mie_scattering_input = sky_shader_program->get_input("mie_scattering"); mie_anisotropy_input = sky_shader_program->get_input("mie_anisotropy"); + ozone_absorption_input = sky_shader_program->get_input("ozone_absorption"); atmosphere_radii_input = sky_shader_program->get_input("atmosphere_radii"); } } @@ -480,9 +486,20 @@ void sky_pass::set_observer_altitude(float altitude) observer_altitude_tween[1] = altitude; } -void sky_pass::set_scale_heights(float rayleigh, float mie) +void sky_pass::set_particle_distributions(float rayleigh_scale_height, float mie_scale_height, float ozone_lower_limit, float ozone_upper_limit, float ozone_mode) { - scale_height_rm = {rayleigh, mie}; + distribution_rm = + { + -1.0f / rayleigh_scale_height, + -1.0f / mie_scale_height + }; + + distribution_o = + { + 1.0f / (ozone_lower_limit - ozone_mode), + 1.0f / (ozone_upper_limit - ozone_mode), + ozone_mode + }; } void sky_pass::set_scattering_coefficients(const float3& r, const float3& m) @@ -496,6 +513,11 @@ void sky_pass::set_mie_anisotropy(float g) mie_anisotropy = {g, g * g}; } +void sky_pass::set_absorption_coefficients(const float3& o) +{ + ozone_absorption = o; +} + void sky_pass::set_atmosphere_radii(float inner, float outer) { atmosphere_radii.x = inner; diff --git a/src/render/passes/sky-pass.hpp b/src/render/passes/sky-pass.hpp index 7c68277..7798554 100644 --- a/src/render/passes/sky-pass.hpp +++ b/src/render/passes/sky-pass.hpp @@ -69,10 +69,11 @@ public: void set_sun_illuminance(const float3& illuminance); void set_sun_angular_radius(float radius); void set_observer_altitude(float altitude); - void set_scale_heights(float rayleigh, float mie); void set_scattering_coefficients(const float3& r, const float3& m); void set_mie_anisotropy(float g); + void set_absorption_coefficients(const float3& o); void set_atmosphere_radii(float inner, float outer); + void set_particle_distributions(float rayleigh_scale_height, float mie_scale_height, float ozone_lower_limit, float ozone_upper_limit, float ozone_mode); void set_moon_position(const float3& position); void set_moon_rotation(const math::quaternion& rotation); @@ -96,11 +97,14 @@ private: const gl::shader_input* sun_luminance_input; const gl::shader_input* sun_illuminance_input; const gl::shader_input* sun_angular_radius_input; - const gl::shader_input* scale_height_rm_input; + const gl::shader_input* distribution_rm_input; + const gl::shader_input* distribution_o_input; const gl::shader_input* rayleigh_scattering_input; const gl::shader_input* mie_scattering_input; const gl::shader_input* mie_anisotropy_input; + const gl::shader_input* ozone_absorption_input; const gl::shader_input* atmosphere_radii_input; + gl::shader_program* moon_shader_program; const gl::shader_input* moon_model_input; @@ -173,11 +177,13 @@ private: tween moon_planetlight_illuminance_tween; float sun_angular_radius; - float2 scale_height_rm; float3 rayleigh_scattering; float3 mie_scattering; float2 mie_anisotropy; + float3 ozone_absorption; float3 atmosphere_radii; + float2 distribution_rm; + float3 distribution_o; float magnification; }; diff --git a/src/resources/entity-archetype-loader.cpp b/src/resources/entity-archetype-loader.cpp index 49054c3..2f71601 100644 --- a/src/resources/entity-archetype-loader.cpp +++ b/src/resources/entity-archetype-loader.cpp @@ -56,6 +56,16 @@ static bool load_component_atmosphere(entity::archetype& archetype, const json& component.mie_anisotropy = element["mie_anisotropy"].get(); if (element.contains("airglow")) component.airglow = element["airglow"].get(); + if (element.contains("air_concentration")) + component.air_concentration = element["air_concentration"].get(); + if (element.contains("ozone_concentration")) + component.ozone_concentration = element["ozone_concentration"].get(); + if (element.contains("ozone_lower_limit")) + component.ozone_lower_limit = element["ozone_lower_limit"].get(); + if (element.contains("ozone_upper_limit")) + component.ozone_upper_limit = element["ozone_upper_limit"].get(); + if (element.contains("ozone_mode")) + component.ozone_mode = element["ozone_mode"].get(); archetype.set(component);