Add celestial time functions

CMakeLists.txt

@@ -4,7 +4,6 @@ option(VERSION_STRING "Project version string" "0.0.0")
 
 project(antkeeper VERSION ${VERSION_STRING} LANGUAGES CXX)
 
-
 # Find dependency packages
 find_package(dr_wav REQUIRED CONFIG)
 find_package(stb REQUIRED CONFIG)

src/game/astronomy/celestial-coordinates.cpp

@@ -58,14 +58,14 @@ double3x3 ecliptic_to_equatorial(double ecl)
 	};
 }
 
-double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst)
+double3x3 ecliptic_to_horizontal(double ecl, double lat, double lmst)
 {
 	const double c_ecl = std::cos(ecl);
 	const double s_ecl = std::sin(ecl);
 	const double c_lat = std::cos(lat);
 	const double s_lat = std::sin(lat);
-	const double c_lst = std::cos(lst);
-	const double s_lst = std::sin(lst);
+	const double c_lst = std::cos(lmst);
+	const double s_lst = std::sin(lmst);
 	
 	return double3x3
 	{
@@ -88,12 +88,12 @@ double3x3 equatorial_to_ecliptic(double ecl)
 	};
 }
 
-double3x3 equatorial_to_horizontal(double lat, double lst)
+double3x3 equatorial_to_horizontal(double lat, double lmst)
 {
 	const double c_lat = std::cos(lat);
 	const double s_lat = std::sin(lat);
-	const double c_lst = std::cos(lst);
-	const double s_lst = std::sin(lst);
+	const double c_lst = std::cos(lmst);
+	const double s_lst = std::sin(lmst);
 	
 	return double3x3
 	{
@@ -103,12 +103,12 @@ double3x3 equatorial_to_horizontal(double lat, double lst)
 	};
 }
 
-double3x3 horizontal_to_equatorial(double lat, double lst)
+double3x3 horizontal_to_equatorial(double lat, double lmst)
 {
 	const double c_lat = std::cos(lat);
 	const double s_lat = std::sin(lat);
-	const double c_lst = std::cos(lst);
-	const double s_lst = std::sin(lst);
+	const double c_lst = std::cos(lmst);
+	const double s_lst = std::sin(lmst);
 	
 	return double3x3
 	{
@@ -118,14 +118,14 @@ double3x3 horizontal_to_equatorial(double lat, double lst)
 	};
 }
 
-double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst)
+double3x3 horizontal_to_ecliptic(double ecl, double lat, double lmst)
 {
 	const double c_ecl = std::cos(ecl);
 	const double s_ecl = std::sin(ecl);
 	const double c_lat = std::cos(lat);
 	const double s_lat = std::sin(lat);
-	const double c_lst = std::cos(lst);
-	const double s_lst = std::sin(lst);
+	const double c_lst = std::cos(lmst);
+	const double s_lst = std::sin(lmst);
 	
 	return double3x3
 	{

src/game/astronomy/celestial-coordinates.hpp

@@ -54,10 +54,10 @@ double3x3 ecliptic_to_equatorial(double ecl);
  *
  * @param ecl Obliquity of the ecliptic, in radians.
  * @param lat Observer's latitude, in radians.
- * @param lst Local sidereal time.
+ * @param lmst Local mean sidereal time.
  * @return Transformation matrix.
  */
-double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst);
+double3x3 ecliptic_to_horizontal(double ecl, double lat, double lmst);
 
 /**
  * Produces a matrix which transforms rectangular coordinates from equatorial space to ecliptic space.
@@ -71,29 +71,37 @@ double3x3 equatorial_to_ecliptic(double ecl);
  * Produces a matrix which transforms rectangular coordinates from equatorial space to horizontal space.
  *
  * @param lat Observer's latitude, in radians.
- * @param lst Local sidereal time.
+ * @param lmst Local mean sidereal time.
  * @return Transformation matrix.
  */
-double3x3 equatorial_to_horizontal(double lat, double lst);
+double3x3 equatorial_to_horizontal(double lat, double lmst);
 
 /**
  * Produces a matrix which transforms rectangular coordinates from horizontal space to equatorial space.
  *
  * @param lat Observer's latitude, in radians.
- * @param lst Local sidereal time.
+ * @param lmst Local mean sidereal time.
  * @return Transformation matrix.
  */
-double3x3 horizontal_to_equatorial(double lat, double lst);
+double3x3 horizontal_to_equatorial(double lat, double lmst);
 
 /**
  * Produces a matrix which transforms rectangular coordinates from horizontal space to ecliptic space.
  *
  * @param ecl Obliquity of the ecliptic, in radians.
  * @param lat Observer's latitude, in radians.
- * @param lst Local sidereal time.
+ * @param lmst Local mean sidereal time.
  * @return Transformation matrix.
  */
-double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst);
+double3x3 horizontal_to_ecliptic(double ecl, double lat, double lmst);
+
+/// Matrix which transforms coordinates from horizontal space to a right-handed coordinate system.
+constexpr double3x3 horizontal_to_right_handed =
+{
+	0.0,  0.0, 1.0,
+	1.0,  0.0, 0.0,
+	0.0, -1.0, 0.0
+};
 
 } // namespace ast
 

