Examples/pyerfa

examples/pyerfa/README.md

@@ -0,0 +1,18 @@
+# pyerfa Examples
+
+Each sub-directory contains a self-contained example. The order in
+which the examples are to appear is specified in `order.json` (an
+array of directory names in the expected order).
+
+In each example directory you'll find:
+
+* `config.toml` - must conform to the specification outlined here:
+  https://docs.pyscript.net/latest/user-guide/configuration/ This is
+  parsed and ultimately turned into a JSON representation as part of
+  the package's API object.
+* `setup.py` - Python code for contextual and environmental setup,
+  NOT SEEN BY THE END USER, but is run before the `code.py` code is
+  evaluated. Allows us to create useful (IPython) shims, avoid
+  repeating boilerplate and whatnot.
+* `code.py` - the actual code added to the editor which forms the
+  practical example of using the package.

examples/pyerfa/angles_and_separations/code.py

@@ -0,0 +1,71 @@
+# ---------------------------------------------------------------------
+# Sexagesimal parsing, angular separations, and great-circle distances.
+# ---------------------------------------------------------------------
+#
+# ERFA includes a rich set of angle utilities. Here we'll parse
+# right-ascension and declination strings, then use erfa.seps to
+# compute the angular separation between bright stars.
+import numpy as np
+import erfa
+
+
+heading("Parsing RA/Dec strings into radians")
+note(
+    "erfa.tf2a converts a sign + hours/minutes/seconds tuple to "
+    "radians (hours form), and erfa.af2a does the same for "
+    "degrees/arcminutes/arcseconds."
+)
+
+stars = [
+    # (name, RA h/m/s, Dec sign/d/m/s)
+    ("Sirius",   ( 6,  45,  8.92), ("-", 16, 42, 58.0)),
+    ("Canopus",  ( 6, 23, 57.11), ("-", 52, 41, 44.4)),
+    ("Vega",     (18, 36, 56.34), ("+", 38, 47,  1.3)),
+    ("Betelgeuse", (5, 55, 10.31), ("+",  7, 24, 25.4)),
+    ("Rigel",    ( 5, 14, 32.27), ("-",  8, 12,  6.0)),
+]
+
+# Build arrays of RA and Dec in radians.
+names = [s[0] for s in stars]
+ra_rad = np.array([erfa.tf2a("+", h, m, s) for _, (h, m, s), _ in stars])
+dec_rad = np.array([erfa.af2a(sign, d, m, s) for _, _, (sign, d, m, s) in stars])
+
+note("Star catalogue (RA, Dec converted to degrees for display):")
+rows = ["<tr><th>Name</th><th>RA (deg)</th><th>Dec (deg)</th></tr>"]
+for name, ra, dec in zip(names, np.degrees(ra_rad), np.degrees(dec_rad)):
+    rows.append(f"<tr><td>{name}</td><td>{ra:.4f}</td><td>{dec:+.4f}</td></tr>")
+display(HTML("<table>" + "".join(rows) + "</table>"), append=True)
+
+
+heading("Pairwise angular separations")
+note(
+    "erfa.seps takes (lon1, lat1, lon2, lat2) in radians and "
+    "returns the great-circle separation in radians. We broadcast "
+    "it over all pairs at once."
+)
+
+# Broadcast: compare every star against every other star.
+sep_rad = erfa.seps(
+    ra_rad[:, None], dec_rad[:, None],
+    ra_rad[None, :], dec_rad[None, :],
+)
+sep_deg = np.degrees(sep_rad)
+
+# Render as an HTML table.
+header = "<tr><th></th>" + "".join(f"<th>{n}</th>" for n in names) + "</tr>"
+body = []
+for i, name in enumerate(names):
+    cells = "".join(f"<td>{sep_deg[i, j]:6.2f}</td>" for j in range(len(names)))
+    body.append(f"<tr><th>{name}</th>{cells}</tr>")
+display(HTML("<table>" + header + "".join(body) + "</table>"), append=True)
+
+
+heading("Closest pair in the catalogue")
+# Mask the diagonal (self-separations are zero) and find the minimum.
+mask = sep_deg + np.eye(len(names)) * 1e6
+i, j = np.unravel_index(np.argmin(mask), mask.shape)
+note(
+    f"Closest pair: <strong>{names[i]}</strong> and "
+    f"<strong>{names[j]}</strong>, separated by "
+    f"<strong>{sep_deg[i, j]:.2f}°</strong>."
+)

examples/pyerfa/angles_and_separations/config.toml

@@ -0,0 +1 @@
+packages = ["pyerfa", "numpy"]

examples/pyerfa/angles_and_separations/setup.py

