Add more blackbody-related functions, add functions related to refraction, improve blackbody and atmosphere-related calculations in the astronomy system

src/ecs/components/atmosphere-component.hpp

@@ -45,8 +45,8 @@ struct atmosphere_component
 	/// Mie scale height, in meters.
 	double mie_scale_height;
 	
-	/// Mie phase function asymmetry factor.
-	double mie_asymmetry;
+	/// Mie phase function anisotropy factor.
+	double mie_anisotropy;
 	
 	/// (Dependent) Rayleigh scattering coefficients at sea level.
 	double3 rayleigh_scattering;

src/ecs/components/blackbody-component.hpp

@@ -31,17 +31,8 @@ struct blackbody_component
 	/// 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;
+	/// (Dependent) RGB luminous intensity, in candela.
+	double3 luminous_intensity;
 };
 
 } // namespace ecs

src/ecs/systems/astronomy-system.cpp

@@ -30,7 +30,9 @@
 #include "physics/light/blackbody.hpp"
 #include "physics/light/photometry.hpp"
 #include "physics/light/luminosity.hpp"
+#include "physics/light/refraction.hpp"
 #include "physics/atmosphere.hpp"
+#include "math/quadrature.hpp"
 #include "geom/cartesian.hpp"
 #include <iostream>
 
@@ -40,7 +42,7 @@ template
 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_r = beta_r * depth_r;
-	math::vector3<T> transmittance_m = beta_m * depth_m;
+	math::vector3<T> transmittance_m = beta_m * 1.1 * depth_m;
 	math::vector3<T> transmittance_o = {0, 0, 0};
 	
 	math::vector3<T> t = transmittance_r + transmittance_m + transmittance_o;
@@ -61,6 +63,10 @@ astronomy_system::astronomy_system(ecs::registry& registry):
 	sun_light(nullptr),
 	sky_pass(nullptr)
 {
+	// RGB wavelengths determined by matching wavelengths to XYZ, transforming XYZ to ACEScg, then selecting the max wavelengths for R, G, and B.
+	rgb_wavelengths_nm = {602.224, 541.069, 448.143};
+	rgb_wavelengths_m = rgb_wavelengths_nm * 1e-9;
+	
 	registry.on_construct<ecs::blackbody_component>().connect<&astronomy_system::on_blackbody_construct>(this);
 	registry.on_replace<ecs::blackbody_component>().connect<&astronomy_system::on_blackbody_replace>(this);
 	
@@ -128,13 +134,13 @@ void astronomy_system::update(double t, double dt)
 		// Calculate distance from observer to blackbody
 		double blackbody_distance = math::length(blackbody_position_topocentric);
 		
-		// Calculate blackbody illuminance according to distance
-		double blackbody_illuminance = blackbody.luminous_intensity / (blackbody_distance * blackbody_distance);
+		// Calculate blackbody distance attenuation
+		double distance_attenuation = 1.0 / (blackbody_distance * blackbody_distance);
 		
-		// Get blackbody color
-		double3 blackbody_color = blackbody.color;
+		// Init atmospheric transmittance
+		double3 atmospheric_transmittance = {1.0, 1.0, 1.0};
 		
-		// Get atmosphere component of reference body, if any
+		// Get atmosphere component of reference body (if any)
 		if (this->registry.has<ecs::atmosphere_component>(reference_body))
 		{
 			const ecs::atmosphere_component& atmosphere = this->registry.get<ecs::atmosphere_component>(reference_body);
@@ -159,10 +165,7 @@ void astronomy_system::update(double t, double dt)
 				double optical_depth_k = physics::atmosphere::optical_depth(sample_start, sample_end, earth_radius, atmosphere.mie_scale_height, 32);
 				double optical_depth_o = 0.0;
 				
-				double3 attenuation = transmittance(optical_depth_r, optical_depth_k, optical_depth_o, atmosphere.rayleigh_scattering, atmosphere.mie_scattering);
-				
-				// Attenuate blackbody color
-				blackbody_color *= attenuation;
+				atmospheric_transmittance = transmittance(optical_depth_r, optical_depth_k, optical_depth_o, atmosphere.rayleigh_scattering, atmosphere.mie_scattering);
 			}
 		}
 		