src/game/astronomy/celestial-mechanics.cpp

@@ -18,11 +18,115 @@
  */
 
 #include "celestial-mechanics.hpp"
+#include "math/angles.hpp"
 #include <cmath>
 
 namespace ast
 {
 
+double approx_ecliptic_obliquity(double jd)
+{
+	return math::radians<double>(23.4393 - 3.563e-7 * (jd - 2451545.0));
+}
+
+double3 approx_sun_ecliptic(double jd)
+{
+	const double t = (jd - 2451545.0) / 36525.0;
+	const double m = 6.24 + 628.302 * t;
+	
+	const double longitude = 4.895048 + 628.331951 * t + (0.033417 - 0.000084 * t) * std::sin(m) + 0.000351 * std::sin(m * 2.0);
+	const double latitude = 0.0;
+	const double distance = 1.000140 - (0.016708 - 0.000042 * t) * std::cos(m) - 0.000141 * std::cos(m * 2.0);
+	
+	double3 ecliptic;
+	ecliptic.x = distance * std::cos(longitude) * std::cos(latitude);
+	ecliptic.y = distance * std::sin(longitude) * std::cos(latitude);
+	ecliptic.z = distance * std::sin(latitude);
+	
+	return ecliptic;
+}
 
+double3 approx_moon_ecliptic(double jd)
+{
+	const double t = (jd - 2451545.0) / 36525.0;
+	const double l1 = 3.8104 + 8399.7091 * t;
+	const double m1 = 2.3554 + 8328.6911 * t;
+	const double m = 6.2300 + 628.3019 * t;
+	const double d = 5.1985 + 7771.3772 * t;
+	const double d2 = d * 2.0;
+	const double f = 1.6280 + 8433.4663 * t;
+	
+	const double longitude = l1
+		+ 0.1098 * std::sin(m1)
+		+ 0.0222 * std::sin(d2 - m1)
+		+ 0.0115 * std::sin(d2)
+		+ 0.0037 * std::sin(m1 * 2.0)
+		- 0.0032 * std::sin(m)
+		- 0.0020 * std::sin(d2)
+		+ 0.0010 * std::sin(d2 - m1 * 2.0)
+		+ 0.0010 * std::sin(d2 - m - m1)
+		+ 0.0009 * std::sin(d2 + m1)
+		+ 0.0008 * std::sin(d2 - m)
+		+ 0.0007 * std::sin(m1 - m)
+		- 0.0006 * std::sin(d)
+		- 0.0005 * std::sin(m + m1);
+
+	const double latitude = 0.0895 * sin(f)
+		+ 0.0049 * std::sin(m1 + f)
+		+ 0.0048 * std::sin(m1 - f)
+		+ 0.0030 * std::sin(d2 - f)
+		+ 0.0010 * std::sin(d2 + f - m1)
+		+ 0.0008 * std::sin(d2 - f - m1)
+		+ 0.0006 * std::sin(d2 + f);
+	
+	const double r = 1.0 / (0.016593
+		+ 0.000904 * std::cos(m1)
+		+ 0.000166 * std::cos(d2 - m1)
+		+ 0.000137 * std::cos(d2)
+		+ 0.000049 * std::cos(m1 * 2.0)
+		+ 0.000015 * std::cos(d2 + m1)
+		+ 0.000009 * std::cos(d2 - m));
+	
+	double3 ecliptic;
+	ecliptic.x = r * std::cos(longitude) * std::cos(latitude);
+	ecliptic.y = r * std::sin(longitude) * std::cos(latitude);
+	ecliptic.z = r * std::sin(latitude);
+	
+	return ecliptic;
+}
+
+double3x3 approx_moon_ecliptic_rotation(double jd)
+{
+	const double t = (jd - 2451545.0) / 36525.0;
+	const double l1 = 3.8104 + 8399.7091 * t;
+	const double f = 1.6280 + 8433.4663 * t;
+	
+	const double az0 = f + math::pi<double>;
+	const double ax  = 0.026920;
+	const double az1 = l1 - f;
+	
+	double3x3 rz0 =
+	{
+		cos(az0), -sin(az0), 0,
+		sin(az0), cos(az0), 0,
+		0, 0, 1
+	};
+	
+	double3x3 rx =
+	{
+		1, 0, 0,
+		0, cos(ax), -sin(ax),
+		0, sin(ax), cos(ax)
+	};
+	
+	double3x3 rz1 =
+	{
+		cos(az1), -sin(az1), 0,
+		sin(az1), cos(az1), 0,
+		0, 0, 1
+	};
+	
+	return rz0 * rx * rz1;
+}
 
 } // namespace ast