@@ -0,0 +1,20 @@
+"""Lightweight setup for example 3 (no IPython shim needed)."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+

examples/pyerfa/calendar_and_constants/code.py

@@ -0,0 +1,59 @@
+"""
+A first look at PyERFA: the Python wrapper for the ERFA C library
+of fundamental astronomy routines (the open-source counterpart of
+the IAU's SOFA library).
+
+Every ERFA function is exposed as a NumPy universal function, so
+they accept scalars or arrays interchangeably. See the docs at
+https://pyerfa.readthedocs.io/ for the full catalogue.
+"""
+from IPython.core.display import display, HTML
+
+import numpy as np
+import erfa
+
+
+# ERFA exposes a number of useful astronomical constants from erfam.h.
+heading("Some ERFA constants")
+note(
+    "These come straight from the ERFA C library and are handy "
+    "when working with time and angles."
+)
+constants = {
+    "DAYSEC (seconds in a day)": erfa.DAYSEC,
+    "DJY (days in a Julian year)": erfa.DJY,
+    "DAU (astronomical unit, metres)": erfa.DAU,
+    "CMPS (speed of light, m/s)": erfa.CMPS,
+    "DR2D (radians to degrees)": erfa.DR2D,
+}
+rows = "".join(f"<tr><td>{k}</td><td>{v}</td></tr>" for k, v in constants.items())
+display(HTML(f"<table>{rows}</table>"), append=True)
+
+
+# Julian Date 2460000.5 corresponds to 2023-02-25 00:00 UTC. Let's
+# convert four consecutive JDs to calendar dates with erfa.jd2cal.
+heading("Julian dates to calendar dates")
+note(
+    "erfa.jd2cal is a NumPy ufunc, so it broadcasts over array "
+    "inputs and returns aligned arrays of year, month, day, and "
+    "fractional day."
+)
+
+jd_base = 2460000.0
+jd_offsets = np.array([0, 1, 2, 3])
+year, month, day, frac = erfa.jd2cal(jd_base, jd_offsets)
+
+note(f"Year:  {year.tolist()}")
+note(f"Month: {month.tolist()}")
+note(f"Day:   {day.tolist()}")
+note(f"Frac:  {frac.tolist()}")
+
+
+# Going the other direction: calendar -> two-part Julian Date.
+heading("Calendar dates to Julian dates")
+note(
+    "erfa.cal2jd returns a (DJM0, DJM) pair where DJM0 is the MJD "
+    "zero-point (2400000.5) and DJM is the Modified Julian Date."
+)
+djm0, djm = erfa.cal2jd(2024, 7, 20)
+note(f"For 2024-07-20: DJM0 = {djm0}, MJD = {djm}, JD = {djm0 + djm}")

examples/pyerfa/calendar_and_constants/config.toml

@@ -0,0 +1 @@
+packages = ["pyerfa", "numpy"]

examples/pyerfa/calendar_and_constants/setup.py

@@ -0,0 +1,40 @@
+"""
+Shim IPython's display API onto PyScript so example code written in a
+Jupyter/IPython idiom runs unmodified in the browser.
+"""
+
+import sys
+import types
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+ipython = types.ModuleType("IPython")
+core = types.ModuleType("IPython.core")
+core_display = types.ModuleType("IPython.core.display")
+core_display.display = display
+core_display.HTML = HTML
+ipython.core = core
+core.display = core_display
+ipython.get_ipython = lambda: None
+ipython.display = core_display
+sys.modules["IPython"] = ipython
+sys.modules["IPython.core"] = core
+sys.modules["IPython.core.display"] = core_display
+sys.modules["IPython.display"] = core_display
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)

examples/pyerfa/order.json

@@ -0,0 +1,5 @@
+[
+    "calendar_and_constants",
+    "planet_positions",
+    "angles_and_separations"
+]

examples/pyerfa/planet_positions/code.py

@@ -0,0 +1,61 @@
+# ---------------------------------------------------------------------
+# Plotting the inner planets' orbits using erfa.plan94.
+# ---------------------------------------------------------------------
+#
+# erfa.plan94 gives heliocentric position and velocity (in AU and
+# AU/day) for planets 1..8 (Mercury through Neptune) using an
+# analytical approximation accurate enough for many purposes.
+
+import numpy as np
+import matplotlib.pyplot as plt
+import erfa
+
+
+heading("Tracing the inner planets over one Earth-year")
+note(
+    "We sample plan94 every few days for a full year starting "
+    "on JD 2460000.5 and project the heliocentric positions onto "
+    "the ecliptic plane (X, Y in AU)."
+)
+
+planets = {
+    1: ("Mercury", "tab:gray"),
+    2: ("Venus",   "tab:orange"),
+    3: ("Earth",   "tab:blue"),
+    4: ("Mars",    "tab:red"),
+}
+
+jd_base = 2460000.5
+days = np.arange(0, 700, 2.0)  # enough to close Mars' orbit too
+
+fig, ax = plt.subplots(figsize=(6, 6))
+ax.plot(0, 0, marker="*", color="gold", markersize=18, label="Sun")
+
+for planet_id, (name, color) in planets.items():
+    pv = erfa.plan94(jd_base, days, planet_id)
+    # pv is a structured array with fields 'p' (position) and 'v' (velocity).
+    x = pv["p"][:, 0]
+    y = pv["p"][:, 1]
+    ax.plot(x, y, color=color, linewidth=1.2, label=name)
+    # Mark the starting position.
+    ax.plot(x[0], y[0], "o", color=color, markersize=5)
+
+ax.set_aspect("equal")
+ax.set_xlabel("X (AU)")
+ax.set_ylabel("Y (AU)")
+ax.set_title("Inner planets, heliocentric ecliptic plane")
+ax.legend(loc="upper right", fontsize=9)
+ax.grid(alpha=0.3)
+fig.tight_layout()
+display(fig, append=True)
+
+# How far is each planet from the Sun on day zero?
+heading("Heliocentric distances on the start date")
+pv0 = erfa.plan94(jd_base, 0.0, 3)  # Earth as a check: should be ~1 AU
+note(f"Earth's distance from the Sun on JD {jd_base}: "
+     f"<strong>{np.linalg.norm(pv0['p']):.4f} AU</strong>")
+
+for planet_id, (name, _) in planets.items():
+    pv = erfa.plan94(jd_base, 0.0, planet_id)
+    r = np.linalg.norm(pv["p"])
+    note(f"{name}: {r:.4f} AU")

examples/pyerfa/planet_positions/config.toml

@@ -0,0 +1 @@
+packages = ["pyerfa", "numpy", "matplotlib"]

examples/pyerfa/planet_positions/setup.py

@@ -0,0 +1,19 @@
+"""Lightweight setup for example 2 (no IPython shim needed)."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)