@@ -180,14 +183,14 @@ void astronomy_system::update(double t, double dt)
 			);
 			
 			// Update blackbody light color and intensity
-			sun_light->set_color(math::type_cast<float>(blackbody_color));	
-			sun_light->set_intensity(static_cast<float>(blackbody_illuminance));
+			sun_light->set_color(math::type_cast<float>(blackbody.luminous_intensity * atmospheric_transmittance));	
+			sun_light->set_intensity(static_cast<float>(distance_attenuation));
 			
 			// Upload blackbody params to sky pass
 			if (this->sky_pass)
 			{
 				this->sky_pass->set_sun_position(math::type_cast<float>(blackbody_position_topocentric));
-				this->sky_pass->set_sun_color(math::type_cast<float>(blackbody.color * blackbody_illuminance));
+				this->sky_pass->set_sun_color(math::type_cast<float>(blackbody.luminous_intensity * distance_attenuation));
 				
 				double blackbody_angular_radius = std::asin((blackbody.radius * 2.0) / (blackbody_distance * 2.0));
 				this->sky_pass->set_sun_angular_radius(static_cast<float>(blackbody_angular_radius));
@@ -219,7 +222,7 @@ void astronomy_system::update(double t, double dt)
 			
 			sky_pass->set_scale_heights(atmosphere.rayleigh_scale_height, atmosphere.mie_scale_height);
 			sky_pass->set_scattering_coefficients(math::type_cast<float>(atmosphere.rayleigh_scattering), math::type_cast<float>(atmosphere.mie_scattering));
-			sky_pass->set_mie_asymmetry(atmosphere.mie_asymmetry);
+			sky_pass->set_mie_anisotropy(atmosphere.mie_anisotropy);
 			sky_pass->set_atmosphere_radii(earth_radius, earth_radius + atmosphere.exosphere_altitude);
 		}
 	}
@@ -289,43 +292,31 @@ void astronomy_system::on_blackbody_construct(ecs::registry& registry, ecs::enti
 
 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<double> * blackbody.radius * blackbody.radius;
-	
-	// Calculate radiant flux
-	blackbody.radiant_flux = physics::light::blackbody::radiant_flux(blackbody.temperature, surface_area);
+	// Calculate the surface area of a spherical blackbody
+	const double surface_area = 4.0 * math::pi<double> * blackbody.radius * blackbody.radius;
 	
-	// Blackbody spectral power distribution function
-	auto spd = [blackbody](double x) -> double
+	// Construct a lambda function which calculates the blackbody's RGB luminous intensity of a given wavelength
+	auto rgb_luminous_intensity = [blackbody, surface_area](double wavelength_nm) -> double3
 	{
-		// Convert nanometers to meters
-		x *= double(1e-9);
+		// Convert wavelength from nanometers to meters
+		const double wavelength_m = wavelength_nm * 1e-9;
 		
-		return physics::light::blackbody::spectral_radiance<double>(blackbody.temperature, x, physics::constants::speed_of_light<double>);
-	};
-	
-	// Luminous efficiency function (photopic)
-	auto lef = [](double x) -> double
-	{
-		return physics::light::luminosity::photopic<double>(x);
+		// Calculate the spectral intensity of the wavelength
+		const double spectral_intensity = physics::light::blackbody::spectral_intensity<double>(blackbody.temperature, surface_area, wavelength_m);
+		
+		// Calculate the ACEScg color of the wavelength using CIE color matching functions
+		double3 spectral_color = color::xyz::to_acescg(color::xyz::match(wavelength_nm));
+		
+		// Scale the spectral color by spectral intensity
+		return spectral_color * spectral_intensity * 1e-9 * physics::light::max_luminous_efficacy<double>;
 	};
 	
-	// Construct range of spectral sample points
-	std::vector<double> samples(10000);
-	std::iota(samples.begin(), samples.end(), 10);
-	
-	// Calculate luminous efficiency
-	const double efficiency = physics::light::luminous_efficiency<double>(spd, lef, samples.begin(), samples.end());
+	// Construct a range of sample wavelengths in the visible spectrum
+	std::vector<double> samples(780 - 280);
+	std::iota(samples.begin(), samples.end(), 280);
 	
-	// Convert radiant flux to luminous flux
-	blackbody.luminous_flux = physics::light::watts_to_lumens<double>(blackbody.radiant_flux, efficiency);
-	
-	// Calculate luminous intensity from luminous flux
-	blackbody.luminous_intensity = blackbody.luminous_flux / (4.0 * math::pi<double>);
-	
-	// Calculate blackbody color from temperature
-	double3 color_xyz = color::cct::to_xyz(blackbody.temperature);
-	blackbody.color = color::xyz::to_acescg(color_xyz);
+	// Integrate the blackbody RGB luminous intensity over wavelengths in the visible spectrum
+	blackbody.luminous_intensity = math::quadrature::simpson(rgb_luminous_intensity, samples.begin(), samples.end());
 }
 
 void astronomy_system::on_atmosphere_construct(ecs::registry& registry, ecs::entity entity, ecs::atmosphere_component& atmosphere)