src/game/astronomy/celestial-mechanics.hpp

@@ -25,6 +25,38 @@
 namespace ast
 {
 
+/**
+ * Approximates the obliquity of the ecliptic.
+ *
+ * @param jd Julian date.
+ * @return Obliquity of the ecliptic, in radians.
+ */
+double approx_ecliptic_obliquity(double jd);
+
+/**
+ * Approximates the ecliptic coordinates of the Earth's sun.
+ *
+ * @param jd Julian date.
+ * @return Ecliptic rectangular geocentric coordinates of the Earth's sun, with distance in AU.
+ */
+double3 approx_sun_ecliptic(double jd);
+
+/**
+ * Approximates the ecliptic coordinates of the Earth's moon.
+ *
+ * @param jd Julian date.
+ * @return Ecliptic rectangular geocentric coordinates of the Earth's moon, with distance in Earth radii.
+ */
+double3 approx_moon_ecliptic(double jd);
+
+/**
+ * Approximates the ecliptic rotation of the Earth's moon.
+ *
+ * @param jd Julian date.
+ * @return Rotation matrix representing the moon's rotation in ecliptic space.
+ */
+double3x3 approx_moon_ecliptic_rotation(double jd);
+
 struct orbital_elements
 {
 	double a; ///< Semi-major axis (km)

src/game/astronomy/celestial-time.cpp

@@ -0,0 +1,62 @@
+/*
+ * Copyright (C) 2020  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/>.
+ */
+
+#include "celestial-time.hpp"
+#include "math/angles.hpp"
+#include <cmath>
+
+namespace ast
+{
+
+double ut_to_jd(int year, int month, int day, int hour, int minute, double second)
+{
+	if (month < 3)
+	{
+		month += 12;
+		year -= 1;
+	}
+	
+	const signed long Y = year;
+	const signed long M = month;
+	const signed long D = day;	
+	const signed long JDN = (1461*(Y+4800+(M-14)/12))/4+(367*(M-2-12*((M-14)/12)))/12-(3*((Y+4900+(M-14)/12)/100))/4+D-32075;
+	
+	const double h = static_cast<double>(hour - 12) / 24.0;
+	const double m = static_cast<double>(minute) / 1440.0;
+	const double s = second / 86400.0;
+	
+	return static_cast<double>(JDN) + h + m + s;
+}
+
+double jd_to_gmst(double jd)
+{
+	return math::wrap_radians<double>(4.894961212 + 6.300388098 * (jd - 2451545.0));
+}
+
+double jd_to_lmst(double jd, double longitude)
+{
+	return gmst_to_lmst(jd_to_gmst(jd), longitude);
+}
+
+double gmst_to_lmst(double gmst, double longitude)
+{
+	return gmst + longitude;
+}
+
+} // namespace ast

src/game/astronomy/celestial-time.hpp

@@ -0,0 +1,66 @@
+/*
+ * Copyright (C) 2020  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_CELESTIAL_TIME_HPP
+#define ANTKEEPER_CELESTIAL_TIME_HPP
+
+namespace ast
+{
+
+/**
+ * Calculates the Julian date from a universal time.
+ *
+ * @param year Gregorian year, with 1 BC as 0.
+ * @param month Gregorian month, on [1, 12].
+ * @param day Gregorian day, on [1, 31].
+ * @param hour Hour, on [0, 23].
+ * @param minute Minute, on [0, 59].
+ * @param second Second on [0, 60).
+ * @return Julian date with fractional day.
+ */
+double ut_to_jd(int year, int month, int day, int hour, int minute, double second);
+
+/**
+ * Calculates the Greenwich mean sidereal time (GMST) from a Julian date.
+ *
+ * @param jd Julian date.
+ * @return GMST, in radians.
+ */
+double jd_to_gmst(double jd);
+
+/**
+ * Calculates local mean sidereal time (LMST) from a Julian date.
+ *
+ * @param jd Julian date.
+ * @param longitude Longitude of the observer, in radians.
+ */
+double jd_to_lmst(double jd, double longitude);
+
+/**
+ * Calculates local mean sidereal time (LMST) from Greenwich mean sidereal time (GMST).
+ *
+ * @param gmst Greenwich mean sidereal time, in radians.
+ * @param longitude Longitude of the observer, in radians.
+ * @return Local mean sidereal time, in radians.
+ */
+double gmst_to_lmst(double gmst, double longitude);
+
+} // namespace ast
+
+#endif // ANTKEEPER_CELESTIAL_TIME_HPP

src/game/states/play-state.cpp

@@ -81,7 +81,7 @@ void play_state_enter(game_context* ctx)
 	sky_pass->set_moon_model(ctx->resource_manager->load<model>("moon.mdl"));
 	
 	ctx->weather_system->set_location(ctx->biome->location[0], ctx->biome->location[1], ctx->biome->location[2]);
-	ctx->weather_system->set_time(2017, 8, 21, 5, 0, 0, -7.0);
+	ctx->weather_system->set_time(2017, 8, 21, 5, 0, 0.0, -7.0);
 	ctx->weather_system->set_sky_palette(ctx->biome->sky_palette);
 	ctx->weather_system->set_sun_palette(ctx->biome->sun_palette);
 	ctx->weather_system->set_ambient_palette(ctx->biome->ambient_palette);
@@ -118,7 +118,7 @@ void play_state_enter(game_context* ctx)
 	marker_archetype->assign(ecs_registry, ctx->marker_entity);
 	container_archetype->assign(ecs_registry, ctx->container_entity);
 	twig_archetype->assign(ecs_registry, ctx->twig_entity);
-		
+	
 	// Create flashlight and light cone, set light cone parent to flashlight, and move both to underworld scene
 	flashlight_archetype->assign(ecs_registry, ctx->flashlight_entity);
 	auto flashlight_light_cone = flashlight_light_cone_archetype->create(ecs_registry);
@@ -137,7 +137,7 @@ void play_state_enter(game_context* ctx)
 	//ec::bind_transform(ecs_registry, lens_light_cone, ctx->lens_entity);
 	ec::parent(ecs_registry, lens_light_cone, ctx->lens_entity);
 	
-
+	
 	
 	// Hide inactive tools
 	ec::assign_render_layers(ecs_registry, ctx->forceps_entity, 0);
@@ -241,19 +241,19 @@ void play_state_enter(game_context* ctx)
 	ecs_registry.assign_or_replace<ecs::transform_component>(ctx->focal_point_entity, focal_point_transform);
 	ecs_registry.assign_or_replace<ecs::camera_follow_component>(ctx->focal_point_entity, focal_point_follow);
 	ecs_registry.assign_or_replace<ecs::snap_component>(ctx->focal_point_entity, focal_point_snap);
-
+	
 	// Setup camera
 	ctx->overworld_camera->look_at({0, 0, 1}, {0, 0, 0}, {0, 1, 0});
 	ctx->camera_system->set_camera(ctx->overworld_camera);
-
+	
 	auto ant_head = ant_head_archetype->create(ecs_registry);
 	ec::place(ecs_registry, ant_head, {50, 0});
-
+	
 	ctx->overworld_scene->update_tweens();
-
+	
 	// Allocate a nest
 	nest* nest = new ::nest();
-
+	
 	// Setup initial nest parameters
 	float tunnel_radius = 1.15f;
 	nest->set_tunnel_radius(tunnel_radius);
@@ -275,7 +275,7 @@ void play_state_enter(game_context* ctx)
 		chamber.outer_radius = 10.0f;
 		central_shaft->chambers.push_back(chamber);
 	}
