Add photometric and radiometric functions to the physics namespace, add more blackbody functions, add quadrature namespace with trapezoid rule integral approximation function, make astronomy system capable of calculating the luminous flux of a blackbody

src/color/xyz.hpp

@@ -24,7 +24,11 @@
 
 namespace color {
 
-/// Functions which operate in the CIE XYZ colorspace.
+/**
+ * Functions which operate in the CIE XYZ colorspace.
+ *
+ * @see https://en.wikipedia.org/wiki/CIE_1931_color_space
+ */
 namespace xyz {
 
 /**
@@ -36,6 +40,54 @@ namespace xyz {
 template <class T>
 T luminance(const math::vector3<T>& x);
 
+/**
+ * Fitted piecewise gaussian approximation to the CIE 1931 standard observer color matching function.
+ *
+ * @param lambda Wavelength of light, in nanometers.
+ * @return Matching CIE XYZ color.
+ *
+ * @see match_x(T)
+ * @see match_y(T)
+ * @see match_z(T)
+ *
+ * @see Wyman, C., Sloan, P.J., & Shirley, P. (2013). Simple Analytic Approximations to the CIE XYZ Color Matching Functions.
+ */
+template <class T>
+math::vector3<T> match(T lambda);
+
+/**
+ * CIE 1931 standard observer color matching function for the X tristimulus value.
+ *
+ * @param lambda Wavelength of light, in nanometers.
+ * @return Matching X tristimulus value.
+ *
+ * @see match(T)
+ */
+template <class T>
+T match_x(T lambda);
+
+/**
+ * CIE 1931 standard observer color matching function for the Y tristimulus value.
+ *
+ * @param lambda Wavelength of light, in nanometers.
+ * @return Matching Y tristimulus value.
+ *
+ * @see match(T)
+ */
+template <class T>
+T match_y(T lambda);
+
+/**
+ * CIE 1931 standard observer color matching function for the Z tristimulus value.
+ *
+ * @param lambda Wavelength of light, in nanometers.
+ * @return Matching Z tristimulus value.
+ *
+ * @see match(T)
+ */
+template <class T>
+T match_z(T lambda);
+
 /**
  * Transforms a CIE XYZ color into the ACEScg colorspace.
  *
@@ -92,6 +144,55 @@ inline T luminance(const math::vector3& x)
 	return x[1];
 }
 
+template <class T>
+math::vector3<T> match(T lambda)
+{
+	return math::vector3<T>
+	{
+		match_x<T>(lambda),
+		match_y<T>(lambda),
+		match_z<T>(lambda)
+	};
+}
+
+template <class T>
+T match_x(T lambda)
+{
+	const T t0 = (lambda - T(442.0)) * ((lambda < T(442.0)) ? T(0.0624) : T(0.0374));
+	const T t1 = (lambda - T(599.8)) * ((lambda < T(599.8)) ? T(0.0264) : T(0.0323));
+	const T t2 = (lambda - T(501.1)) * ((lambda < T(501.1)) ? T(0.0490) : T(0.0382));
+	
+	const T x0 = T( 0.362) * std::exp(T(-0.5) * t0 * t0);
+	const T x1 = T( 1.056) * std::exp(T(-0.5) * t1 * t1);
+	const T x2 = T(-0.056) * std::exp(T(-0.5) * t2 * t2);
+	
+	return x0 + x1 + x2;
+}
+
+template <class T>
+T match_y(T lambda)
+{
+	const T t0 = (lambda - T(568.8)) * ((lambda < T(568.8)) ? T(0.0213) : T(0.0247));
+	const T t1 = (lambda - T(530.9)) * ((lambda < T(530.9)) ? T(0.0613) : T(0.0322));
+	
+	const T y0 = T(0.821) * std::exp(T(-0.5) * t0 * t0);
+	const T y1 = T(0.286) * std::exp(T(-0.5) * t1 * t1);
+	
+	return y0 + y1;
+}
+
+template <class T>
+T match_z(T lambda)
+{
+	const T t0 = (lambda - T(437.0)) * ((lambda < T(437.0)) ? T(0.0845) : T(0.0278));
+	const T t1 = (lambda - T(459.0)) * ((lambda < T(459.0)) ? T(0.0385) : T(0.0725));
+	
+	const T z0 = T(1.217) * std::exp(T(-0.5) * t0 * t0);
+	const T z1 = T(0.681) * std::exp(T(-0.5) * t1 * t1);
+	
+	return z0 + z1;
+}
+
 template <class T>
 math::vector3<T> to_acescg(const math::vector3<T>& x)
 {

src/ecs/components/atmosphere-component.hpp

@@ -20,16 +20,27 @@
 #ifndef ANTKEEPER_ECS_ATMOSPHERE_COMPONENT_HPP
 #define ANTKEEPER_ECS_ATMOSPHERE_COMPONENT_HPP
 
+#include "utility/fundamental-types.hpp"
+
 namespace ecs {
 
 /// Atmosphere
 struct atmosphere_component
 {
+	/// Radius of the outer atmosphere, in meters.
+	double exosphere_radius;
+	
 	/// Rayleigh scale height
-	double scale_rayleigh;
+	double rayleigh_scale_height;
 	
 	/// Mie scale height
-	double scale_mie;
+	double mie_scale_height;
+	
+	/// Rayleigh scattering coefficients
+	double3 rayleigh_scattering_coefficients;
+	
+	/// Mie scattering coefficients
+	double3 mie_scattering_coefficients;
 };
 
 } // namespace ecs

src/ecs/systems/astronomy-system.cpp

@@ -23,14 +23,60 @@
 #include "ecs/components/blackbody-component.hpp"
 #include "ecs/components/atmosphere-component.hpp"
 #include "ecs/components/transform-component.hpp"
+#include "geom/intersection.hpp"
 #include "color/color.hpp"
 #include "physics/orbit/orbit.hpp"
 #include "physics/time/ut1.hpp"
+#include "physics/light/blackbody.hpp"
+#include "physics/light/photometry.hpp"
+#include "physics/light/luminosity.hpp"
 #include "geom/cartesian.hpp"
 #include <iostream>
 
 namespace ecs {
 
+/**
+ * Calculates the optical depth between two points.
+ *
+ * @param start Start point.
+ * @param end End point.
+ * @param sample_count Number of samples to take between the start and end points.
+ * @param scale_height Scale height of the atmospheric particles to measure.
+ */
+template <class T>
+T transmittance(math::vector3<T> start, math::vector3<T> end, std::size_t sample_count, T scale_height)
+{
+	const T inverse_scale_height = T(1) / -scale_height;
+	
+	math::vector3<T> direction = end - start;
+	T distance = length(direction);
+	direction /= distance;
+	
+	// Calculate the distance between each sample point
+	T step_distance = distance / T(sample_count);
+	
+	// Sum the atmospheric particle densities at each sample point
+	T total_density = 0.0;
+	math::vector3<T> sample_point = start;
+	for (std::size_t i = 0; i < sample_count; ++i)
+	{
+		// Determine altitude of sample point
+		T altitude = length(sample_point);
+		
+		// Calculate atmospheric particle density at sample altitude
+		T density = exp(altitude * inverse_scale_height);
+		
+		// Add density to the total density
+		total_density += density;
+		
+		// Advance to next sample point
+		sample_point += direction * step_distance;
+	}
+	
+	// Scale the total density by the step distance
+	return total_density * step_distance;
+}
+
 astronomy_system::astronomy_system(ecs::registry& registry):
 	entity_system(registry),
 	universal_time(0.0),
@@ -118,16 +164,75 @@ void astronomy_system::update(double t, double dt)
 		
 		//std::cout << "el: " << math::degrees(sun_az_el.y) << "; az: " << math::degrees(sun_az_el.z) << std::endl;
 		
-		// Calculate sun color
-		float cct = 3000.0f + std::sin(sun_az_el.y) * 5000.0f;
-		float3 color_xyz = color::cct::to_xyz(cct);
-		float3 color_acescg = color::xyz::to_acescg(color_xyz);
-		sun_light->set_color(color_acescg);
-
-		// Calculate sun intensity (in lux)
-		const float illuminance_zenith = 108000.0f;
-		float illuminance = std::max(0.0, std::sin(sun_az_el.y) * illuminance_zenith);
-		sun_light->set_intensity(illuminance);
+		// If the reference body has an atmosphere
+		if (registry.has<ecs::atmosphere_component>(reference_body))
+		{
+			// Get the atmosphere component of the reference body
+			const auto& atmosphere = registry.get<ecs::atmosphere_component>(reference_body);
+			
+			const double meters_per_au = 1.496e+11;
+			const double earth_radius_au = 4.26352e-5;
+			const double earth_radius_m = earth_radius_au * meters_per_au;
+			const double observer_altitude_m = (observer_location[0] - earth_radius_au) * meters_per_au;
+			
+			// Altitude of observer in meters			
+			geom::ray<double> sample_ray;
+			sample_ray.origin = {0, observer_altitude_m, 0};
+			sample_ray.direction = math::normalize(sun_position_topocentric);
+			
+			geom::sphere<double> exosphere;
+			exosphere.center = {0, -earth_radius_m, 0};
+			exosphere.radius = atmosphere.exosphere_radius;
+			
+			auto intersection_result = geom::ray_sphere_intersection(sample_ray, exosphere);
+			
+			if (std::get<0>(intersection_result))
+			{
+				double3 sample_start = sample_ray.origin;
+				double3 sample_end = sample_ray.extrapolate(std::get<2>(intersection_result));
+				
+				double transmittance_rayleigh = transmittance(sample_start, sample_end, 16, atmosphere.rayleigh_scale_height);
+				double transmittance_mie = transmittance(sample_start, sample_end, 16, atmosphere.mie_scale_height);
+				
+				// Calculate attenuation due to atmospheric scattering
+				double3 scattering_attenuation =
+					atmosphere.rayleigh_scattering_coefficients * transmittance_rayleigh +
+					atmosphere.mie_scattering_coefficients * transmittance_mie;
+				scattering_attenuation.x = std::exp(-scattering_attenuation.x);
+				scattering_attenuation.y = std::exp(-scattering_attenuation.y);
+				scattering_attenuation.z = std::exp(-scattering_attenuation.z);
+				
+				double scattering_mean = (scattering_attenuation.x + scattering_attenuation.y + scattering_attenuation.z) / 3.0;
+				
+				const double sun_temperature = 5777.0;
+				const double sun_radius = 6.957e+8;
+				const double sun_surface_area = 4.0 * math::pi<double> * sun_radius * sun_radius;
+				
+				// Calculate distance attenuation
+				double sun_distance_m = math::length(sun_position_topocentric) * meters_per_au;
+				double distance_attenuation = 1.0 / (sun_distance_m * sun_distance_m);
+				
+				double sun_luminous_flux = blackbody_luminous_flux(sun_temperature, sun_radius);
+				double sun_luminous_intensity = sun_luminous_flux / (4.0 * math::pi<double>);
+				double sun_illuminance = sun_luminous_intensity / (sun_distance_m * sun_distance_m);
+				
+				std::cout << "distance atten: " << distance_attenuation << std::endl;
+				std::cout << "scatter atten: " << scattering_attenuation << std::endl;
+				
+				std::cout << "luminous flux: " << sun_luminous_flux << std::endl;
+				std::cout << "luminous intensity: " << sun_luminous_intensity << std::endl;
+				std::cout << "illuminance: " << sun_illuminance * scattering_mean << std::endl;
+				
+				
+				// Calculate sun color
+				double3 color_xyz = color::cct::to_xyz(sun_temperature);
+				double3 color_acescg = color::xyz::to_acescg(color_xyz);
+				
+				sun_light->set_color(math::type_cast<float>(color_acescg * scattering_attenuation));
+				
+				sun_light->set_intensity(sun_illuminance);
+			}
+		}
 	}
 	
 	if (sky_pass != nullptr)
@@ -202,4 +307,40 @@ void astronomy_system::set_sky_pass(::sky_pass* pass)
 	this->sky_pass = pass;
 }
 
