Add functions to generate celestial coordinate conversion matrices, and improve astronomical calculations

CMakeLists.txt

@@ -4,6 +4,7 @@ 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)
@@ -14,7 +15,6 @@ find_package(SDL2 REQUIRED COMPONENTS SDL2::SDL2-static SDL2::SDL2main CONFIG)
 find_package(OpenAL REQUIRED CONFIG)
 find_library(physfs REQUIRED NAMES physfs-static PATHS "${CMAKE_PREFIX_PATH}/lib")
 
-
 # Determine dependencies
 set(STATIC_LIBS
 	dr_wav

src/game/astronomy/celestial-coordinates.cpp

@@ -0,0 +1,138 @@
+/*
+ * 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-coordinates.hpp"
+#include <cmath>
+
+namespace ast
+{
+
+double3 rectangular_to_spherical(const double3& rectangular)
+{
+	const double xx_yy = rectangular.x * rectangular.x + rectangular.y * rectangular.y;
+	
+	return double3
+	{
+		std::sqrt(xx_yy + rectangular.z * rectangular.z),
+		std::atan2(rectangular.z, std::sqrt(xx_yy)),
+		std::atan2(rectangular.y, rectangular.x)
+	};
+}
+
+double3 spherical_to_rectangular(const double3& spherical)
+{
+	return double3
+	{
+		spherical[0] * std::sin(spherical[1]) * std::cos(spherical[2]),
+		spherical[0] * std::sin(spherical[1]) * std::sin(spherical[2]),
+		spherical[0] * std::cos(spherical[1]),
+	};
+}
+
+double3x3 ecliptic_to_equatorial(double ecl)
+{
+	const double c_ecl = std::cos(ecl);
+	const double s_ecl = std::sin(ecl);
+	
+	return double3x3
+	{
+		1.0, 0.0, 0.0,
+		0.0, c_ecl, s_ecl,
+		0.0, -s_ecl, c_ecl
+	};
+}
+
+double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst)
+{
+	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);
+	
+	return double3x3
+	{
+		s_lat * c_lst, s_lst, c_lat * c_lst,
+		s_lat * s_lst * c_ecl - c_lat * s_ecl, -c_lst * c_ecl, c_lat * s_lst * c_ecl + s_lat * s_ecl,
+		s_lat * s_lst * -s_ecl - c_lat * c_ecl,  c_lst * s_ecl, c_lat * s_lst * -s_ecl + s_lat * c_ecl
+	};
+}
+
+double3x3 equatorial_to_ecliptic(double ecl)
+{
+	const double c_ecl = std::cos(ecl);
+	const double s_ecl = std::sin(ecl);
+	
+	return double3x3
+	{
+		1.0, 0.0, 0.0,
+		0.0, c_ecl, -s_ecl,
+		0.0, s_ecl, c_ecl
+	};
+}
+
+double3x3 equatorial_to_horizontal(double lat, double lst)
+{
+	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);
+	
+	return double3x3
+	{
+		s_lat * c_lst, s_lst, c_lat * c_lst,
+		s_lat * s_lst, -c_lst, c_lat * s_lst,
+		-c_lat, 0.0, s_lat
+	};
+}
+
+double3x3 horizontal_to_equatorial(double lat, double lst)
+{
+	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);
+	
+	return double3x3
+	{
+		c_lst * s_lat, s_lst * s_lat, -c_lat,
+		s_lst, c_lst, 0.0,
+		c_lst * c_lat, s_lst * c_lat, s_lat
+	};
+}
+
+double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst)
+{
+	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);
+	
+	return double3x3
+	{
+		s_lat * c_lst, s_lat * s_lst * c_ecl - c_lat * s_ecl, s_lat * s_lst * -s_ecl - c_lat * c_ecl,
+		s_lst, -c_lst * c_ecl, c_lst * s_ecl,
+		c_lat * c_lst, c_lat * s_lst * c_ecl + s_lat * s_ecl, c_lat * s_lst * -s_ecl + s_lat * c_ecl
+	};
+}
+
+} // namespace ast

src/game/astronomy/celestial-coordinates.hpp

@@ -0,0 +1,100 @@
+/*
+ * 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_COORDINATES_HPP
+#define ANTKEEPER_CELESTIAL_COORDINATES_HPP
+
+#include "utility/fundamental-types.hpp"
+
+namespace ast
+{
+
+/**
+ * Converts rectangular coordinates to spherical coordinates.
+ *
+ * @param rectangular Rectangular coordinates.
+ * @return Equivalent spherical coordinates, in the ISO order of radial distance, polar angle (radians), and azimuthal angle (radians).
+ */
+double3 rectangular_to_spherical(const double3& rectangular);
+
+/**
+ * Convert spherical coordinates to rectangular coordinates.
+ *
+ * @param spherical Spherical coordinates, in the ISO order of radial distance, polar angle (radians), and azimuthal angle (radians).
+ * @return Equivalent rectangular coordinates.
+ */
+double3 spherical_to_rectangular(const double3& spherical);
+
+/**
+ * Produces a matrix which transforms rectangular coordinates from ecliptic space to equatorial space.
+ *
+ * @param ecl Obliquity of the ecliptic, in radians.
+ * @return Transformation matrix.
+ */
+double3x3 ecliptic_to_equatorial(double ecl);
+
+/**
+ * Produces a matrix which transforms coordinates from ecliptic space to horizontal space.
+ *
+ * @param ecl Obliquity of the ecliptic, in radians.
+ * @param lat Observer's latitude, in radians.
+ * @param lst Local sidereal time.
+ * @return Transformation matrix.
+ */
+double3x3 ecliptic_to_horizontal(double ecl, double lat, double lst);
+
+/**
+ * Produces a matrix which transforms rectangular coordinates from equatorial space to ecliptic space.
+ *
+ * @param ecl Obliquity of the ecliptic, in radians.
+ * @return Transformation matrix.
+ */
+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.
+ * @return Transformation matrix.
+ */
+double3x3 equatorial_to_horizontal(double lat, double lst);
+
+/**
+ * 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.
+ * @return Transformation matrix.
+ */
+double3x3 horizontal_to_equatorial(double lat, double lst);
+
+/**
+ * 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.
+ * @return Transformation matrix.
+ */
+double3x3 horizontal_to_ecliptic(double ecl, double lat, double lst);
+
+} // namespace ast
+
+#endif // ANTKEEPER_CELESTIAL_COORDINATES_HPP