-
+	
 	// Dig nest shafts
 	float shift = 0.1f;
 	for (int i = 0; i < 800; ++i)

src/game/systems/weather-system.cpp

@@ -26,6 +26,8 @@
 #include "utility/gamma.hpp"
 #include "resources/image.hpp"
 #include "game/astronomy/celestial-coordinates.hpp"
+#include "game/astronomy/celestial-mechanics.hpp"
+#include "game/astronomy/celestial-time.hpp"
 #include <cmath>
 #include <iostream>
 
@@ -33,190 +35,6 @@ static constexpr double hours_per_day = 24.0;
 static constexpr double minutes_per_day = hours_per_day * 60.0;
 static constexpr double seconds_per_day = minutes_per_day * 60.0;
 
-/**
- *
- * @param year Gregorian year
- * @param month Month (1 = January, 12 = December)
- * @param day Day (1-31)
- * @param time Universal time in decimal hours.
- */
-static double julian_day(int year, int month, int day, double time)
-{
-	if (month < 3)
-	{
-		month += 12;
-		year -= 1;
-	}
-	
-	double y = static_cast<double>(year);
-	double m = static_cast<double>(month);
-	double d = static_cast<double>(day);
-	
-	return std::floor(365.25 * y) + std::floor(30.6001 * (m + 1.0)) - 15.0 + 1720996.5 + d + time;
-}
-
-/**
- * Calculates the ecliptic rectangular geocentric coordinates of the sun, with distance in AU.
- */
-double3 calculate_sun_ecliptic(double jd)
-{
-	const double t = (jd - 2451545.0) / 36525.0;
-	const double m = 6.24 + 628.302 * t;
-	
-	const double longitude = 4.895048 + 628.331951 * t + (0.033417 - 0.000084 * t) * std::sin(m) + 0.000351 * std::sin(m * 2.0);
-	const double latitude = 0.0;
-	const double distance = 1.000140 - (0.016708 - 0.000042 * t) * std::cos(m) - 0.000141 * std::cos(m * 2.0);
-	
-	double3 ecliptic;
-	ecliptic.x = distance * std::cos(longitude) * std::cos(latitude);
-	ecliptic.y = distance * std::sin(longitude) * std::cos(latitude);
-	ecliptic.z = distance * std::sin(latitude);
-	
-	return ecliptic;
-}
-
-/**
- * Calculates the ecliptic rectangular geocentric coordinates of the moon, with distance in Earth radii.
- */
-double3 calculate_moon_ecliptic(double jd)
-{
-	const double t = (jd - 2451545.0) / 36525.0;
-	const double l1 = 3.8104 + 8399.7091 * t;
-	const double m1 = 2.3554 + 8328.6911 * t;
-	const double m = 6.2300 + 628.3019 * t;
-	const double d = 5.1985 + 7771.3772 * t;
-	const double d2 = d * 2.0;
-	const double f = 1.6280 + 8433.4663 * t;
-	
-	const double longitude = l1
-		+ 0.1098 * std::sin(m1)
-		+ 0.0222 * std::sin(d2 - m1)
-		+ 0.0115 * std::sin(d2)
-		+ 0.0037 * std::sin(m1 * 2.0)
-		- 0.0032 * std::sin(m)
-		- 0.0020 * std::sin(d2)
-		+ 0.0010 * std::sin(d2 - m1 * 2.0)
-		+ 0.0010 * std::sin(d2 - m - m1)
-		+ 0.0009 * std::sin(d2 + m1)
-		+ 0.0008 * std::sin(d2 - m)
-		+ 0.0007 * std::sin(m1 - m)
-		- 0.0006 * std::sin(d)
-		- 0.0005 * std::sin(m + m1);
-
-	const double latitude = 0.0895 * sin(f)
-		+ 0.0049 * std::sin(m1 + f)
-		+ 0.0048 * std::sin(m1 - f)
-		+ 0.0030 * std::sin(d2 - f)
-		+ 0.0010 * std::sin(d2 + f - m1)
-		+ 0.0008 * std::sin(d2 - f - m1)
-		+ 0.0006 * std::sin(d2 + f);
-	
-	const double r = 1.0 / (0.016593
-		+ 0.000904 * std::cos(m1)
-		+ 0.000166 * std::cos(d2 - m1)
-		+ 0.000137 * std::cos(d2)
-		+ 0.000049 * std::cos(m1 * 2.0)
-		+ 0.000015 * std::cos(d2 + m1)
-		+ 0.000009 * std::cos(d2 - m));
-	
-	double3 ecliptic;
-	ecliptic.x = r * std::cos(longitude) * std::cos(latitude);
-	ecliptic.y = r * std::sin(longitude) * std::cos(latitude);
-	ecliptic.z = r * std::sin(latitude);
-	
-	return ecliptic;
-}
-
-double3x3 find_moon_ecliptic_rotation(double jd)
-{
-	const double t = (jd - 2451545.0) / 36525.0;
-	const double l1 = 3.8104 + 8399.7091 * t;
-	const double f = 1.6280 + 8433.4663 * t;
-	
-	const double az0 = f + math::pi<double>;
-	const double ax  = 0.026920;
-	const double az1 = l1 - f;
-	
-	double3x3 rz0 =
-	{
-		cos(az0), -sin(az0), 0,
-		sin(az0), cos(az0), 0,
-		0, 0, 1
-	};
-	
-	double3x3 rx =
-	{
-		1, 0, 0,
-		0, cos(ax), -sin(ax),
-		0, sin(ax), cos(ax)
-	};
-	
-	double3x3 rz1 =
-	{
-		cos(az1), -sin(az1), 0,
-		sin(az1), cos(az1), 0,
-		0, 0, 1
-	};
-	
-	return rz0 * rx * rz1;
-}
-
-/// @see http://www.stjarnhimlen.se/comp/ppcomp.html
-/// @see http://www.geoastro.de/elevazmoon/basics/index.htm
-void ecliptic_to_equatorial(double longitude, double latitude, double ecl, double* right_ascension, double* declination)
-{
-	double eclip_x = std::cos(longitude) * std::cos(latitude);
-	double eclip_y = std::sin(longitude) * std::cos(latitude);
-	double eclip_z = std::sin(latitude);
-	
-	double equat_x = eclip_x;
-	double equat_y = eclip_y * std::cos(ecl) - eclip_z * std::sin(ecl);
-	double equat_z = eclip_y * std::sin(ecl) + eclip_z * std::cos(ecl);
-	
-    *right_ascension = std::atan2(equat_y, equat_x);
-    *declination = std::atan2(equat_z, sqrt(equat_x * equat_x + equat_y * equat_y));
-}
-
-/// @see http://www.stjarnhimlen.se/comp/ppcomp.html
-/// @see http://www.geoastro.de/elevazmoon/basics/index.htm
-void equatorial_to_horizontal(double right_ascension, double declination, double lmst, double latitude, double* azimuth, double* elevation)
-{
-	double hour_angle = lmst - right_ascension;
-	
-	double x = std::cos(hour_angle) * std::cos(declination);
-	double y = std::sin(hour_angle) * std::cos(declination);
-	double z = std::sin(declination);
-	
-	double horiz_x = x * std::cos(math::half_pi<double> - latitude) - z * std::sin(math::half_pi<double> - latitude);
-	double horiz_y = y;
-	double horiz_z = x * std::sin(math::half_pi<double> - latitude) + z * std::cos(math::half_pi<double> - latitude);
-	
-	*azimuth = math::wrap_radians<double>(std::atan2(horiz_y, horiz_x) + math::pi<double>);
-	*elevation = math::wrap_radians<double>(std::atan2(horiz_z, std::sqrt(horiz_x * horiz_x + horiz_y * horiz_y)));
-}
-
-
-double3x3 horizontal_to_right_handed()
-{
-	return double3x3
-	{
-		 0.0,  0.0,  1.0,
-		 1.0,  0.0,  0.0,
-		 0.0, -1.0,  0.0
-	};
-}
-
-/**
- * Calculates the Greenwich mean sidereal time (GMST) from a Julian day.
- *
- * @param jd Julian day.
- * @return GMST, in radians.
- */
-static double jd_to_gmst(double jd)
-{	 
-	return math::wrap_radians<double>(4.894961212 + 6.300388098 * (jd - 2451545.0));
-}
-
 weather_system::weather_system(entt::registry& registry):
 	entity_system(registry),
 	ambient_light(nullptr),