+double astronomy_system::blackbody_luminous_flux(double t, double r)
+{
+	// Blackbody spectral power distribution function
+	auto spd = [t](double x) -> double
+	{
+		// Convert nanometers to meters
+		x *= double(1e-9);
+		
+		return physics::light::blackbody::spectral_radiance<double>(t, x, physics::constants::speed_of_light<double>);
+	};
+	
+	// Luminous efficiency function (photopic)
+	auto lef = [](double x) -> double
+	{
+		return physics::light::luminosity::photopic<double>(x);
+	};
+	
+	// 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());
+	
+	// Calculate surface area of spherical blackbody
+	const double a = double(4) * math::pi<double> * r * r;
+	
+	// Calculate radiant flux
+	const double radiant_flux = physics::light::blackbody::radiant_flux(t, a);
+	
+	// Convert radiant flux to luminous flux
+	const double luminous_flux = physics::light::watts_to_lumens<double>(radiant_flux, efficiency);
+	
+	return luminous_flux;
+}
+
 } // namespace ecs

src/ecs/systems/astronomy-system.hpp

@@ -86,6 +86,15 @@ public:
 	void set_sky_pass(sky_pass* pass);
 	
 private:
+	/**
+	 * Calculates the luminous flux of a spherical blackbody in vacuum using the CIE 1931 standard observer photopic luminosity function.
+	 *
+	 * @param t Temperature of the blackbody, in kelvin.
+	 * @param r Radius of the blackbody, in meters.
+	 * @return Luminous flux, in lumens.
+	 */
+	static double blackbody_luminous_flux(double t, double r);
+	
 	double universal_time;
 	double time_scale;
 	ecs::entity reference_body;