src/game/astronomy/celestial-mechanics.cpp

@@ -0,0 +1,28 @@
+/*
+ * 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-mechanics.hpp"
+#include <cmath>
+
+namespace ast
+{
+
+
+
+} // namespace ast

src/game/astronomy/celestial-mechanics.hpp

@@ -0,0 +1,48 @@
+/*
+ * 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_MECHANICS_HPP
+#define ANTKEEPER_CELESTIAL_MECHANICS_HPP
+
+#include "utility/fundamental-types.hpp"
+
+namespace ast
+{
+
+struct orbital_elements
+{
+	double a; ///< Semi-major axis (km)
+	double e; ///< Eccentricity
+	double w; ///< Argument of periapsis (radians)
+	double m; ///< Mean anomaly (radians)
+	double i; ///< Inclination to the ecliptic (radians)
+	double n; ///< Longitude of the ascending node (radians)
+};
+
+//constexpr orbital_elements j2k_orbit_moon = {384400.0, 0.0554, 318.15, 135.275.16, 125.08};
+
+struct orbital_state
+{
+	double3 r; ///< Cartesian position
+	double3 v; ///< Cartesian velocity
+};
+
+} // namespace ast
+
+#endif // ANTKEEPER_CELESTIAL_MECHANICS_HPP

src/game/biome.hpp

@@ -36,17 +36,6 @@ struct biome
 	material* terrain_material;
 	
 	// Weather
-	float3 ambient_color;
-	float ambient_intensity;
-	float sun_azimuth;
-	float sun_elevation;
-	float3 sun_color;
-	float sun_intensity;
-	float sun_angular_radius;
-	float3 horizon_color;
-	float3 zenith_color;
-	float wind_speed;
-	float wind_direction;
 	image* sky_palette;
 	image* sun_palette;
 	image* moon_palette;

src/game/states/play-state.cpp

@@ -75,24 +75,13 @@ void play_state_enter(game_context* ctx)
 	}
 	
 	// Apply biome parameters to scene
-	ctx->sun_indirect->set_color(ctx->biome->ambient_color);
-	ctx->sun_indirect->set_intensity(ctx->biome->ambient_intensity);
-	
-	math::quaternion<float> sun_azimuth_rotation = math::angle_axis(ctx->biome->sun_azimuth, float3{0, 1, 0});
-	math::quaternion<float> sun_elevation_rotation = math::angle_axis(ctx->biome->sun_elevation, float3{-1, 0, 0});
-	math::quaternion<float> sun_rotation = math::normalize(sun_azimuth_rotation * sun_elevation_rotation);
-	
-	ctx->sun_direct->set_rotation(sun_rotation);
-	ctx->sun_direct->set_color(ctx->biome->sun_color);
-	ctx->sun_direct->set_intensity(ctx->biome->sun_intensity);
-	
 	sky_pass* sky_pass = ctx->overworld_sky_pass;
 	sky_pass->set_enabled(true);
 	sky_pass->set_sky_model(ctx->resource_manager->load<model>("sky-dome.mdl"));
 	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(2020, 6, 1, 5, 0, 0, -7.0);
+	ctx->weather_system->set_time(2017, 8, 21, 5, 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);

src/game/systems/weather-system.cpp

@@ -25,6 +25,7 @@
 #include "renderer/passes/material-pass.hpp"
 #include "utility/gamma.hpp"
 #include "resources/image.hpp"
+#include "game/astronomy/celestial-coordinates.hpp"
 #include <cmath>
 #include <iostream>
 
@@ -54,31 +55,30 @@ static double julian_day(int year, int month, int day, double time)
 	return std::floor(365.25 * y) + std::floor(30.6001 * (m + 1.0)) - 15.0 + 1720996.5 + d + time;
 }
 
-/// @see A Physically-Based Night Sky Model
-/// @see http://www.powerfromthesun.net/Book/chapter03/chapter03.html
-void find_sun_ecliptic(double jd, double* longitude, double* latitude, double* distance)
+/**
+ * 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;
 	
-	*longitude = 4.895048 + 628.331951 * t + (0.033417 - 0.000084 * t) * std::sin(m) + 0.000351 * std::sin(m * 2.0);
-	*latitude = 0.0;
-	*distance = 1.000140 - (0.016708 - 0.000042 * t) * std::cos(m) - 0.000141 * std::cos(m * 2.0);
+	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 geocentric coordinates of the moon, given a Julian day.
- *
- * @param[in] jd Julian day.
- * @param[out] longitude Ecliptic longitude of the moon, in radians.
- * @param[out] latitude Ecliptic latitude of the moon, in radians.
- * @param[out] distance Distance to the moon, in Earth radii.
- 
- * @return Array containing the ecliptic longitude and latitude of the moon, in radians.
- *
- * @see A Physically-Based Night Sky Model
+ * Calculates the ecliptic rectangular geocentric coordinates of the moon, with distance in Earth radii.
  */