@@ -339,19 +330,16 @@ void astronomy_system::on_atmosphere_replace(ecs::registry& registry, ecs::entit
 	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);
 	
-	// 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 = {600.0e-9, 540.0e-9, 450.0e-9};
-	
 	// Calculate Rayleigh scattering coefficients
 	atmosphere.rayleigh_scattering =
 	{
-		physics::atmosphere::scatter_rayleigh(acescg_wavelengths.x, atmosphere.rayleigh_density, rayleigh_polarization),
-		physics::atmosphere::scatter_rayleigh(acescg_wavelengths.y, atmosphere.rayleigh_density, rayleigh_polarization),
-		physics::atmosphere::scatter_rayleigh(acescg_wavelengths.z, atmosphere.rayleigh_density, rayleigh_polarization)
+		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)
 	};
 	
 	// Calculate Mie scattering coefficients
-	const double mie_scattering = physics::atmosphere::scatter_mie(atmosphere.mie_density, mie_polarization);
+	const double mie_scattering = physics::atmosphere::scattering_mie(atmosphere.mie_density, mie_polarization);
 	atmosphere.mie_scattering = 
 	{
 		mie_scattering,

src/ecs/systems/astronomy-system.hpp

@@ -109,6 +109,9 @@ private:
 	physics::frame<double> inertial_to_topocentric;
 	physics::frame<double> sez_to_ezs;
 	physics::frame<double> ezs_to_sez;
+	
+	double3 rgb_wavelengths_nm;
+	double3 rgb_wavelengths_m;
 };
 
 } // namespace ecs

src/ecs/systems/tool-system.cpp