src/game/states/play-state.cpp

@@ -33,6 +33,7 @@
 #include "ecs/components/camera-follow-component.hpp"
 #include "ecs/components/orbit-component.hpp"
 #include "ecs/components/blackbody-component.hpp"
+#include "ecs/components/atmosphere-component.hpp"
 #include "ecs/components/light-component.hpp"
 #include "ecs/commands.hpp"
 #include "game/game-context.hpp"
@@ -136,11 +137,20 @@ void play_state_enter(game_context* ctx)
 		orbit.elements.w = longitude_periapsis - orbit.elements.raan;
 		orbit.elements.ta = math::radians(100.46457166) - longitude_periapsis;
 		
+		const double earth_radius_m = 6378e3;
+		ecs::atmosphere_component atmosphere;
+		atmosphere.exosphere_radius = earth_radius_m + 100e3;
+		atmosphere.rayleigh_scale_height = 8000.0;
+		atmosphere.mie_scale_height = 1200.0;
+		atmosphere.rayleigh_scattering_coefficients = {5.8e-6, 1.35e-5, 3.31e-5};
+		atmosphere.mie_scattering_coefficients = {2e-6, 2e-6, 2e-6};
+		
 		ecs::transform_component transform;
 		transform.local = math::identity_transform<float>;
 		transform.warp = true;
 		
 		ecs_registry.assign<ecs::orbit_component>(earth_entity, orbit);