@@ -239,42 +57,38 @@ void weather_system::update(double t, double dt)
 	const float latitude = location[0];
 	const float longitude = location[1];
 	
-	// Calculate local mean sidereal time (LMST)
+	// Calculate local time
 	double time_correction = longitude / (math::two_pi<double> / 24.0);
 	double local_jd = jd + time_correction / 24.0 - 0.5;
 	double local_time = (local_jd - std::floor(local_jd)) * 24.0;
-	double hour = local_time;
-	double gmst = jd_to_gmst(jd);
-	double lmst = gmst + longitude;
-	
-	// Obliquity of the ecliptic
-	double ecl = math::radians<double>(23.4393 - 3.563e-7 * (jd - 2451545.0));
 	
 	// Solar distance in AU
 	//double sr = ...
 	// Apparent radius in degrees
 	//double sradius = 0.2666 / sr;
 	
+	double lmst = ast::jd_to_lmst(jd, longitude);
+	double ecl = ast::approx_ecliptic_obliquity(jd);
 	double3x3 ecliptic_to_horizontal = ast::ecliptic_to_horizontal(ecl, latitude, lmst);
+
 	
-	double3 sun_ecliptic = calculate_sun_ecliptic(jd);
+	double3 sun_ecliptic = ast::approx_sun_ecliptic(jd);
 	double3 sun_horizontal = ecliptic_to_horizontal * sun_ecliptic;