@@ -292,6 +292,9 @@ void tool_system::set_active_tool(ecs::entity entity)
 
 void tool_system::set_tool_active(bool active)
 {
+	if (active_tool == entt::null)
+		return;
+	
 	tool_active = active;
 	
 	if (active)

src/game/states/play-state.cpp

@@ -130,15 +130,15 @@ void play_state_enter(game_context* ctx)
 		orbit.elements.ta = math::radians(100.46457166) - longitude_periapsis;
 		
 		ecs::atmosphere_component atmosphere;
-		atmosphere.exosphere_altitude = 100e3;
+		atmosphere.exosphere_altitude = 65e3;
 		
-		atmosphere.index_of_refraction = 1.0003;
+		atmosphere.index_of_refraction = 1.000293;
 		atmosphere.rayleigh_density = 2.545e25;
 		atmosphere.mie_density = 14.8875;
 		
 		atmosphere.rayleigh_scale_height = 8000.0;
 		atmosphere.mie_scale_height = 1200.0;
-		atmosphere.mie_asymmetry = 0.8;
+		atmosphere.mie_anisotropy = 0.8;
 		
 		ecs::transform_component transform;
 		transform.local = math::identity_transform<float>;
@@ -230,7 +230,7 @@ void play_state_enter(game_context* ctx)
 	ecs::command::assign_render_layers(ecs_registry, ctx->twig_entity, 0);
 	
 	// Activate brush tool
-	ctx->tool_system->set_active_tool(ctx->brush_entity);
+	//ctx->tool_system->set_active_tool(ctx->brush_entity);
 
 	// Create ant-hill
 	auto ant_hill_entity = ant_hill_archetype->create(ecs_registry);
@@ -291,6 +291,7 @@ void play_state_enter(game_context* ctx)
 	
 	// Setup camera
 	ctx->overworld_camera->look_at({0, 0, 1}, {0, 0, 0}, {0, 1, 0});
+	ctx->overworld_camera->set_exposure(-14.5f);
 	ctx->camera_system->set_camera(ctx->overworld_camera);
 	
 	ctx->overworld_scene->update_tweens();

src/math/quadrature.hpp

@@ -63,7 +63,7 @@ typename std::invoke_result::val
 	typedef decltype(*last - *first) difference_type;
 	
 	if (first == last)
-		return output_type(0);
+		return output_type{0};
 	
 	output_type f_a = f(*first);
 	
@@ -102,7 +102,7 @@ typename std::invoke_result::val
 	typedef decltype(*last - *first) difference_type;
 	
 	if (first == last)
-		return output_type(0);
+		return output_type{0};
 	
 	output_type f_a = f(*first);
 	

src/physics/atmosphere.hpp

@@ -46,19 +46,20 @@ T density(T d0, T z, T sh)
 }
 
 /**
- * Calculates a particle polarization factor used in computing scattering coefficients.
+ * Calculates a particle polarizability factor used in computing scattering coefficients.
  *
  * @param ior Atmospheric index of refraction.
  * @param density Molecular density.
  * @return Polarization factor.
  *
  * @see Elek, Oskar. (2009). Rendering Parametrizable Planetary Atmospheres with Multiple Scattering in Real-Time.
+ * @see Real-Time Spectral Scattering in Large-Scale Natural Participating Media.
  */
 template <class T>
 T polarization(T ior, T density)
 {
-	const T ior2 = ior * ior;
-	const T num = T(2) * math::pi<T> * math::pi<T> * ((ior2 - T(1.0)) * (ior2 - T(1.0)));
+	const T ior2m1 = ior * ior - T(1.0);
+	const T num = T(2) * math::pi<T> * math::pi<T> * ior2m1 * ior2m1;
 	const T den = T(3) * density * density;
 	return num / den;
 }
@@ -73,9 +74,10 @@ T polarization(T ior, T density)
  * @see atmosphere::polarization
  *
  * @see Elek, Oskar. (2009). Rendering Parametrizable Planetary Atmospheres with Multiple Scattering in Real-Time.
+ * @see Real-Time Spectral Scattering in Large-Scale Natural Participating Media.
  */
 template <class T>