+		ecs_registry.assign<ecs::atmosphere_component>(earth_entity, atmosphere);
 		ecs_registry.assign<ecs::transform_component>(earth_entity, transform);
 	}
 	

src/geom/intersection.hpp

@@ -24,6 +24,7 @@
 #include "geom/mesh.hpp"
 #include "geom/plane.hpp"
 #include "geom/ray.hpp"
+#include "geom/sphere.hpp"
 #include "utility/fundamental-types.hpp"
 #include <tuple>
 
@@ -44,6 +45,25 @@ std::tuple ray_aabb_intersection(const ray& ray, cons
 
 std::tuple<bool, float, float, std::size_t, std::size_t> ray_mesh_intersection(const ray<float>& ray, const mesh& mesh);
 
+/**
+ * Ray-sphere intersection test.
+ */
+template <class T>
+std::tuple<bool, T, T> ray_sphere_intersection(const ray<T>& ray, const sphere<T>& sphere)
+{
+    const auto d = ray.origin - sphere.center;
+    const T b = math::dot(d, ray.direction);
+	const T c = math::dot(d, d) - sphere.radius * sphere.radius;
+	T h = b * b - c;
+	
+	if (h < T(0))
+		return {false, T(0), T(0)};
+	
+	h = std::sqrt(h);
+	
+	return {true, -b - h, -b + h};
+}
+
 bool aabb_sphere_intersection(const aabb<float>& aabb, const float3& center, float radius);
 
 bool aabb_aabb_intersection(const aabb<float>& a, const aabb<float>& b);

src/math/math.hpp

@@ -43,6 +43,7 @@ namespace math {}
 
 #include "math/angles.hpp"
 #include "math/constants.hpp"
+#include "math/quadrature.hpp"
 #include "math/interpolation.hpp"
 #include "math/random.hpp"
 

src/math/quadrature.hpp

@@ -0,0 +1,73 @@
+/*
+ * 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_MATH_QUADRATURE_HPP
+#define ANTKEEPER_MATH_QUADRATURE_HPP
+
+#include <iterator>
+#include <type_traits>
+
+namespace math {
+
+/// Numerical integration functions.
+namespace quadrature {
+
+/**
+ * Approximates the definite integral of a function using the trapezoidal rule.
+ *
+ * @param f Unary function object to integrate.
+ * @param first,last Range of sample points on `[first, last)`.
+ * @return Approximated integral of @p f.
+ *
+ * @see https://en.wikipedia.org/wiki/Trapezoidal_rule
+ */
+template<class UnaryOp, class InputIt>
+typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type trapezoid(UnaryOp f, InputIt first, InputIt last)
+{
+	typedef typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type result_type;
+	
+	if (first == last)
+		return result_type(0);
+	
+	result_type f_a = f(*first);
+	
+	InputIt second = first;
+	++second;
+	
+	if (second == last)
+		return f_a;
+	
+	result_type f_b = f(*second);
+	result_type sum = (f_a + f_b) / result_type(2) * (*second - *first);
+	
+	for (first = second++; second != last; first = second++)
+	{
+		f_a = f_b;
+		f_b = f(*second);
+		sum += (f_a + f_b) / result_type(2) * (*second - *first);
+	}
+	
+	return sum;
+}
+
+} // namespace quadrature
+
+} // namespace math
+
+#endif // ANTKEEPER_MATH_QUADRATURE_HPP