+	sun_horizontal.z -= 4.25875e-5; // Subtract one earth radius (in AU), for positon of observer
 	double3 sun_spherical = ast::rectangular_to_spherical(sun_horizontal);
-	double3 sun_positiond = horizontal_to_right_handed() * sun_horizontal;
+	double3 sun_positiond = ast::horizontal_to_right_handed * sun_horizontal;
 	float2 sun_az_el = {static_cast<float>(sun_spherical.z) - math::pi<float>, static_cast<float>(sun_spherical.y)};
 	float3 sun_position = math::normalize(float3{static_cast<float>(sun_positiond.x), static_cast<float>(sun_positiond.y), static_cast<float>(sun_positiond.z)});
 	
-	double3 moon_ecliptic = calculate_moon_ecliptic(jd);
+	double3 moon_ecliptic = ast::approx_moon_ecliptic(jd);
 	double3 moon_horizontal = ecliptic_to_horizontal * moon_ecliptic;
 	moon_horizontal.z -= 1.0; // Subtract one earth radius, for position of observer
 	double3 moon_spherical = ast::rectangular_to_spherical(moon_horizontal);
-	double3 moon_positiond = horizontal_to_right_handed() * moon_horizontal;
+	double3 moon_positiond = ast::horizontal_to_right_handed * moon_horizontal;
 	float2 moon_az_el = {static_cast<float>(moon_spherical.z) - math::pi<float>, static_cast<float>(moon_spherical.y)};
 	float3 moon_position = math::normalize(float3{static_cast<float>(moon_positiond.x), static_cast<float>(moon_positiond.y), static_cast<float>(moon_positiond.z)});