-T scatter_rayleigh(T wavelength, T density, T polarization)
+T scattering_rayleigh(T wavelength, T density, T polarization)
 {
 	const T wavelength2 = wavelength * wavelength;
 	return T(4) * math::pi<T> * density / (wavelength2 * wavelength2) * polarization;
@@ -90,24 +92,25 @@ T scatter_rayleigh(T wavelength, T density, T polarization)
  * @see atmosphere::polarization
  *
  * @see Elek, Oskar. (2009). Rendering Parametrizable Planetary Atmospheres with Multiple Scattering in Real-Time.
+ * @see Real-Time Spectral Scattering in Large-Scale Natural Participating Media.
  */
 template <class T>
-T scatter_mie(T density, T polarization)
+T scattering_mie(T density, T polarization)
 {
 	return T(4) * math::pi<T> * density * polarization;
 }
 
 /**
- * Calculates an extinction coefficient given a scattering coefficient and an absorbtion coefficient.
+ * Calculates attenuation due to extinction (absorption + out-scattering).
  *
- * @param s Scattering coefficient.
- * @param a Absorbtion coefficient.
- * @return Extinction coefficient.
+ * @param ec Extinction coefficient (absorption coefficient + scattering coefficient).
+ * @param s Scale factor.
+ * @return Attenuation factor.
  */
 template <class T>
-T extinction(T s, T a)
+T extinction(T ec, T s)
 {
-	return s + a;
+	return std::exp(-(ec * s));
 }
 
 /**

src/physics/light/blackbody.hpp

@@ -22,13 +22,14 @@
 
 #include "math/constants.hpp"
 #include "physics/constants.hpp"
-#include "physics/planck.hpp"
 
 namespace physics {
 namespace light {
 
 /**
  * Blackbody radiation functions.
+ *
+ * @see https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law
  */
 namespace blackbody {
 
@@ -37,9 +38,6 @@ namespace blackbody {
  *
  * @param t Temperature of the blackbody, in kelvin.
  * @return Radiant exitance of the blackbody, in watt per square meter.
- *
- * @see https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law
- * @see https://en.wikipedia.org/wiki/Radiant_exitance
  */
 template <class T>
 T radiant_exitance(T t);
@@ -49,25 +47,66 @@ T radiant_exitance(T t);
  *
  * @param t Temperature of the blackbody, in kelvin.
  * @param a Surface area of the blackbody, in square meters.
- * @return Radiant power of the blackbody, in watts.
- *
- * @see https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law
+ * @return Radiant flux of the blackbody, in watt.
  */
 template <class T>
 T radiant_flux(T t, T a);
 
 /**
- * Calculates the spectral radiance of a blackbody at a given wavelength.
+ * Calculates the radiant intensity of a blackbody.
+ *
+ * @param t Temperature of the blackbody, in kelvin.
+ * @param a Surface area of the blackbody, in square meters.
+ * @return Radiant intensity of the blackbody, in watt per steradian.
+ */
+template <class T>
+T radiant_intensity(T t, T a);
+
+/**
+ * Calculates the spectral exitance of a blackbody for the given wavelength.
  *
  * @param t Temperature of the blackbody, in kelvin.
  * @param lambda Wavelength of light, in meters.
  * @param c Speed of light in medium.
- * @return Spectral radiance, in watt per steradian per square meter per meter.
+ * @return Spectral exitance, in watt per square meter, per meter.
+ */
+template <class T>
+T spectral_exitance(T t, T lambda, T c = constants::speed_of_light<T>);
+
+/**
+ * Calculates the spectral flux of a blackbody for the given wavelength.
+ *
+ * @param t Temperature of the blackbody, in kelvin.
+ * @param a Surface area of the blackbody, in square meters.
+ * @param lambda Wavelength of light, in meters.
+ * @param c Speed of light in medium.
+ * @return Spectral flux of the blackbody, in watt per meter.
+ */
+template <class T>
+T spectral_flux(T t, T a, T lambda, T c  = constants::speed_of_light<T>);
+
+/**
+ * Calculates the spectral intensity of a blackbody for the given wavelength.
+ *
+ * @param t Temperature of the blackbody, in kelvin.
+ * @param a Surface area of the blackbody, in square meters.
+ * @param lambda Wavelength of light, in meters.
+ * @param c Speed of light in medium.
+ * @return Spectral intensity of the blackbody for the given wavelength, in watt per steradian per meter.
+ */
+template <class T>
+T spectral_intensity(T t, T a, T lambda, T c  = constants::speed_of_light<T>);
+
+/**
+ * Calculates the spectral radiance of a blackbody for the given wavelength.
  *
- * @see physics::planck::wavelength
+ * @param t Temperature of the blackbody, in kelvin.
+ * @param lambda Wavelength of light, in meters.
+ * @param c Speed of light in medium.
+ * @return Spectral radiance, in watt per steradian per square meter per meter.
  */
 template <class T>
-constexpr auto spectral_radiance = planck::wavelength<T>;
+T spectral_radiance(T t, T lambda, T c = constants::speed_of_light<T>);
 
 template <class T>
 T radiant_exitance(T t)
@@ -82,6 +121,54 @@ T radiant_flux(T t, T a)
 	return a * radiant_exitance(t);
 }
 
+template <class T>
+T radiant_intensity(T t, T a)
+{
+	return radiant_flux(t, a) / (T(4) * math::pi<T>);
+}
+
+template <class T>
+T spectral_exitance(T t, T lambda, T c)
+{
+	const T hc = constants::planck<T> * c;
+	const T lambda2 = lambda * lambda;
+	
+	// First radiation constant (c1)
+	const T c1 = T(2) * math::pi<T> * hc * c;
+	
+	// Second radiation constant (c2)
+	const T c2 = hc / constants::boltzmann<T>;
+	
+	return (c1 / (lambda2 * lambda2 * lambda)) / std::expm1(c2 / (lambda * t));
+}
+
+template <class T>
+T spectral_flux(T t, T a, T lambda, T c)
+{
+	return a * spectral_exitance(t, lambda, c);
+}
+
+template <class T>
+T spectral_intensity(T t, T a, T lambda, T c)
+{
+	return spectral_flux(t, a, lambda, c) / (T(4) * math::pi<T>);
+}
+
+template <class T>
+T spectral_radiance(T t, T lambda, T c)
+{
+	const T hc = constants::planck<T> * c;
+	const T lambda2 = lambda * lambda;
+	
+	// First radiation constant (c1L)
+	const T c1l = T(2) * hc * c;
+	
+	// Second radiation constant (c2)
+	const T c2 = hc / constants::boltzmann<T>;
+	
+	return (c1l / (lambda2 * lambda2 * lambda)) / std::expm1(c2 / (lambda * t));
+}
+
 } // namespace blackbody
 
 } // namespace light