src/physics/constants.hpp

@@ -0,0 +1,72 @@
+/*
+ * 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_CONSTANTS_HPP
+#define ANTKEEPER_PHYSICS_CONSTANTS_HPP
+
+namespace physics {
+
+/// Physical constants
+namespace constants {
+
+/**
+ * Boltzmann constant (k), in joule per kelvin.
+ *
+ * @see https://physics.nist.gov/cgi-bin/cuu/Value?k|search_for=universal_in!
+ */
+template <class T>
+constexpr T boltzmann = T(1.380649e-23);
+
+/**
+ * Gravitational constant (G), in cubic meter per second squared per kilogram.
+ *
+ * @see https://physics.nist.gov/cgi-bin/cuu/Value?G|search_for=universal_in!
+ */
+template <class T>
+constexpr T gravitational = T(6.67430e-11);
+
+/**
+ * Planck constant (h), in joule per hertz.
+ *
+ * @see https://physics.nist.gov/cgi-bin/cuu/Value?h|search_for=universal_in!
+ */
+template <class T>
+constexpr T planck = T(6.62607015e-34);
+
+/**
+ * Stefan-Boltzmann constant (sigma), in watt per square meter kelvin to the fourth.
+ *
+ * @see https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_constant
+ */
+template <class T>
+constexpr T stefan_boltzmann = T(5.67037441918442945397099673188923087584012297029130e-8);
+
+/**
+ * Speed of light in vacuum (c), in meters per second.
+ *
+ * @see https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=universal_in!
+ */
+template <class T>
+constexpr T speed_of_light = T(299792458);
+
+} // namespace constant
+
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_CONSTANTS_HPP

