Modeling Earth

Modeling Earth

Model the earth from horizontal coordinates with day/night cycle with stars and sun modeled, add atmosphere later.
Integrate various kinds of earth datasets.
Visualizations.
Spaceflight simulator.

Celestial coordinate systems

Spherical coordinates

import math

def sphericalcoord(r, theta, phi): # altitude, longitude, latitude
    x = r * math.sin(theta) * math.cos(phi)
    y = r * math.sin(theta) * math.sin(phi)
    z = r * math.cos(theta)
    return sphere_x, sphere_y, sphere_z

def cartesian_to_spherical(x, y, z): # z is up
    r = math.sqrt(math.pow(x, 2) + math.pow(y, 2) + math.pow(z, 2))
    theta = math.acos(z / r)
    phi = math.atan(y / x)
    return r, theta, phi

Using Numpy implementing from The sun's position in the sky.

import numpy as np
import math

spherecoord = np.array([r * math.sin(theta) * math.cos(phi),
                        r * math.sin(theta) * math.sin(phi), 
                        r * math.cos(theta)])

Ecliptic

$$ \vec{r}_{\odot}(t) = \begin{bmatrix}\cos{\phi_{\odot}}(t) \\ \sin{\phi_{\odot}}(t) \\ 0\end{bmatrix} $$

t = 0 # time

ecliptic = np.array([math.cos(t), math.sin(t), 0])
ecliptic.reshape(3, 1)
ecliptic = np.matrix(ecliptic)

Equatorial

$$ \varepsilon = 23.44 = 0.40910517666747087 $$

$$ \begin{bmatrix}1 & 0 & 0 \\ 0 & \cos{\varepsilon} & -\sin{\varepsilon} \\ 0 & \sin{\varepsilon} & \cos{\varepsilon}\end{bmatrix} $$

epsilon = math.radians(23.44) # earths tilt in radians

axialtilt = np.array([1,                 0,                 0,
                      0, math.cos(epsilon), math.sin(epsilon),
                      0, math.sin(epsilon), math.cos(epsilon)])
axialtilt.reshape(3, 3)
axialtilt = np.matrix(axialtilt)

equatorial = axialtilt * ecliptic

$$ \vec{r}'_{\odot} = \underbrace{ \begin{bmatrix}1 & 0 & 0 \\ 0 & \cos{\varepsilon} & -\sin{\varepsilon} \\ 0 & \sin{\varepsilon} & \cos{\varepsilon}\end{bmatrix} }_{\text{Equatorial}} \underbrace{ \begin{bmatrix}\cos{\phi_{\odot}} \\ \sin{\phi_{\odot}} \\ 0\end{bmatrix} }_{\text{Ecliptic}} = \begin{bmatrix} \cos{\phi_\odot} \\ \cos{\varepsilon} \sin{\phi_\odot} \\ \sin{\varepsilon} \sin{\phi_\odot}\end{bmatrix} $$

Terrestrial frame

Sidereal day

A sidereal day is the time it takes for earth to make one revolution around it's axis of rotation so the stars line up in exactly the same place they started in the sky from one sidereal day ago.

$$ \omega = \frac{2 \pi}{23.9345 \text{ hours}} $$

$$ \vec{r}''_{\odot} = \begin{bmatrix} x''_{\odot} \\ y''_{\odot} \\ z''_{\odot} \end{bmatrix} = \begin{bmatrix}\sin{\theta''_{\odot}} \cos{\phi''_{\odot}} \\ \sin{\theta''_{\odot}} \sin{\phi''_{\odot}} \\ \cos{\theta''_{\odot}}\end{bmatrix} = \underbrace{ \begin{bmatrix}1 & 0 & 0 \\ 0 & \sin{L} & -\cos{L} \\ 0 & \cos{L} & \sin{L}\end{bmatrix} }_{\text{Latitude / co-latitude rotation}} \underbrace{ \begin{bmatrix}\cos{[\omega(t - t_0)]} & \sin{[\omega(t - t_0)]} & 0 \\ \sin{[\omega(t - t_0)]} & \cos{[\omega(t - t_0)]} & 0 \\ 0 & 0 & 1\end{bmatrix} }_{\text{Rotation of Earth}} \begin{bmatrix} \cos{\phi_\odot} \\ \cos{\varepsilon} \sin{\phi_\odot} \\ \sin{\varepsilon} \sin{\phi_\odot}\end{bmatrix} $$

Using set-builder notation we define the interval for $ L $ as

$$ \{ L \in \mathbb{R} \mid -\frac{\pi}{2} \leq L \leq \frac{\pi}{2} \} $$

We can also define the interval with this notation

$$ L \in \left[ -\frac{\pi}{2},\ \frac{\pi}{2} \right] $$

$$ \theta'_{\odot} = \frac{\pi}{2} - L $$

Place	Latitude	Colatitude
North pole	$\frac{\pi}{2}$	$0$
Equator	$0$	$\frac{\pi}{2}$
South pole	$-\frac{\pi}{2}$	$\pi$

$L = 0$: Equator
$L > 0$: Northern Hemisphere (up to the North Pole)
$L < 0$: Southern Hemisphere (down to the South Pole)
$L = +\frac{\pi}{2}$: North Pole
$L = -\frac{\pi}{2}$: South Pole