src/physics/light/light.hpp

@@ -31,5 +31,6 @@ namespace light {}
 #include "luminosity.hpp"
 #include "phase.hpp"
 #include "photometry.hpp"
+#include "refraction.hpp"
 
 #endif // ANTKEEPER_PHYSICS_LIGHT_HPP

src/physics/light/phase.hpp

@@ -69,7 +69,7 @@ T rayleigh(T mu);
 template <class T>
 T cornette_shanks(T mu, T g)
 {
-	const T k = T(3) / (T(8) * math::pi<T>);
+	static const T k = T(3) / (T(8) * math::pi<T>);
 	const T gg = g * g;
 	const T num = (T(1) - gg) * (T(1) + mu * mu);
 	const T den = (T(2) + gg) * std::pow(T(1) + gg - T(2) * g * mu, T(1.5));
@@ -86,7 +86,8 @@ T henyey_greenstein(T mu, T g)
 template <class T>
 T rayleigh(T mu)
 {
-	return (T(1) + mu * mu) * (T(3) / (T(16) * math::pi<T>));
+	static const T k = T(3) / (T(16) * math::pi<T>);
+	return k * (1.0 + mu * mu);
 }
 
 } // namespace phase

src/physics/light/photometry.hpp

@@ -26,6 +26,10 @@
 namespace physics {
 namespace light {
 
+/// Maximum luminous efficacy of an ideal monochromatic source, in lumen per watt.
+template <class T>
+constexpr T max_luminous_efficacy = T(683.002);
+
 /**
  * Calculates the luminous efficiency of a light source.
  *
@@ -60,20 +64,20 @@ T luminous_efficiency(UnaryOp1 spd, UnaryOp2 lef, InputIt first, InputIt last)
 template <class T>
 T luminous_efficacy(T efficiency)
 {
-	return T(683.002) * efficiency;
+	return max_luminous_efficacy<T> * efficiency;
 }
 
 /**
  * Converts watts (radiant flux) to lumens (luminous flux).
  *
- * @param radient_flux Radiance flux, in watts.
+ * @param radiant_flux Radiant flux, in watts.
  * @param efficiency Luminous efficiency, on `[0, 1]`.
  * @return Luminous flux, in lumens.
  */
 template <class T>
-T watts_to_lumens(T radient_flux, T efficiency)
+T watts_to_lumens(T radiant_flux, T efficiency)
 {
-	return radient_flux * luminous_efficacy(efficiency);
+	return radiant_flux * luminous_efficacy(efficiency);
 }
 
 } // namespace light

src/physics/light/refraction.hpp

@@ -0,0 +1,117 @@
+/*
+ * 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_LIGHT_REFRACTION_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_REFRACTION_HPP
+
+namespace physics {
+namespace light {
+
+/// Index of refraction formulas.
+namespace ior {
+
+/**
+ * Two-term form of Cauchy's transmission equation.
+ *
+ * @param lambda Wavelength of light, in micrometers.
+ * @param a First coefficient of the medium.
+ * @param b Second coefficient of the medium.
+ * @return Refractive index.
+ *
+ * @see https://en.wikipedia.org/wiki/Cauchy%27s_equation
+ */
+template <class T>
+T cauchy(T lambda, T a, T b)
+{
+	return a + b / (lambda * lambda);
+}
+
+/**
+ * Three-term form of Cauchy's transmission equation.
+ *
+ * @param lambda Wavelength of light, in micrometers.
+ * @param a First coefficient.
+ * @param b Second coefficient.
+ * @param c Third coefficient.
+ * @return Refractive index.
+ *
+ * @see https://en.wikipedia.org/wiki/Cauchy%27s_equation
+ */
+template <class T>
+T cauchy(T lambda, T a, T b, T c)
+{
+	const T lambda2 = lambda * lambda;
+	return a + b / lambda2 + c / (lambda2 * lambda2);
+}
+
+/**
+ * Two-term form of the Sellmeier equation.
+ *
+ * @param lambda Wavelength of light, in micrometers.
+ * @param b1 B1 coefficient.
+ * @param c1 C1 coefficient.
+ * @param b2 B2 coefficient.
+ * @param c2 C2 coefficient.
+ * @return Refractive index.
+ *
+ * @see https://en.wikipedia.org/wiki/Sellmeier_equation
+ */
+template <class T>
+T sellmeier(T lambda, T b1, T c1, T b2, T c2)
+{
+	const T lambda2 = lambda * lambda;
+
+	const T t1 = (b1 * lambda2) / (lambda2 - c1);
+	const T t2 = (b2 * lambda2) / (lambda2 - c2);
+	
+	return std::sqrt(T(1) + t1 + t2);
+}
+
+/**
+ * Three-term form of the Sellmeier equation.
+ *
+ * @param lambda Wavelength of light, in micrometers.
+ * @param b1 B1 coefficient.
+ * @param c1 C1 coefficient.
+ * @param b2 B2 coefficient.
+ * @param c2 C2 coefficient.
+ * @param b3 B3 coefficient.
+ * @param c3 C3 coefficient.
+ * @return Refractive index.
+ *
+ * @see https://en.wikipedia.org/wiki/Sellmeier_equation
+ */
+template <class T>
+T sellmeier(T lambda, T b1, T c1, T b2, T c2, T b3, T c3)
+{
+	const T lambda2 = lambda * lambda;
+
+	const T t1 = (b1 * lambda2) / (lambda2 - c1);
+	const T t2 = (b2 * lambda2) / (lambda2 - c2);
+	const T t3 = (b3 * lambda2) / (lambda2 - c3);
+	
+	return std::sqrt(T(1) + t1 + t2 + t3);
+}
+
+} // namespace ior
+
+} // namespace light
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_REFRACTION_HPP