src/physics/light/blackbody.hpp

@@ -0,0 +1,90 @@
+/*
+ * 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_BLACKBODY_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_BLACKBODY_HPP
+
+#include "math/constants.hpp"
+#include "physics/constants.hpp"
+#include "physics/planck.hpp"
+
+namespace physics {
+namespace light {
+
+/**
+ * Blackbody radiation functions.
+ */
+namespace blackbody {
+
+/**
+ * Calculates the radiant exitance of a 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);
+
+/**
+ * Calculates the radiant flux of a blackbody.
+ *
+ * @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
+ */
+template <class T>
+T radiant_flux(T t, T a);
+
+/**
+ * Calculates the spectral radiance of a blackbody at a 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.
+ *
+ * @see physics::planck::wavelength
+ */
+template <class T>
+constexpr auto spectral_radiance = planck::wavelength<T>;
+
+template <class T>
+T radiant_exitance(T t)
+{
+	const T tt = t * t;
+	return constants::stefan_boltzmann<T> * (tt * tt);
+}
+
+template <class T>
+T radiant_flux(T t, T a)
+{
+	return a * radiant_exitance(t);
+}
+
+} // namespace blackbody
+
+} // namespace light
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_BLACKBODY_HPP

src/physics/light/light.hpp

@@ -0,0 +1,35 @@
+/*
+ * 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_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_HPP
+
+namespace physics {
+
+/// Light-related functions.
+namespace light {}
+
+} // namespace physics
+
+#include "blackbody.hpp"
+#include "luminosity.hpp"
+#include "phase.hpp"
+#include "photometry.hpp"
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_HPP

src/physics/light/luminosity.hpp

@@ -0,0 +1,56 @@
+/*
+ * 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_LUMINOSITY_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_LUMINOSITY_HPP
+
+#include <cmath>
+
+namespace physics {
+namespace light {
+
+/// Luminous efficiency functions.
+namespace luminosity {
+
+/**
+ * Fitted Gaussian approximation to the CIE 1931 standard observer photopic luminosity function.
+ *
+ * @param lambda Wavelength of light, in nanometers.
+ * @return Luminous efficiency, on `[0, 1]`.
+ *
+ * @see Wyman, C., Sloan, P.J., & Shirley, P. (2013). Simple Analytic Approximations to the CIE XYZ Color Matching Functions.
+ */
+template <class T>
+T photopic(T lambda)
+{
+	const T t0 = (lambda - T(568.8)) * ((lambda < T(568.8)) ? T(0.0213) : T(0.0247));
+	const T t1 = (lambda - T(530.9)) * ((lambda < T(530.9)) ? T(0.0613) : T(0.0322));
+	
+	const T y0 = T(0.821) * std::exp(T(-0.5) * t0 * t0);
+	const T y1 = T(0.286) * std::exp(T(-0.5) * t1 * t1);
+	
+	return y0 + y1;
+}
+
+} // namespace luminosity
+
+} // namespace light
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_LUMINOSITY_HPP

src/physics/light/phase.hpp

@@ -0,0 +1,78 @@
+/*
+ * 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_PHASE_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_PHASE_HPP
+
+#include "math/constants.hpp"
+#include <cmath>
+
+namespace physics {
+namespace light {
+
+/// Light-scattering phase functions.
+namespace phase {
+
+/**
+ * Henyey–Greenstein phase function.
+ *
+ * @param mu Cosine of the angle between the light and view directions.
+ * @param g Asymmetry factor, on `[-1, 1]`. Positive values cause forward scattering, negative values cause back scattering.
+ *
+ * @see http://www.pbr-book.org/3ed-2018/Volume_Scattering/Phase_Functions.html
+ */
+template <class T>
+T henyey_greenstein(T mu, T g);
+
+/**
+ * Isotropic phase function. Causes light to be scattered uniformly in all directions.
+ */
+template <class T>
+constexpr T isotropic()
+{
+	return T(1) / (T(4) * math::pi<T>);
+}
+
+/**
+ * Rayleigh phase function.
+ *
+ * @param mu Cosine of the angle between the light and view directions.
+ */
+template <class T>
+T rayleigh(T mu);
+
+template <class T>
+T henyey_greenstein(T mu, T g)
+{
+	const T gg = g * g;
+	return (T(1) - gg) / (T(4) * math::pi<T> * std::pow(T(1) + gg - T(2) * g * mu, T(1.5)));
+}
+
+template <class T>
+T rayleigh(T mu)
+{
+	return (T(1) + mu * mu) * (T(3) / (T(16) * math::pi<T>));
+}
+
+} // namespace phase
+
+} // namespace light
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_PHASE_HPP