-	
-	//std::cout << "new moon: " << math::degrees(moon_az_el[0]) << ", " << math::degrees(moon_az_el[1]) << std::endl;
-	
-	double3x3 moon_rotation_matrix = horizontal_to_right_handed() * ecliptic_to_horizontal;
+		
+	double3x3 moon_rotation_matrix = ast::horizontal_to_right_handed * ecliptic_to_horizontal;
 	math::quaternion<double> moon_rotationd = math::normalize(math::quaternion_cast(moon_rotation_matrix) * math::angle_axis(math::half_pi<double>, double3{0, 1, 0}) * math::angle_axis(-math::half_pi<double>, double3{0, 0, -1}));
 	math::quaternion<float> moon_rotation =
 	{
@@ -300,9 +114,6 @@ void weather_system::update(double t, double dt)
 		moon_light->set_rotation(moon_az_el_rotation);
 	}
 	
-	std::size_t hour_index = static_cast<std::size_t>(hour);
-	float lerp_factor = hour - std::floor(hour);
-	
 	float sun_gradient_position = static_cast<float>(std::max<double>(0.0, ((sun_az_el[1] + math::half_pi<double>) / math::pi<double>)));
 	float moon_gradient_position = static_cast<float>(std::max<double>(0.0, ((moon_az_el[1] + math::half_pi<double>) / math::pi<double>)));
 	float sky_gradient_position = sun_gradient_position;