src/renderer/passes/shadow-map-pass.cpp

@@ -103,7 +103,8 @@ void shadow_map_pass::render(render_context* context) const
 	glDepthMask(GL_TRUE);
 	
 	// Disable face culling
-	glDisable(GL_CULL_FACE);
+	glEnable(GL_CULL_FACE);
+	glCullFace(GL_FRONT);
 	
 	// For half-z buffer
 	//glDepthRange(-1.0f, 1.0f);

src/renderer/passes/sky-pass.cpp

@@ -122,14 +122,14 @@ sky_pass::sky_pass(gl::rasterizer* rasterizer, const gl::framebuffer* framebuffe
 		// Transform XYZ color to ACEScg colorspace
 		double3 color_acescg = color::xyz::to_acescg(color_xyz);
 		
-		// Convert apparent magnitude to irradiance W/m2
-		double vmag_lux = astro::vmag_to_lux(vmag);
+		// Convert apparent magnitude to irradiance (W/m^2)
+		double vmag_irradiance = std::pow(10.0, 0.4 * (-vmag - 19.0 + 0.4));
 		
-		// Convert irradiance to illuminance (using luminous efficiency of sun)
-		double illuminance = physics::light::watts_to_lumens<double>(vmag_lux, 0.13);
+		// Convert irradiance to illuminance
+		double vmag_illuminance = vmag_irradiance * (683.0 * 0.14);
 		
 		// Scale color by illuminance
-		double3 scaled_color = color_acescg * illuminance;
+		double3 scaled_color = color_acescg * vmag_illuminance;
 		
 		// Build vertex
 		*(star_vertex++) = static_cast<float>(position_inertial.x);
@@ -244,8 +244,8 @@ void sky_pass::render(render_context* context) const
 			rayleigh_scattering_input->upload(rayleigh_scattering);
 		if (mie_scattering_input)
 			mie_scattering_input->upload(mie_scattering);
-		if (mie_asymmetry_input)
-			mie_asymmetry_input->upload(mie_asymmetry);
+		if (mie_anisotropy_input)
+			mie_anisotropy_input->upload(mie_anisotropy);
 		if (atmosphere_radii_input)
 			atmosphere_radii_input->upload(atmosphere_radii);
 		
@@ -351,7 +351,7 @@ void sky_pass::set_sky_model(const model* model)
 				scale_height_rm_input = sky_shader_program->get_input("scale_height_rm");
 				rayleigh_scattering_input = sky_shader_program->get_input("rayleigh_scattering");
 				mie_scattering_input = sky_shader_program->get_input("mie_scattering");
-				mie_asymmetry_input = sky_shader_program->get_input("mie_asymmetry");
+				mie_anisotropy_input = sky_shader_program->get_input("mie_anisotropy");
 				atmosphere_radii_input = sky_shader_program->get_input("atmosphere_radii");
 			}
 		}