src/physics/light/photometry.hpp

@@ -0,0 +1,82 @@
+/*
+ * 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_PHOTOMETRY_HPP
+#define ANTKEEPER_PHYSICS_LIGHT_PHOTOMETRY_HPP
+
+#include "math/constants.hpp"
+#include "math/quadrature.hpp"
+
+namespace physics {
+namespace light {
+
+/**
+ * Calculates the luminous efficiency of a light source.
+ *
+ * @param spd Unary function object that returns spectral radiance given a wavelength.
+ * @param lef Unary function object that returns luminous efficiency given a wavelength.
+ * @param first,last Range of sample wavelengths.
+ * @return Luminous efficiency, on `[0, 1]`.
+ *
+ * @see physics::light::blackbody::spectral_radiance
+ * @see physics::light::luminosity::photopic
+ */
+template <class T, class UnaryOp1, class UnaryOp2, class InputIt>
+T luminous_efficiency(UnaryOp1 spd, UnaryOp2 lef, InputIt first, InputIt last)
+{
+	auto spd_lef = [spd, lef](T x) -> T
+	{
+		return spd(x) * lef(x);
+	};
+	
+	const T num = math::quadrature::trapezoid(spd_lef, first, last);
+	const T den = math::quadrature::trapezoid(spd, first, last);
+	
+	return num / den;
+}
+
+/**
+ * Calculates luminous efficacy given luminous efficiency.
+ *
+ * @param efficiency Luminous efficiency, on `[0, 1]`.
+ * @return Luminous flux, in lumens.
+ */
+template <class T>
+T luminous_efficacy(T efficiency)
+{
+	return T(683.002) * efficiency;
+}
+
+/**
+ * Converts watts (radiant flux) to lumens (luminous flux).
+ *
+ * @param radient_flux Radiance 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)
+{
+	return radient_flux * luminous_efficacy(efficiency);
+}
+
+} // namespace light
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_LIGHT_PHOTOMETRY_HPP

src/physics/physics.hpp

@@ -23,7 +23,11 @@
 /// Physics
 namespace physics {}
 
+#include "constants.hpp"
 #include "frame.hpp"
+#include "planck.hpp"
+
+#include "light/light.hpp"
 #include "orbit/orbit.hpp"
 #include "time/time.hpp"
 

src/physics/planck.hpp

@@ -0,0 +1,61 @@
+/*
+ * 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_PLANCK_HPP
+#define ANTKEEPER_PHYSICS_PLANCK_HPP
+
+#include "physics/constants.hpp"
+#include <cmath>
+
+namespace physics {
+
+/// Various forms of Planck's law.
+namespace planck {
+
+/**
+ * Wavelength variant of Planck's law.
+ *
+ * @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>
+T wavelength(T t, T lambda, T c = constants::speed_of_light<T>);
+
+template <class T>
+T wavelength(T t, T lambda, T c)
+{
+	const T hc = constants::planck<T> * c;
+	const T lambda2 = lambda * lambda;
+	
+	// First radiation constant (c1L)
+	const T c1 = T(2) * hc * c;
+	
+	// Second radiation constant (c2)
+	const T c2 = hc / constants::boltzmann<T>;
+	
+	return (c1 / (lambda2 * lambda2 * lambda)) / std::expm1(c2 / (lambda * t));
+}
+
+} // namespace planck
+
+} // namespace physics
+
+#endif // ANTKEEPER_PHYSICS_PLANCK_HPP