@@ -325,7 +136,7 @@ void weather_system::update(double t, double dt)
 		
 		sky_pass->set_horizon_color(horizon_color);
 		sky_pass->set_zenith_color(zenith_color);
-		sky_pass->set_time_of_day(static_cast<float>(hour * 60.0 * 60.0));
+		sky_pass->set_time_of_day(static_cast<float>(local_time * 60.0 * 60.0));
 		sky_pass->set_observer_location(location[0], location[1], location[2]);
 		sky_pass->set_sun_coordinates(sun_position, sun_az_el);
 		sky_pass->set_moon_coordinates(moon_position, moon_az_el);
@@ -399,10 +210,9 @@ void weather_system::set_material_pass(::material_pass* pass)
 	material_pass = pass;
 }
 
-void weather_system::set_time(int year, int month, int day, int hour, int minute, int second, double tc)
+void weather_system::set_time(int year, int month, int day, int hour, int minute, double second, double tc)
 {
-	double time = ((static_cast<double>(hour) - tc) + ((static_cast<double>(minute) + static_cast<double>(second) / 60.0) / 60.0)) / 24.0;
-	jd = julian_day(year, month, day, time);
+	jd = ast::ut_to_jd(year, month, day, hour, minute, second) - tc / 24.0;
 }
 
 void weather_system::set_time_scale(float scale)

src/game/systems/weather-system.hpp

@@ -53,7 +53,7 @@ public:
 	void set_material_pass(::material_pass* pass);
 	
 	/// @param tc Timezone correction, in hours
-	void set_time(int year, int month, int day, int hour, int minute, int second, double tc);
+	void set_time(int year, int month, int day, int hour, int minute, double second, double tc);
 	void set_time_scale(float scale);
 	
 	void set_sky_palette(const ::image* image);