more

http.go

@@ -29,9 +29,11 @@ func TelemRouter(log *slog.Logger, broker *JBroker) http.Handler {
 		w.Write([]byte(skylab.SkylabDefinitions))
 	})
 
+	// heartbeat request.
 	r.Get("/ping", func(w http.ResponseWriter, r *http.Request) {
 		w.Write([]byte("pong"))
 	})
+
 	r.Mount("/api/v1", apiV1(broker))
 
 	// To future residents - you can add new API calls/systems in /api/v2

py/pytelem/optimus.py

@@ -1,11 +1,12 @@
 # hyperspeed forward and backwards analytics engine
 
-import numpy.typing as npt
-from scipy.integrate import solve_bvp, solve_ivp
 
-from jax import jit, grad, vmap
+import numpy as np
 
-import jax.numpy as np
+from numba import jit
+
+
+# import jax.numpy as np
 
 
 # TODO: define 3d vector space - x,y,z oriented around car/world?
@@ -484,16 +485,18 @@ def helio_vector(vec, jme):
     )
 
 
-def solar_position(timestamp, latitude, longitude):
+def solar_position(timestamp, latitude, longitude, elevation):
     """Calculate the position of the sun at a given location and time.
 
     Args:
-        timestamp (array-like): The timestamp(s) to calculate.
+        timestamp (array-like): The timestamp(s) at each point.
+        latitude (array-like): The latitude(s) of each point.
+        longitude (array-like): The longitude(s) of each point.
+        elevation (array-like): The elevation of each point.
 
     Returns:
-        ndarray: An array containing the altitude and azimuth for each timestamp.
+        ndarray: An array containing the altitude and azimuth for each point.
     """
-    timestamp = np.array(timestamp)
     jd = timestamp / 86400.0 + 2440587.5
     jc = (jd - 2451545) / 36525
 
@@ -544,13 +547,13 @@ def solar_position(timestamp, latitude, longitude):
 
     nut = NUTATION_ABCD_ARRAY
 
-    ## TODO: these are gross - use loops instead of broadcasting?
+    # TODO: these are gross - use loops instead of broadcasting?
+    # FIXME: use guvectorize, treat jce as a scalar.
     d_psi = (nut[:, 0] + jce[..., np.newaxis] * nut[:, 1]) * np.sin(
         np.sum(X[:, np.newaxis, :] * NUTATION_YTERM_ARRAY[np.newaxis, ...], axis=2)
     )
     d_psi = np.sum(d_psi, axis=-1) / 36, 000, 000
 
-
     d_epsilon = (nut[:, 2] + jce[..., np.newaxis] * nut[:, 3]) * np.cos(
         np.sum(X[:, np.newaxis, :] * NUTATION_YTERM_ARRAY[np.newaxis, ...], axis=2)
     )
@@ -576,3 +579,38 @@ def solar_position(timestamp, latitude, longitude):
     d_tau = -20.4898 / (3600 * r_deg)
     sun_longitude = theta + d_psi + d_tau
 
+    v_0 = 280.46061837 + 360.98564736629 * (jd - 2451545) + 0.000387933 * jc ** 2 - jc ** 3 / 38710000
+    v_0 = v_0 % 360
+
+    v = v_0 + d_psi * np.cos(np.deg2rad(epsilon))
+
+    alpha = np.arctan2(np.sin(np.radians(sun_longitude)) *
+                       np.cos(np.radians(epsilon)) -
+                       np.tan(np.radians(beta)) *
+                       np.sin(np.radians(epsilon)),
+                       np.cos(np.radians(sun_longitude)))
+    alpha_deg = np.rad2deg(alpha) % 360
+    delta = np.arcsin(
+        np.sin(np.radians(beta)) *
+        np.cos(np.radians(epsilon)) +
+        np.cos(np.radians(beta)) *
+        np.sin(np.radians(epsilon)) *
+        np.cos(np.radians(sun_longitude))
+    )
+    delta_deg = np.rad2deg(delta) % 360
+
+    h = v + latitude - alpha_deg
+
+    xi_deg = 8.794 / (3600 * r_deg)
+    u = np.arctan(0.99664719 * np.tan(latitude))
+
+    x = np.cos(u) + elevation / 6378140 * np.cos(latitude)
+
+    y = 0.99664719 * np.sin(u) + elevation / 6378140 * np.sin(latitude)
+
+    d_alpha = np.arctan2(-1 * x * np.sin(np.radians(xi_deg)) * np.sin(np.radians(h)), np.cos(delta))
+    d_alpha = np.rad2deg(d_alpha)
+    alpha_prime = alpha_deg + d_alpha
+    delta_prime = np.arctan2((np.sin(delta) - y * np.sin(np.radians(xi_deg))) * np.cos(np.radians(d_alpha)),
+                             np.cos(delta) - x * np.sin(np.radians(xi_deg)) * np.cos(np.radians(h)))
+    topo_local_hour_angle_deg = h - d_alpha