@@ -449,9 +449,9 @@ void sky_pass::set_scattering_coefficients(const float3& r, const float3& m)
 	mie_scattering = m;
 }
 
-void sky_pass::set_mie_asymmetry(float g)
+void sky_pass::set_mie_anisotropy(float g)
 {
-	mie_asymmetry = {g, g * g};
+	mie_anisotropy = {g, g * g};
 }
 
 void sky_pass::set_atmosphere_radii(float inner, float outer)

src/renderer/passes/sky-pass.hpp

@@ -64,7 +64,7 @@ public:
 	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_asymmetry(float g);
+	void set_mie_anisotropy(float g);
 	void set_atmosphere_radii(float inner, float outer);
 
 private:
@@ -84,7 +84,7 @@ private:
 	const gl::shader_input* scale_height_rm_input;
 	const gl::shader_input* rayleigh_scattering_input;
 	const gl::shader_input* mie_scattering_input;
-	const gl::shader_input* mie_asymmetry_input;
+	const gl::shader_input* mie_anisotropy_input;
 	const gl::shader_input* atmosphere_radii_input;
 	
 	gl::shader_program* moon_shader_program;
@@ -132,7 +132,7 @@ private:
 	float2 scale_height_rm;
 	float3 rayleigh_scattering;
 	float3 mie_scattering;
-	float2 mie_asymmetry;
+	float2 mie_anisotropy;
 	float3 atmosphere_radii;
 };