💿🐜 Antkeeper source code https://antkeeper.com
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

628 lines
24 KiB

  1. /*
  2. * Copyright (C) 2021 Christopher J. Howard
  3. *
  4. * This file is part of Antkeeper source code.
  5. *
  6. * Antkeeper source code is free software: you can redistribute it and/or modify
  7. * it under the terms of the GNU General Public License as published by
  8. * the Free Software Foundation, either version 3 of the License, or
  9. * (at your option) any later version.
  10. *
  11. * Antkeeper source code is distributed in the hope that it will be useful,
  12. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14. * GNU General Public License for more details.
  15. *
  16. * You should have received a copy of the GNU General Public License
  17. * along with Antkeeper source code. If not, see <http://www.gnu.org/licenses/>.
  18. */
  19. #include "entity/systems/astronomy.hpp"
  20. #include "astro/apparent-size.hpp"
  21. #include "entity/components/blackbody.hpp"
  22. #include "entity/components/transform.hpp"
  23. #include "entity/components/diffuse-reflector.hpp"
  24. #include "geom/intersection.hpp"
  25. #include "geom/cartesian.hpp"
  26. #include "color/color.hpp"
  27. #include "physics/orbit/orbit.hpp"
  28. #include "physics/time/ut1.hpp"
  29. #include "physics/time/jd.hpp"
  30. #include "physics/light/photometry.hpp"
  31. #include "physics/light/luminosity.hpp"
  32. #include "physics/light/refraction.hpp"
  33. #include "physics/gas/atmosphere.hpp"
  34. #include "geom/cartesian.hpp"
  35. #include "astro/apparent-size.hpp"
  36. #include "geom/solid-angle.hpp"
  37. #include "math/polynomial.hpp"
  38. #include <iostream>
  39. namespace entity {
  40. namespace system {
  41. astronomy::astronomy(entity::registry& registry):
  42. updatable(registry),
  43. time_days(0.0),
  44. time_centuries(0.0),
  45. time_scale(1.0),
  46. observer_eid(entt::null),
  47. reference_body_eid(entt::null),
  48. transmittance_samples(0),
  49. sun_light(nullptr),
  50. sky_light(nullptr),
  51. moon_light(nullptr),
  52. bounce_light(nullptr),
  53. bounce_albedo{0, 0, 0},
  54. sky_pass(nullptr),
  55. starlight_illuminance{0, 0, 0}
  56. {
  57. // Construct ENU to EUS transformation
  58. enu_to_eus = math::transformation::se3<double>
  59. {
  60. {0, 0, 0},
  61. math::quaternion<double>::rotate_x(-math::half_pi<double>)
  62. };
  63. registry.on_construct<entity::component::observer>().connect<&astronomy::on_observer_modified>(this);
  64. registry.on_replace<entity::component::observer>().connect<&astronomy::on_observer_modified>(this);
  65. registry.on_destroy<entity::component::observer>().connect<&astronomy::on_observer_destroyed>(this);
  66. registry.on_construct<entity::component::celestial_body>().connect<&astronomy::on_celestial_body_modified>(this);
  67. registry.on_replace<entity::component::celestial_body>().connect<&astronomy::on_celestial_body_modified>(this);
  68. registry.on_destroy<entity::component::celestial_body>().connect<&astronomy::on_celestial_body_destroyed>(this);
  69. registry.on_construct<entity::component::orbit>().connect<&astronomy::on_orbit_modified>(this);
  70. registry.on_replace<entity::component::orbit>().connect<&astronomy::on_orbit_modified>(this);
  71. registry.on_destroy<entity::component::orbit>().connect<&astronomy::on_orbit_destroyed>(this);
  72. registry.on_construct<entity::component::atmosphere>().connect<&astronomy::on_atmosphere_modified>(this);
  73. registry.on_replace<entity::component::atmosphere>().connect<&astronomy::on_atmosphere_modified>(this);
  74. registry.on_destroy<entity::component::atmosphere>().connect<&astronomy::on_atmosphere_destroyed>(this);
  75. }
  76. astronomy::~astronomy()
  77. {
  78. registry.on_construct<entity::component::observer>().disconnect<&astronomy::on_observer_modified>(this);
  79. registry.on_replace<entity::component::observer>().disconnect<&astronomy::on_observer_modified>(this);
  80. registry.on_destroy<entity::component::observer>().disconnect<&astronomy::on_observer_destroyed>(this);
  81. registry.on_construct<entity::component::celestial_body>().disconnect<&astronomy::on_celestial_body_modified>(this);
  82. registry.on_replace<entity::component::celestial_body>().disconnect<&astronomy::on_celestial_body_modified>(this);
  83. registry.on_destroy<entity::component::celestial_body>().disconnect<&astronomy::on_celestial_body_destroyed>(this);
  84. registry.on_construct<entity::component::orbit>().disconnect<&astronomy::on_orbit_modified>(this);
  85. registry.on_replace<entity::component::orbit>().disconnect<&astronomy::on_orbit_modified>(this);
  86. registry.on_destroy<entity::component::orbit>().disconnect<&astronomy::on_orbit_destroyed>(this);
  87. registry.on_construct<entity::component::atmosphere>().disconnect<&astronomy::on_atmosphere_modified>(this);
  88. registry.on_replace<entity::component::atmosphere>().disconnect<&astronomy::on_atmosphere_modified>(this);
  89. registry.on_destroy<entity::component::atmosphere>().disconnect<&astronomy::on_atmosphere_destroyed>(this);
  90. }
  91. void astronomy::update(double t, double dt)
  92. {
  93. double3 sky_light_illuminance = {0.0, 0.0, 0.0};
  94. // Add scaled timestep to current time
  95. set_time(time_days + dt * time_scale);
  96. // Abort if no valid observer entity or reference body entity
  97. if (observer_eid == entt::null || reference_body_eid == entt::null)
  98. return;
  99. // Get pointer to observer component
  100. const auto observer = registry.try_get<component::observer>(observer_eid);
  101. // Abort if no observer component
  102. if (!observer)
  103. return;
  104. // Get pointers to reference body components
  105. const auto
  106. [
  107. reference_body,
  108. reference_orbit,
  109. reference_atmosphere
  110. ] = registry.try_get<component::celestial_body, component::orbit, component::atmosphere>(reference_body_eid);
  111. // Abort if no reference body or reference orbit
  112. if (!reference_body || !reference_orbit)
  113. return;
  114. // Update ICRF to EUS transformation
  115. update_icrf_to_eus(*reference_body, *reference_orbit);
  116. // Set the transform component translations of orbiting bodies to their topocentric positions
  117. registry.view<component::celestial_body, component::orbit, component::transform>().each(
  118. [&](entity::id entity_id, const auto& body, const auto& orbit, auto& transform)
  119. {
  120. // Skip reference body entity
  121. if (entity_id == reference_body_eid)
  122. return;
  123. // Transform orbital Cartesian position (r) from the ICRF frame to the EUS frame
  124. const double3 r_eus = icrf_to_eus * orbit.position;
  125. // Evaluate body orientation polynomials
  126. const double body_pole_ra = math::polynomial::horner(body.pole_ra.begin(), body.pole_ra.end(), time_centuries);
  127. const double body_pole_dec = math::polynomial::horner(body.pole_dec.begin(), body.pole_dec.end(), time_centuries);
  128. const double body_prime_meridian = math::polynomial::horner(body.prime_meridian.begin(), body.prime_meridian.end(), time_days);
  129. // Determine body orientation in the ICRF frame
  130. math::quaternion<double> rotation_icrf = physics::orbit::frame::bcbf::to_bci
  131. (
  132. body_pole_ra,
  133. body_pole_dec,
  134. body_prime_meridian
  135. ).r;
  136. // Transform body orientation from the ICRF frame to the EUS frame.
  137. math::quaternion<double> rotation_eus = math::normalize(icrf_to_eus.r * rotation_icrf);
  138. // Update local transform
  139. if (orbit.parent != entt::null)
  140. {
  141. transform.local.translation = math::normalize(math::type_cast<float>(r_eus));
  142. transform.local.rotation = math::type_cast<float>(rotation_eus);
  143. transform.local.scale = {1.0f, 1.0f, 1.0f};
  144. }
  145. });
  146. constexpr double3 bounce_normal = {0, 1, 0};
  147. double3 bounce_illuminance = {0, 0, 0};
  148. // Update blackbody lighting
  149. registry.view<component::celestial_body, component::orbit, component::blackbody>().each(
  150. [&](entity::id entity_id, const auto& blackbody_body, const auto& blackbody_orbit, const auto& blackbody)
  151. {
  152. // Transform blackbody position from ICRF frame to EUS frame
  153. const double3 blackbody_position_eus = icrf_to_eus * blackbody_orbit.position;
  154. // Measure distance and direction, in EUS frame, from observer to blackbody
  155. const double observer_blackbody_distance = math::length(blackbody_position_eus);
  156. const double3 observer_blackbody_direction_eus = blackbody_position_eus / observer_blackbody_distance;
  157. // Measure blackbody solid angle as seen by observer
  158. const double observer_blackbody_angular_radius = astro::angular_radius(blackbody_body.radius, observer_blackbody_distance);
  159. const double observer_blackbody_solid_angle = geom::solid_angle::cone(observer_blackbody_angular_radius);
  160. // Calculate illuminance from blackbody reaching observer
  161. const double3 observer_blackbody_illuminance = blackbody.luminance * observer_blackbody_solid_angle;
  162. // Calculate illuminance from blackbody reaching observer after atmospheric extinction
  163. double3 observer_blackbody_transmitted_illuminance = observer_blackbody_illuminance;
  164. if (reference_atmosphere)
  165. {
  166. // Construct ray at observer pointing towards the blackbody
  167. const geom::ray<double> ray = {{0, 0, 0}, observer_blackbody_direction_eus};
  168. // Integrate atmospheric spectral transmittance factor between observer and blackbody
  169. const double3 transmittance = integrate_transmittance(*observer, *reference_body, *reference_atmosphere, ray);
  170. // Attenuate illuminance from blackbody reaching observer by spectral transmittance factor
  171. observer_blackbody_transmitted_illuminance *= transmittance;
  172. }
  173. // Update sun light
  174. if (sun_light != nullptr)
  175. {
  176. const double3 blackbody_up_eus = icrf_to_eus.r * double3{0, 0, 1};
  177. sun_light->set_rotation
  178. (
  179. math::look_rotation
  180. (
  181. math::type_cast<float>(-observer_blackbody_direction_eus),
  182. math::type_cast<float>(blackbody_up_eus)
  183. )
  184. );
  185. sun_light->set_color(math::type_cast<float>(observer_blackbody_transmitted_illuminance));
  186. // Bounce sun light
  187. bounce_illuminance += std::max(0.0, math::dot(bounce_normal, -observer_blackbody_direction_eus)) * observer_blackbody_transmitted_illuminance * bounce_albedo;
  188. }
  189. // Update sky light
  190. if (sky_light != nullptr)
  191. {
  192. // Calculate sky illuminance
  193. double3 blackbody_position_enu_spherical = physics::orbit::frame::enu::spherical(enu_to_eus.inverse() * blackbody_position_eus);
  194. const double sky_illuminance = 25000.0 * std::max<double>(0.0, std::sin(blackbody_position_enu_spherical.y));
  195. // Add sky illuminance to sky light illuminance
  196. sky_light_illuminance += {sky_illuminance, sky_illuminance, sky_illuminance};
  197. // Add starlight illuminance to sky light illuminance
  198. sky_light_illuminance += starlight_illuminance;
  199. // Update sky light
  200. sky_light->set_color(math::type_cast<float>(sky_light_illuminance));
  201. // Bounce sky light
  202. bounce_illuminance += sky_light_illuminance * bounce_albedo;
  203. }
  204. // Upload blackbody params to sky pass
  205. if (this->sky_pass)
  206. {
  207. this->sky_pass->set_sun_position(math::type_cast<float>(blackbody_position_eus));
  208. this->sky_pass->set_sun_luminance(math::type_cast<float>(blackbody.luminance));
  209. this->sky_pass->set_sun_illuminance(math::type_cast<float>(observer_blackbody_illuminance), math::type_cast<float>(observer_blackbody_transmitted_illuminance));
  210. this->sky_pass->set_sun_angular_radius(static_cast<float>(observer_blackbody_angular_radius));
  211. }
  212. // Update diffuse reflectors
  213. this->registry.view<component::celestial_body, component::orbit, component::diffuse_reflector, component::transform>().each(
  214. [&](entity::id entity_id, const auto& reflector_body, const auto& reflector_orbit, const auto& reflector, const auto& transform)
  215. {
  216. // Transform reflector position from ICRF frame to EUS frame
  217. const double3 reflector_position_eus = icrf_to_eus * reflector_orbit.position;
  218. // Measure distance and direction, in EUS frame, from observer to reflector
  219. const double observer_reflector_distance = math::length(reflector_position_eus);
  220. const double3 observer_reflector_direction_eus = reflector_position_eus / observer_reflector_distance;
  221. // Measure distance and direction, in EUS frame, from reflector to blackbody
  222. double3 reflector_blackbody_direction_eus = blackbody_position_eus - reflector_position_eus;
  223. const double reflector_blackbody_distance = math::length(reflector_blackbody_direction_eus);
  224. reflector_blackbody_direction_eus /= reflector_blackbody_distance;
  225. // Measure blackbody solid angle as seen by reflector
  226. const double reflector_blackbody_angular_radius = astro::angular_radius(blackbody_body.radius, reflector_blackbody_distance);
  227. const double reflector_blackbody_solid_angle = geom::solid_angle::cone(reflector_blackbody_angular_radius);
  228. // Calculate blackbody illuminance reaching reflector
  229. const double3 reflector_blackbody_illuminance = blackbody.luminance * reflector_blackbody_solid_angle;
  230. // Measure reflector solid angle as seen by observer
  231. const double observer_reflector_angular_radius = astro::angular_radius(reflector_body.radius, observer_reflector_distance);
  232. const double observer_reflector_solid_angle = geom::solid_angle::cone(observer_reflector_angular_radius);
  233. // Determine phase factor of reflector as seen by observer
  234. const double observer_reflector_phase_factor = dot(observer_reflector_direction_eus, -reflector_blackbody_direction_eus) * 0.5 + 0.5;
  235. // Measure observer reference body solid angle as seen by reflector
  236. const double reflector_observer_angular_radius = astro::angular_radius(reference_body->radius, observer_reflector_distance);
  237. const double reflector_observer_solid_angle = geom::solid_angle::cone(reflector_observer_angular_radius);
  238. // Determine phase factor of observer reference body as by reflector
  239. const double reflector_observer_phase_factor = dot(-observer_reflector_direction_eus, -observer_blackbody_direction_eus) * 0.5 + 0.5;
  240. // Calculate spectral transmittance between observer and reflector factor due to atmospheric extinction
  241. double3 observer_reflector_transmittance = {1, 1, 1};
  242. if (reference_atmosphere)
  243. {
  244. const geom::ray<double> ray = {{0, 0, 0}, observer_reflector_direction_eus};
  245. observer_reflector_transmittance = integrate_transmittance(*observer, *reference_body, *reference_atmosphere, ray);
  246. }
  247. // Measure luminance of observer reference body as seen by reflector
  248. const double3 reflector_observer_luminance = observer_blackbody_illuminance * reference_body->albedo * observer_reflector_transmittance * reflector_observer_phase_factor * math::inverse_pi<double>;
  249. // Measure illuminance from observer reference body reaching reflector
  250. const double3 reflector_observer_illuminance = reflector_observer_luminance * reflector_observer_solid_angle;
  251. // Measure luminance of reflector as seen by observer
  252. const double3 observer_reflector_luminance = (reflector_blackbody_illuminance * observer_reflector_phase_factor + reflector_observer_illuminance) * reflector.albedo * observer_reflector_transmittance * math::inverse_pi<double>;
  253. // Measure illuminance from reflector reaching observer
  254. const double3 observer_reflector_illuminance = observer_reflector_luminance * observer_reflector_solid_angle;
  255. if (this->sky_pass)
  256. {
  257. this->sky_pass->set_moon_position(transform.local.translation);
  258. this->sky_pass->set_moon_rotation(transform.local.rotation);
  259. this->sky_pass->set_moon_angular_radius(static_cast<float>(observer_reflector_angular_radius));
  260. this->sky_pass->set_moon_sunlight_direction(math::type_cast<float>(-reflector_blackbody_direction_eus));
  261. this->sky_pass->set_moon_sunlight_illuminance(math::type_cast<float>(reflector_blackbody_illuminance * observer_reflector_transmittance));
  262. this->sky_pass->set_moon_planetlight_direction(math::type_cast<float>(observer_reflector_direction_eus));
  263. this->sky_pass->set_moon_planetlight_illuminance(math::type_cast<float>(reflector_observer_illuminance * observer_reflector_transmittance));
  264. this->sky_pass->set_moon_illuminance(math::type_cast<float>(observer_reflector_illuminance / observer_reflector_transmittance), math::type_cast<float>(observer_reflector_illuminance));
  265. }
  266. if (this->moon_light)
  267. {
  268. const float3 reflector_up_eus = math::type_cast<float>(icrf_to_eus.r * double3{0, 0, 1});
  269. this->moon_light->set_color(math::type_cast<float>(observer_reflector_illuminance));
  270. this->moon_light->set_rotation
  271. (
  272. math::look_rotation
  273. (
  274. math::type_cast<float>(-observer_reflector_direction_eus),
  275. reflector_up_eus
  276. )
  277. );
  278. // Bounce moon light
  279. bounce_illuminance += std::max(0.0, math::dot(bounce_normal, -observer_reflector_direction_eus)) * observer_reflector_illuminance * bounce_albedo;
  280. }
  281. });
  282. });
  283. if (bounce_light)
  284. {
  285. bounce_light->set_color(math::type_cast<float>(bounce_illuminance));
  286. }
  287. }
  288. void astronomy::set_time(double t)
  289. {
  290. time_days = t;
  291. time_centuries = time_days * physics::time::jd::centuries_per_day<double>;
  292. }
  293. void astronomy::set_time_scale(double scale)
  294. {
  295. time_scale = scale;
  296. }
  297. void astronomy::set_observer(entity::id eid)
  298. {
  299. if (observer_eid != eid)
  300. {
  301. observer_eid = eid;
  302. if (observer_eid != entt::null)
  303. observer_modified();
  304. else
  305. reference_body_eid = entt::null;
  306. }
  307. }
  308. void astronomy::set_transmittance_samples(std::size_t samples)
  309. {
  310. transmittance_samples = samples;
  311. }
  312. void astronomy::set_sun_light(scene::directional_light* light)
  313. {
  314. sun_light = light;
  315. }
  316. void astronomy::set_sky_light(scene::ambient_light* light)
  317. {
  318. sky_light = light;
  319. }
  320. void astronomy::set_moon_light(scene::directional_light* light)
  321. {
  322. moon_light = light;
  323. }
  324. void astronomy::set_bounce_light(scene::directional_light* light)
  325. {
  326. bounce_light = light;
  327. }
  328. void astronomy::set_bounce_albedo(const double3& albedo)
  329. {
  330. bounce_albedo = albedo;
  331. }
  332. void astronomy::set_starlight_illuminance(const double3& illuminance)
  333. {
  334. starlight_illuminance = illuminance;
  335. }
  336. void astronomy::set_sky_pass(::render::sky_pass* pass)
  337. {
  338. this->sky_pass = pass;
  339. if (sky_pass)
  340. {
  341. if (observer_eid != entt::null)
  342. {
  343. // Get pointer to observer
  344. const auto observer = registry.try_get<component::observer>(reference_body_eid);
  345. sky_pass->set_observer_elevation(static_cast<float>(observer->elevation));
  346. }
  347. if (reference_body_eid != entt::null)
  348. {
  349. // Get pointer to reference celestial body
  350. const auto reference_body = registry.try_get<component::celestial_body>(reference_body_eid);
  351. if (reference_body)
  352. sky_pass->set_planet_radius(static_cast<float>(reference_body->radius));
  353. else
  354. sky_pass->set_planet_radius(0.0f);
  355. }
  356. }
  357. }
  358. void astronomy::on_observer_modified(entity::registry& registry, entity::id entity_id, entity::component::observer& component)
  359. {
  360. if (entity_id == observer_eid)
  361. observer_modified();
  362. }
  363. void astronomy::on_observer_destroyed(entity::registry& registry, entity::id entity_id)
  364. {
  365. if (entity_id == observer_eid)
  366. observer_modified();
  367. }
  368. void astronomy::on_celestial_body_modified(entity::registry& registry, entity::id entity_id, entity::component::celestial_body& component)
  369. {
  370. if (entity_id == reference_body_eid)
  371. reference_body_modified();
  372. }
  373. void astronomy::on_celestial_body_destroyed(entity::registry& registry, entity::id entity_id)
  374. {
  375. if (entity_id == reference_body_eid)
  376. reference_body_modified();
  377. }
  378. void astronomy::on_orbit_modified(entity::registry& registry, entity::id entity_id, entity::component::orbit& component)
  379. {
  380. if (entity_id == reference_body_eid)
  381. reference_orbit_modified();
  382. }
  383. void astronomy::on_orbit_destroyed(entity::registry& registry, entity::id entity_id)
  384. {
  385. if (entity_id == reference_body_eid)
  386. reference_orbit_modified();
  387. }
  388. void astronomy::on_atmosphere_modified(entity::registry& registry, entity::id entity_id, entity::component::atmosphere& component)
  389. {
  390. if (entity_id == reference_body_eid)
  391. reference_atmosphere_modified();
  392. }
  393. void astronomy::on_atmosphere_destroyed(entity::registry& registry, entity::id entity_id)
  394. {
  395. if (entity_id == reference_body_eid)
  396. reference_atmosphere_modified();
  397. }
  398. void astronomy::observer_modified()
  399. {
  400. // Get pointer to observer component
  401. const auto observer = registry.try_get<component::observer>(observer_eid);
  402. if (observer)
  403. {
  404. if (reference_body_eid != observer->reference_body_eid)
  405. {
  406. // Reference body changed
  407. reference_body_eid = observer->reference_body_eid;
  408. reference_body_modified();
  409. reference_orbit_modified();
  410. reference_atmosphere_modified();
  411. }
  412. if (reference_body_eid != entt::null)
  413. {
  414. // Get pointer to reference celestial body
  415. const auto reference_body = registry.try_get<component::celestial_body>(reference_body_eid);
  416. // Update BCBF to EUS transformation
  417. if (reference_body)
  418. update_bcbf_to_eus(*observer, *reference_body);
  419. }
  420. // Upload observer elevation to sky pass
  421. if (sky_pass)
  422. sky_pass->set_observer_elevation(static_cast<float>(observer->elevation));
  423. }
  424. }
  425. void astronomy::reference_body_modified()
  426. {
  427. // Get pointer to reference celestial body
  428. const auto reference_body = registry.try_get<component::celestial_body>(reference_body_eid);
  429. if (reference_body)
  430. {
  431. // Get pointer to observer
  432. const auto observer = registry.try_get<component::observer>(observer_eid);
  433. // Update BCBF to EUS transformation
  434. if (observer)
  435. update_bcbf_to_eus(*observer, *reference_body);
  436. }
  437. // Update reference celestial body-related sky pass parameters
  438. if (sky_pass)
  439. {
  440. if (reference_body)
  441. sky_pass->set_planet_radius(static_cast<float>(reference_body->radius));
  442. else
  443. sky_pass->set_planet_radius(0.0f);
  444. }
  445. }
  446. void astronomy::reference_orbit_modified()
  447. {
  448. }
  449. void astronomy::reference_atmosphere_modified()
  450. {
  451. }
  452. void astronomy::update_bcbf_to_eus(const entity::component::observer& observer, const entity::component::celestial_body& body)
  453. {
  454. // Construct BCBF to EUS transformation
  455. bcbf_to_eus = physics::orbit::frame::bcbf::to_enu
  456. (
  457. body.radius + observer.elevation,
  458. observer.latitude,
  459. observer.longitude
  460. ) * enu_to_eus;
  461. }
  462. void astronomy::update_icrf_to_eus(const entity::component::celestial_body& body, const entity::component::orbit& orbit)
  463. {
  464. // Evaluate reference body orientation polynomials
  465. const double body_pole_ra = math::polynomial::horner(body.pole_ra.begin(), body.pole_ra.end(), time_centuries);
  466. const double body_pole_dec = math::polynomial::horner(body.pole_dec.begin(), body.pole_dec.end(), time_centuries);
  467. const double body_prime_meridian = math::polynomial::horner(body.prime_meridian.begin(), body.prime_meridian.end(), time_days);
  468. // Construct ICRF frame to BCBF transformation
  469. math::transformation::se3<double> icrf_to_bcbf = physics::orbit::frame::bci::to_bcbf
  470. (
  471. body_pole_ra,
  472. body_pole_dec,
  473. body_prime_meridian
  474. );
  475. icrf_to_bcbf.t = icrf_to_bcbf.r * -orbit.position;
  476. /// Construct ICRF to EUS transformation
  477. icrf_to_eus = icrf_to_bcbf * bcbf_to_eus;
  478. // Pass ICRF to EUS transformation to sky pass
  479. if (sky_pass)
  480. {
  481. // Upload topocentric frame to sky pass
  482. sky_pass->set_icrf_to_eus
  483. (
  484. math::transformation::se3<float>
  485. {
  486. math::type_cast<float>(icrf_to_eus.t),
  487. math::type_cast<float>(icrf_to_eus.r)
  488. }
  489. );
  490. }
  491. }
  492. double3 astronomy::integrate_transmittance(const entity::component::observer& observer, const entity::component::celestial_body& body, const entity::component::atmosphere& atmosphere, geom::ray<double> ray) const
  493. {
  494. double3 transmittance = {1, 1, 1};
  495. // Make ray height relative to center of reference body
  496. ray.origin.y += body.radius + observer.elevation;
  497. // Construct sphere representing upper limit of the atmosphere
  498. geom::sphere<double> atmosphere_sphere;
  499. atmosphere_sphere.center = {0, 0, 0};
  500. atmosphere_sphere.radius = body.radius + atmosphere.upper_limit;
  501. // Check for intersection between the ray and atmosphere
  502. auto intersection = geom::ray_sphere_intersection(ray, atmosphere_sphere);
  503. if (std::get<0>(intersection))
  504. {
  505. // Get point of intersection
  506. const double3 intersection_point = ray.extrapolate(std::get<2>(intersection));
  507. // Integrate optical of Rayleigh, Mie, and ozone particles
  508. const double optical_depth_r = physics::gas::atmosphere::optical_depth_exp(ray.origin, intersection_point, body.radius, atmosphere.rayleigh_scale_height, transmittance_samples);
  509. const double optical_depth_m = physics::gas::atmosphere::optical_depth_exp(ray.origin, intersection_point, body.radius, atmosphere.mie_scale_height, transmittance_samples);
  510. const double optical_depth_o = physics::gas::atmosphere::optical_depth_tri(ray.origin, intersection_point, body.radius, atmosphere.ozone_lower_limit, atmosphere.ozone_upper_limit, atmosphere.ozone_mode, transmittance_samples);
  511. // Calculate transmittance factor due to scattering and absorption
  512. const double3 extinction_r = atmosphere.rayleigh_scattering * optical_depth_r;
  513. const double extinction_m = atmosphere.mie_extinction * optical_depth_m;
  514. const double3 extinction_o = atmosphere.ozone_absorption * optical_depth_o;
  515. transmittance = extinction_r + double3{extinction_m, extinction_m, extinction_m} + extinction_o;
  516. transmittance.x = std::exp(-transmittance.x);
  517. transmittance.y = std::exp(-transmittance.y);
  518. transmittance.z = std::exp(-transmittance.z);
  519. }
  520. return transmittance;
  521. }
  522. } // namespace system
  523. } // namespace entity