-void find_moon_ecliptic(double jd, double* longitude, double* latitude, double* distance)
+double3 calculate_moon_ecliptic(double jd)
 {
 	const double t = (jd - 2451545.0) / 36525.0;
 	const double l1 = 3.8104 + 8399.7091 * t;
@@ -88,7 +88,7 @@ void find_moon_ecliptic(double jd, double* longitude, double* latitude, double*
 	const double d2 = d * 2.0;
 	const double f = 1.6280 + 8433.4663 * t;
 	
-	*longitude = l1
+	const double longitude = l1
 		+ 0.1098 * std::sin(m1)
 		+ 0.0222 * std::sin(d2 - m1)
 		+ 0.0115 * std::sin(d2)
@@ -103,7 +103,7 @@ void find_moon_ecliptic(double jd, double* longitude, double* latitude, double*
 		- 0.0006 * std::sin(d)
 		- 0.0005 * std::sin(m + m1);
 
-	*latitude = 0.0895 * sin(f)
+	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)
@@ -111,13 +111,54 @@ void find_moon_ecliptic(double jd, double* longitude, double* latitude, double*
 		+ 0.0008 * std::sin(d2 - f - m1)
 		+ 0.0006 * std::sin(d2 + f);
 	
-	*distance = 1.0 / (0.016593
+	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
@@ -154,6 +195,17 @@ void equatorial_to_horizontal(double right_ascension, double declination, 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.
  *
@@ -187,97 +239,72 @@ void weather_system::update(double t, double dt)
 	const float latitude = location[0];
 	const float longitude = location[1];
 	
-	// Time correction
-	double tc = longitude / (math::two_pi<double> / 24.0);
-	
-	//double pst_tc = -7.0;
-	
-	double local_jd = jd + tc / 24.0 - 0.5;
+	// Calculate local mean sidereal time (LMST)
+	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;
-	
-	
-
-	
-	// Calculate equation of time
-	//float eot_b = (360.0f / 365.0f) * (day_of_year - 81.0f);
-	//float eot = 9.87f * std::sin(eot_b * 2.0f) - 7.53f * std::cos(eot_b) - 1.5f * std::sin(eot_b);
-	
-	// Calculate local mean sidereal time (LST)
-	//double tc = longitude / (math::two_pi<double> / 24.0); // Time correction
-	//double ut = local_time + tc;                           // Universal time
 	double gmst = jd_to_gmst(jd);
 	double lmst = gmst + longitude;
 	
-	// Calculate sun position
-	//float local_solar_time = local_time;// + eot / 60.0f;
-	/*
-	float sun_declination = math::radians(23.45f) * std::sin((math::two_pi<float> / 365.0f) * (284.0f + day_of_year));
-	float sun_hour_angle = math::radians(15.0f) * (local_solar_time - 12.0f);
-	sun_elevation = std::asin(std::sin(sun_declination) * std::sin(latitude) + std::cos(sun_declination) * std::cos(sun_hour_angle) * std::cos(latitude));
-	sun_azimuth = std::acos((std::sin(sun_declination) * std::cos(latitude) - std::cos(sun_declination) * std::cos(sun_hour_angle) * std::sin(latitude)) / std::cos(sun_elevation));
-	if (sun_hour_angle > 0.0f)
-		sun_azimuth = math::two_pi<float> - sun_azimuth;
-	*/
-	
-	// J2000 day
-	double d = jd - 2451545.0;
-	
 	// Obliquity of the ecliptic
-	double ecl = math::radians<double>(23.4393 - 3.563e-7 * d);
-	
-	// Calculation sun coordinates
-	double sun_longitude;
-	double sun_latitude;
-	double sun_distance;
-	double sun_right_ascension;
-	double sun_declination;
-	double sun_azimuth;
-	double sun_elevation;
-	find_sun_ecliptic(jd, &sun_longitude, &sun_latitude, &sun_distance);
-	ecliptic_to_equatorial(sun_longitude, sun_latitude, ecl, &sun_right_ascension, &sun_declination);
-	equatorial_to_horizontal(sun_right_ascension, sun_declination, lmst, latitude, &sun_azimuth, &sun_elevation);
-	
-	// Calculate moon coordinates
-	double moon_longitude;
-	double moon_latitude;
-	double moon_distance;
-	double moon_right_ascension;
-	double moon_declination;
-	double moon_azimuth;
-	double moon_elevation;
-	find_moon_ecliptic(jd, &moon_longitude, &moon_latitude, &moon_distance);
-	ecliptic_to_equatorial(moon_longitude, moon_latitude, ecl, &moon_right_ascension, &moon_declination);
-	equatorial_to_horizontal(moon_right_ascension, moon_declination, lmst, latitude, &moon_azimuth, &moon_elevation);
-	
-	float2 sun_az_el = float2{static_cast<float>(sun_azimuth), static_cast<float>(sun_elevation)};
-	math::quaternion<float> sun_azimuth_rotation = math::angle_axis(sun_az_el[0], float3{0, 1, 0});
-	math::quaternion<float> sun_elevation_rotation = math::angle_axis(sun_az_el[1], float3{-1, 0, 0});
-	math::quaternion<float> sun_rotation = math::normalize(sun_azimuth_rotation * sun_elevation_rotation);
-	float3 sun_position = math::normalize(sun_rotation * float3{0, 0, -1});
-	
-	float2 moon_az_el = float2{static_cast<float>(moon_azimuth), static_cast<float>(moon_elevation)};
-	math::quaternion<float> moon_azimuth_rotation = math::angle_axis(moon_az_el[0], float3{0, 1, 0});
-	math::quaternion<float> moon_elevation_rotation = math::angle_axis(moon_az_el[1], float3{-1, 0, 0});
-	math::quaternion<float> moon_rotation = math::normalize(moon_azimuth_rotation * moon_elevation_rotation);
-	float3 moon_position = math::normalize(moon_rotation * float3{0, 0, -1});
+	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;
+	
+	double3x3 ecliptic_to_horizontal = ast::ecliptic_to_horizontal(ecl, latitude, lmst);
+	
+	double3 sun_ecliptic = calculate_sun_ecliptic(jd);
+	double3 sun_horizontal = ecliptic_to_horizontal * sun_ecliptic;
+	double3 sun_spherical = ast::rectangular_to_spherical(sun_horizontal);
+	double3 sun_positiond = 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_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;
+	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;
+	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 =
+	{
+		static_cast<float>(moon_rotationd.w),
+		static_cast<float>(moon_rotationd.x),
+		static_cast<float>(moon_rotationd.y),
+		static_cast<float>(moon_rotationd.z)
+	};
 	
 	if (sun_light)
 	{
-		sun_light->set_rotation(sun_rotation);
+		math::quaternion<float> sun_azimuth_rotation = math::angle_axis(sun_az_el[0], float3{0, 1, 0});
+		math::quaternion<float> sun_elevation_rotation = math::angle_axis(sun_az_el[1], float3{-1, 0, 0});
+		math::quaternion<float> sun_az_el_rotation = math::normalize(sun_azimuth_rotation * sun_elevation_rotation);
+		sun_light->set_rotation(sun_az_el_rotation);
 	}
 	
 	if (moon_light)
 	{
-		moon_light->set_rotation(moon_rotation);
+		math::quaternion<float> moon_azimuth_rotation = math::angle_axis(moon_az_el[0], float3{0, 1, 0});
+		math::quaternion<float> moon_elevation_rotation = math::angle_axis(moon_az_el[1], float3{-1, 0, 0});
+		math::quaternion<float> moon_az_el_rotation = math::normalize(moon_azimuth_rotation * moon_elevation_rotation);
+		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_elevation + math::half_pi<double>) / math::pi<double>)));
-	float moon_gradient_position = static_cast<float>(std::max<double>(0.0, ((moon_elevation + math::half_pi<double>) / math::pi<double>)));
+	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;
 	float ambient_gradient_position = sun_gradient_position;
 	
@@ -290,9 +317,11 @@ void weather_system::update(double t, double dt)
 		float3 ambient_color = interpolate_gradient(ambient_colors, ambient_gradient_position);
 		
 		sun_light->set_color(sun_color);
+		sun_light->set_intensity(1.0f);
 		moon_light->set_color(moon_color);
 		moon_light->set_intensity(1.0f);
 		ambient_light->set_color(ambient_color);
+		ambient_light->set_intensity(0.5f);
 		
 		sky_pass->set_horizon_color(horizon_color);
 		sky_pass->set_zenith_color(zenith_color);
@@ -301,12 +330,13 @@ void weather_system::update(double t, double dt)
 		sky_pass->set_sun_coordinates(sun_position, sun_az_el);
 		sky_pass->set_moon_coordinates(moon_position, moon_az_el);
 		sky_pass->set_julian_day(static_cast<float>(jd));
+		sky_pass->set_moon_rotation(moon_rotation);
 	}
 	
 	shadow_light = sun_light;
 	if (shadow_map_pass)
 	{
-		if (sun_elevation < 0.0f)
+		if (sun_az_el[1] < 0.0f)
 		{
 			shadow_map_pass->set_light(moon_light);
 		}
@@ -367,11 +397,6 @@ void weather_system::set_shadow_map_pass(::shadow_map_pass* pass)
 void weather_system::set_material_pass(::material_pass* pass)
 {
 	material_pass = pass;
-	
-	if (material_pass)
-	{
-		material_pass->set_shadow_strength(0.75f);
-	}
 }
 
 void weather_system::set_time(int year, int month, int day, int hour, int minute, int second, double tc)

src/renderer/passes/sky-pass.cpp

@@ -147,7 +147,12 @@ void sky_pass::render(render_context* context) const
 		float moon_distance = (clip_near + clip_far) * 0.5f;		
 		float moon_radius = moon_angular_radius * moon_distance;
 		
-		model = math::scale(math::translate(math::identity4x4<float>, moon_position * -moon_distance), float3{moon_radius, moon_radius, moon_radius});
+		math::transform<float> moon_transform;
+		moon_transform.translation = moon_position * -moon_distance;
+		moon_transform.rotation = moon_rotation;
+		moon_transform.scale = {moon_radius, moon_radius, moon_radius};
+		
+		model = math::matrix_cast(moon_transform);		
 		model_view = view * model;
 		model_view_projection = projection * model_view;
 		float3x3 normal_model = math::transpose(math::inverse(math::resize<3, 3>(model)));
@@ -299,6 +304,11 @@ void sky_pass::set_moon_coordinates(const float3& position, const float2& az_el)
 	moon_az_el_tween[1] = az_el;
 }
 
+void sky_pass::set_moon_rotation(const math::quaternion<float>& rotation)
+{
+	moon_rotation = rotation;
+}
+
 void sky_pass::set_moon_angular_radius(float radius)
 {
 	moon_angular_radius = radius;

src/renderer/passes/sky-pass.hpp

@@ -25,6 +25,7 @@
 #include "event/event-handler.hpp"
 #include "event/input-events.hpp"
 #include "animation/tween.hpp"
+#include "math/quaternion-type.hpp"
 
 class shader_program;
 class shader_input;
@@ -61,6 +62,7 @@ public:
 	void set_observer_location(float latitude, float longitude, float altitude);
 	void set_sun_coordinates(const float3& position, const float2& az_el);
 	void set_moon_coordinates(const float3& position, const float2& az_el);
+	void set_moon_rotation(const math::quaternion<float>& rotation);
 	
 	void set_moon_angular_radius(float radius);
 	void set_sun_angular_radius(float radius);
@@ -121,6 +123,8 @@ private:
 	tween<float3> horizon_color_tween;
 	tween<float3> zenith_color_tween;
 	
+	math::quaternion<float> moon_rotation;
+	
 	float moon_angular_radius;
 	float cos_moon_angular_radius;
 	float sun_angular_radius;

src/resources/biome-loader.cpp

@@ -87,26 +87,6 @@ biome* resource_loader::load(resource_manager* resource_manager, PHYSFS_F
 	
 	if (auto weather = json.find("weather"); weather != json.end())
 	{
-		load_array(&biome->ambient_color.x, 3, weather.value(), "ambient_color");
-		load_value(&biome->ambient_intensity, weather.value(), "ambient_intensity");
-		
-		
-		load_value(&biome->sun_azimuth, weather.value(), "sun_azimuth");
-		biome->sun_azimuth = math::radians(biome->sun_azimuth);
-		load_value(&biome->sun_elevation, weather.value(), "sun_elevation");
-		biome->sun_elevation = math::radians(biome->sun_elevation);
-		load_array(&biome->sun_color.x, 3, weather.value(), "sun_color");
-		load_value(&biome->sun_intensity, weather.value(), "sun_intensity");
-		load_value(&biome->sun_angular_radius, weather.value(), "sun_angular_radius");
-		biome->sun_angular_radius = math::radians(biome->sun_angular_radius);
-		
-		load_array(&biome->horizon_color.x, 3, weather.value(), "horizon_color");
-		load_array(&biome->zenith_color.x, 3, weather.value(), "zenith_color");
-		
-		load_value(&biome->wind_speed, weather.value(), "wind_speed");
-		load_value(&biome->wind_direction, weather.value(), "wind_direction");
-		biome->wind_direction = math::radians(biome->wind_direction);
-		
 		std::string sky_palette_filename;
 		if (load_value(&sky_palette_filename, weather.value(), "sky_palette"))
 		{