diff --git a/src/color/xyz.hpp b/src/color/xyz.hpp index 976c49d..f54d4f6 100644 --- a/src/color/xyz.hpp +++ b/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 T luminance(const math::vector3& 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 +math::vector3 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 +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 +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 +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 +math::vector3 match(T lambda) +{ + return math::vector3 + { + match_x(lambda), + match_y(lambda), + match_z(lambda) + }; +} + +template +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 +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 +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 math::vector3 to_acescg(const math::vector3& x) { diff --git a/src/ecs/components/atmosphere-component.hpp b/src/ecs/components/atmosphere-component.hpp index 274d64c..16cc1cf 100644 --- a/src/ecs/components/atmosphere-component.hpp +++ b/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 diff --git a/src/ecs/systems/astronomy-system.cpp b/src/ecs/systems/astronomy-system.cpp index e1ff177..687eb55 100644 --- a/src/ecs/systems/astronomy-system.cpp +++ b/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 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 +T transmittance(math::vector3 start, math::vector3 end, std::size_t sample_count, T scale_height) +{ + const T inverse_scale_height = T(1) / -scale_height; + + math::vector3 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 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(reference_body)) + { + // Get the atmosphere component of the reference body + const auto& atmosphere = registry.get(reference_body); + + const double meters_per_au = 1.496e+11; + const double earth_radius_au = 4.26352e-5; + const double earth_radius_m = earth_radius_au * meters_per_au; + const double observer_altitude_m = (observer_location[0] - earth_radius_au) * meters_per_au; + + // Altitude of observer in meters + geom::ray sample_ray; + sample_ray.origin = {0, observer_altitude_m, 0}; + sample_ray.direction = math::normalize(sun_position_topocentric); + + geom::sphere 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 * sun_radius * sun_radius; + + // Calculate distance attenuation + double sun_distance_m = math::length(sun_position_topocentric) * meters_per_au; + double distance_attenuation = 1.0 / (sun_distance_m * sun_distance_m); + + double sun_luminous_flux = blackbody_luminous_flux(sun_temperature, sun_radius); + double sun_luminous_intensity = sun_luminous_flux / (4.0 * math::pi); + double sun_illuminance = sun_luminous_intensity / (sun_distance_m * sun_distance_m); + + std::cout << "distance atten: " << distance_attenuation << std::endl; + std::cout << "scatter atten: " << scattering_attenuation << std::endl; + + std::cout << "luminous flux: " << sun_luminous_flux << std::endl; + std::cout << "luminous intensity: " << sun_luminous_intensity << std::endl; + std::cout << "illuminance: " << sun_illuminance * scattering_mean << std::endl; + + + // 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(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(t, x, physics::constants::speed_of_light); + }; + + // Luminous efficiency function (photopic) + auto lef = [](double x) -> double + { + return physics::light::luminosity::photopic(x); + }; + + // Construct range of spectral sample points + std::vector samples(10000); + std::iota(samples.begin(), samples.end(), 10); + + // Calculate luminous efficiency + const double efficiency = physics::light::luminous_efficiency(spd, lef, samples.begin(), samples.end()); + + // Calculate surface area of spherical blackbody + const double a = double(4) * math::pi * r * r; + + // 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(radiant_flux, efficiency); + + return luminous_flux; +} + } // namespace ecs diff --git a/src/ecs/systems/astronomy-system.hpp b/src/ecs/systems/astronomy-system.hpp index c9c657b..a7f104e 100644 --- a/src/ecs/systems/astronomy-system.hpp +++ b/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; diff --git a/src/game/states/play-state.cpp b/src/game/states/play-state.cpp index 75ca1df..fd8cee1 100644 --- a/src/game/states/play-state.cpp +++ b/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; transform.warp = true; ecs_registry.assign(earth_entity, orbit); + ecs_registry.assign(earth_entity, atmosphere); ecs_registry.assign(earth_entity, transform); } diff --git a/src/geom/intersection.hpp b/src/geom/intersection.hpp index 08ca36b..71f7ee6 100644 --- a/src/geom/intersection.hpp +++ b/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 @@ -44,6 +45,25 @@ std::tuple ray_aabb_intersection(const ray& ray, cons std::tuple ray_mesh_intersection(const ray& ray, const mesh& mesh); +/** + * Ray-sphere intersection test. + */ +template +std::tuple ray_sphere_intersection(const ray& ray, const sphere& 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& aabb, const float3& center, float radius); bool aabb_aabb_intersection(const aabb& a, const aabb& b); diff --git a/src/math/math.hpp b/src/math/math.hpp index 0b7860b..91ecc12 100644 --- a/src/math/math.hpp +++ 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" diff --git a/src/math/quadrature.hpp b/src/math/quadrature.hpp new file mode 100644 index 0000000..0b9fc88 --- /dev/null +++ b/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 . + */ + +#ifndef ANTKEEPER_MATH_QUADRATURE_HPP +#define ANTKEEPER_MATH_QUADRATURE_HPP + +#include +#include + +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 +typename std::invoke_result::value_type>::type trapezoid(UnaryOp f, InputIt first, InputIt last) +{ + typedef typename std::invoke_result::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 diff --git a/src/physics/constants.hpp b/src/physics/constants.hpp new file mode 100644 index 0000000..3316590 --- /dev/null +++ b/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 . + */ + +#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 +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 +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 +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 +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 +constexpr T speed_of_light = T(299792458); + +} // namespace constant + +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_CONSTANTS_HPP diff --git a/src/physics/light/blackbody.hpp b/src/physics/light/blackbody.hpp new file mode 100644 index 0000000..95760cb --- /dev/null +++ b/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 . + */ + +#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 +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 +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 +constexpr auto spectral_radiance = planck::wavelength; + +template +T radiant_exitance(T t) +{ + const T tt = t * t; + return constants::stefan_boltzmann * (tt * tt); +} + +template +T radiant_flux(T t, T a) +{ + return a * radiant_exitance(t); +} + +} // namespace blackbody + +} // namespace light +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_LIGHT_BLACKBODY_HPP diff --git a/src/physics/light/light.hpp b/src/physics/light/light.hpp new file mode 100644 index 0000000..a186497 --- /dev/null +++ b/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 . + */ + +#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 diff --git a/src/physics/light/luminosity.hpp b/src/physics/light/luminosity.hpp new file mode 100644 index 0000000..0a1492c --- /dev/null +++ b/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 . + */ + +#ifndef ANTKEEPER_PHYSICS_LIGHT_LUMINOSITY_HPP +#define ANTKEEPER_PHYSICS_LIGHT_LUMINOSITY_HPP + +#include + +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 +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 diff --git a/src/physics/light/phase.hpp b/src/physics/light/phase.hpp new file mode 100644 index 0000000..6b1c833 --- /dev/null +++ b/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 . + */ + +#ifndef ANTKEEPER_PHYSICS_LIGHT_PHASE_HPP +#define ANTKEEPER_PHYSICS_LIGHT_PHASE_HPP + +#include "math/constants.hpp" +#include + +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 +T henyey_greenstein(T mu, T g); + +/** + * Isotropic phase function. Causes light to be scattered uniformly in all directions. + */ +template +constexpr T isotropic() +{ + return T(1) / (T(4) * math::pi); +} + +/** + * Rayleigh phase function. + * + * @param mu Cosine of the angle between the light and view directions. + */ +template +T rayleigh(T mu); + +template +T henyey_greenstein(T mu, T g) +{ + const T gg = g * g; + return (T(1) - gg) / (T(4) * math::pi * std::pow(T(1) + gg - T(2) * g * mu, T(1.5))); +} + +template +T rayleigh(T mu) +{ + return (T(1) + mu * mu) * (T(3) / (T(16) * math::pi)); +} + +} // namespace phase + +} // namespace light +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_LIGHT_PHASE_HPP diff --git a/src/physics/light/photometry.hpp b/src/physics/light/photometry.hpp new file mode 100644 index 0000000..3ea61af --- /dev/null +++ b/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 . + */ + +#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 +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 +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 +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 diff --git a/src/physics/physics.hpp b/src/physics/physics.hpp index bf71126..1aa33d9 100644 --- a/src/physics/physics.hpp +++ b/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" diff --git a/src/physics/planck.hpp b/src/physics/planck.hpp new file mode 100644 index 0000000..facc549 --- /dev/null +++ b/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 . + */ + +#ifndef ANTKEEPER_PHYSICS_PLANCK_HPP +#define ANTKEEPER_PHYSICS_PLANCK_HPP + +#include "physics/constants.hpp" +#include + +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 +T wavelength(T t, T lambda, T c = constants::speed_of_light); + +template +T wavelength(T t, T lambda, T c) +{ + const T hc = constants::planck * 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; + + return (c1 / (lambda2 * lambda2 * lambda)) / std::expm1(c2 / (lambda * t)); +} + +} // namespace planck + +} // namespace physics + +#endif // ANTKEEPER_PHYSICS_PLANCK_HPP