latitude = 0 # let's place the observation point at the equator for now

terrestrial = np.array([1,                  0,                   0,
                        0, math.sin(latitude), -math.cos(latitude),
                        0, math.cos(latitude), math.sin(latitude)])
terrestrial.reshape(3, 3)
terrestrial = np.matrix(terrestrial)

Zenith - Wikipedia
Azimuth - Wikipedia
Latitude Colatitude - Wikipedia
Horizontal coordinate system - Wikipedia
Geographic coordinate system Axial precession - Wikipedia Equation of the center - Wikipedia Kepler's equation Equation of time Apsis Celestial mechanics

Star data

We can generate a skybox from measured data.

Hipparcos - Wikipedia

The Hipparcos and Tycho Catalogues (1997)

Modular arithmetic

Modulo etc for timekeeping.

Barycentric coordinates/interpolation

Interpolate between three data points.

$$ \vec{s}=u\vec{A}+v\vec{B}+w\vec{C}=u\vec{A}+v\vec{B}+(1-(u+v))\vec{C} $$

Sky models

Sky models - Github A Framework for the Experimental Comparison of Solar and Skydome Illumination

Sources

https://iopscience.iop.org/article/10.1088/0031-9120/45/6/010

GPT DUMP

Ah, perfect — you're likely referring to the axial tilt matrix that transforms the Sun’s position from the ecliptic frame to the equatorial frame.

That’s a key step!

🌍 Where It Fits

The full transformation pipeline should be:

Sun in ecliptic coordinates
Rotate by Earth's axial tilt $\varepsilon$
Rotate by Earth’s rotation angle $\omega(t - t_0)$
Rotate by co-latitude $\frac{\pi}{2} - L$ (or use the matrix that directly encodes this)
→ Result: Sun in terrestrial frame (as seen by observer)

✅ Full Matrix Equation (Including Tilt)

$$ \vec{r}''_{\odot} = \underbrace{ \begin{bmatrix} 1 & 0 & 0 \\ 0 & \sin{L} & -\cos{L} \\ 0 & \cos{L} & \sin{L} \end{bmatrix} }_{\text{Latitude / co-latitude rotation}} \cdot \underbrace{ \begin{bmatrix} \cos{\theta} & \sin{\theta} & 0 \\ -\sin{\theta} & \cos{\theta} & 0 \\ 0 & 0 & 1 \end{bmatrix} }_{\text{Earth rotation, } \theta = \omega(t - t_0)} \cdot \underbrace{ \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos{\varepsilon} & -\sin{\varepsilon} \\ 0 & \sin{\varepsilon} & \cos{\varepsilon} \end{bmatrix} }_{\text{Axial tilt (ecliptic → equatorial)}} \cdot \underbrace{ \begin{bmatrix} \cos{\phi_\odot} \\ \sin{\phi_\odot} \\ 0 \end{bmatrix} }_{\text{Sun in ecliptic frame}} $$

🧭 Summary of Each Matrix

Matrix	Description	Axis of Rotation
Axial tilt $R_{\varepsilon}$	Ecliptic → Equatorial	x-axis
Earth rotation $R_z(\theta)$	Equatorial → rotating Earth	z-axis
Observer tilt $R_x(\pi/2 - L)$	Earth → local horizontal	x-axis (tilts for latitude)

This full composition gives the Sun’s position as seen from a point on Earth at latitude $L$ and time $t$, ready for rendering.

Would you like me to re-write this into code or Blender matrices again with all four steps?

Absolutely!

🌍 Interval Notation for Geographic Latitude $L$

In radians, the geographic latitude $L$ lies within the following range:

$$ L \in \left[ -\frac{\pi}{2},\ \frac{\pi}{2} \right] $$

Or approximately:

$$ L \in [-1.5708,\ 1.5708] $$

🧭 Interpretation:

$L = 0$: Equator
$L > 0$: Northern Hemisphere (up to the North Pole)
$L < 0$: Southern Hemisphere (down to the South Pole)
$L = +\frac{\pi}{2}$: North Pole
$L = -\frac{\pi}{2}$: South Pole

Let me know if you also want this in degrees or want to work with co-latitude instead!

Ah, got it — you're thinking of set-builder notation!

🧮 Set-Builder Notation for Latitude $L$

To express latitude $L$ using curly braces and a condition:

$$ \{ L \in \mathbb{R} \mid -\frac{\pi}{2} \leq L \leq \frac{\pi}{2} \} $$

✅ Explanation

${ \cdot }$ = “the set of all…”
$L \in \mathbb{R}$ = “L is a real number”
$\mid$ = “such that”
The inequality defines the range of valid latitudes in radians.

📖 Name of this Notation

It's called set-builder notation, and it's used to describe a set by specifying a property its members must satisfy.

Let me know if you want the same for co-latitude or for degrees!

Place	Latitude	Colatitude
North pole	\(\frac{\pi}{2}\)	\(0\)
Equator	\(0\)	\(\frac{\pi}{2}\)
South pole	\(-\frac{\pi}{2}\)	\(\pi\)