Add ozone to atmosphere

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

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 <class T>
-math::vector3<T> transmittance(T depth_r, T depth_m, T depth_o, const math::vector3<T>& beta_r, const math::vector3<T>& beta_m)
+math::vector3<T> transmittance(T depth_r, T depth_m, T depth_o, const math::vector3<T>& beta_r, const math::vector3<T>& beta_m, const math::vector3<T>& beta_o)
 {
 	math::vector3<T> transmittance_r = beta_r * depth_r;
-	math::vector3<T> transmittance_m = beta_m * 1.1 * depth_m;
-	math::vector3<T> transmittance_o = {0, 0, 0};
+	math::vector3<T> transmittance_m = beta_m * (T(10)/T(9)) * depth_m;
+	math::vector3<T> transmittance_o = beta_o * depth_o;
 	
 	math::vector3<T> 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<entity::component::atmosphere>(reference_entity);
 
-			// Altitude of observer in meters	
+			// Altitude of observer in meters
 			geom::ray<double> 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<entity::component::atmosphere>(reference_entity);
 			
-			sky_pass->set_scale_heights(reference_atmosphere.rayleigh_scale_height, reference_atmosphere.mie_scale_height);
+			sky_pass->set_particle_distributions
+			(
+				static_cast<float>(reference_atmosphere.rayleigh_scale_height),
+				static_cast<float>(reference_atmosphere.mie_scale_height),
+				static_cast<float>(reference_atmosphere.ozone_lower_limit),
+				static_cast<float>(reference_atmosphere.ozone_upper_limit),
+				static_cast<float>(reference_atmosphere.ozone_mode)
+			);
 			sky_pass->set_scattering_coefficients(math::type_cast<float>(reference_atmosphere.rayleigh_scattering), math::type_cast<float>(reference_atmosphere.mie_scattering));
 			sky_pass->set_mie_anisotropy(reference_atmosphere.mie_anisotropy);
+			sky_pass->set_absorption_coefficients(math::type_cast<float>(reference_atmosphere.ozone_absorption));
 			sky_pass->set_atmosphere_radii(reference_body.radius, reference_body.radius + reference_atmosphere.exosphere_altitude);
 		}
 	}

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<entity::component::atmosphere>().connect<&atmosphere::on_atmosphere_construct>(this);
 	registry.on_replace<entity::component::atmosphere>().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<component::atmosphere>(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)

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

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

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<double> visible_wavelengths_nm;
 };
 

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 <cxxopts.hpp>
 #include <entt/entt.hpp>
 #include <filesystem>
@@ -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<float>(viewport_dimensions[0]), static_cast<float>(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<double>,
+		physics::gas::ozone::cross_section_550nm_293k<double>,
+		physics::gas::ozone::cross_section_440nm_293k<double>
+	};
+	
+	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);

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<float>());
@@ -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<float>());

src/physics/atmosphere.hpp → src/physics/gas/atmosphere.hpp

@@ -17,35 +17,20 @@
  * along with Antkeeper source code.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#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 <algorithm>
 #include <cmath>
 
 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 <class T>
-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<T> * density / (wavelength2 * wavelength2) * polarization;
+	return math::four_pi<T> * 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 <class T>
 T scattering_mie(T density, T polarization)
 {
-	return T(4) * math::pi<T> * density * polarization;
+	return math::four_pi<T> * 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 <class T>
+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 <class T>
-T optical_depth(const math::vector3<T>& a, const math::vector3<T>& b, T r, T sh, std::size_t n)
+T optical_depth_exp(const math::vector3<T>& a, const math::vector3<T>& 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 <class T>
+T optical_depth_tri(const math::vector3<T>& p0, const math::vector3<T>& 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<T> dy = (p1 - p0) / T(n);
+	math::vector3<T> 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 <class T>
+	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 <class T>
+	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

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 <http://www.gnu.org/licenses/>.
+ */
+
+#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

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 <http://www.gnu.org/licenses/>.
+ */
+
+#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 <class T>
+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 <class T>
+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 <class T>
+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 <class T>
+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

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 <http://www.gnu.org/licenses/>.
+ */
+
+#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 <class T>
+T number_density(T c)
+{
+	return physics::constants::avogadro<T> * c;
+}
+
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_NUMBER_DENSITY_HPP

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"

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;

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<float>& 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<float3> 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;
 };

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<double>();
 	if (element.contains("airglow"))
 		component.airglow = element["airglow"].get<double>();
+	if (element.contains("air_concentration"))
+		component.air_concentration = element["air_concentration"].get<double>();
+	if (element.contains("ozone_concentration"))
+		component.ozone_concentration = element["ozone_concentration"].get<double>();
+	if (element.contains("ozone_lower_limit"))
+		component.ozone_lower_limit = element["ozone_lower_limit"].get<double>();
+	if (element.contains("ozone_upper_limit"))
+		component.ozone_upper_limit = element["ozone_upper_limit"].get<double>();
+	if (element.contains("ozone_mode"))
+		component.ozone_mode = element["ozone_mode"].get<double>();
 
 	archetype.set<entity::component::atmosphere>(component);