diff --git a/examples/pyerfa/README.md b/examples/pyerfa/README.md new file mode 100644 index 0000000..7d18dac --- /dev/null +++ b/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. diff --git a/examples/pyerfa/angles_and_separations/code.py b/examples/pyerfa/angles_and_separations/code.py new file mode 100644 index 0000000..31d8010 --- /dev/null +++ b/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 = ["NameRA (deg)Dec (deg)"] +for name, ra, dec in zip(names, np.degrees(ra_rad), np.degrees(dec_rad)): + rows.append(f"{name}{ra:.4f}{dec:+.4f}") +display(HTML("" + "".join(rows) + "
"), 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 = "" + "".join(f"{n}" for n in names) + "" +body = [] +for i, name in enumerate(names): + cells = "".join(f"{sep_deg[i, j]:6.2f}" for j in range(len(names))) + body.append(f"{name}{cells}") +display(HTML("" + header + "".join(body) + "
"), 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: {names[i]} and " + f"{names[j]}, separated by " + f"{sep_deg[i, j]:.2f}°." +) diff --git a/examples/pyerfa/angles_and_separations/config.toml b/examples/pyerfa/angles_and_separations/config.toml new file mode 100644 index 0000000..f2f0312 --- /dev/null +++ b/examples/pyerfa/angles_and_separations/config.toml @@ -0,0 +1 @@ +packages = ["pyerfa", "numpy"] diff --git a/examples/pyerfa/angles_and_separations/setup.py b/examples/pyerfa/angles_and_separations/setup.py new file mode 100644 index 0000000..2d23644 --- /dev/null +++ b/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"{text}"), append=True) + + +def note(text): + display(HTML(f"

{text}

"), append=True) + diff --git a/examples/pyerfa/calendar_and_constants/code.py b/examples/pyerfa/calendar_and_constants/code.py new file mode 100644 index 0000000..ef8cee4 --- /dev/null +++ b/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"{k}{v}" for k, v in constants.items()) +display(HTML(f"{rows}
"), 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}") diff --git a/examples/pyerfa/calendar_and_constants/config.toml b/examples/pyerfa/calendar_and_constants/config.toml new file mode 100644 index 0000000..f2f0312 --- /dev/null +++ b/examples/pyerfa/calendar_and_constants/config.toml @@ -0,0 +1 @@ +packages = ["pyerfa", "numpy"] diff --git a/examples/pyerfa/calendar_and_constants/setup.py b/examples/pyerfa/calendar_and_constants/setup.py new file mode 100644 index 0000000..84faac4 --- /dev/null +++ b/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"{text}"), append=True) + + +def note(text): + display(HTML(f"

{text}

"), append=True) diff --git a/examples/pyerfa/order.json b/examples/pyerfa/order.json new file mode 100644 index 0000000..2d70c9f --- /dev/null +++ b/examples/pyerfa/order.json @@ -0,0 +1,5 @@ +[ + "calendar_and_constants", + "planet_positions", + "angles_and_separations" +] diff --git a/examples/pyerfa/planet_positions/code.py b/examples/pyerfa/planet_positions/code.py new file mode 100644 index 0000000..b0b48bb --- /dev/null +++ b/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"{np.linalg.norm(pv0['p']):.4f} AU") + +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") diff --git a/examples/pyerfa/planet_positions/config.toml b/examples/pyerfa/planet_positions/config.toml new file mode 100644 index 0000000..1cad054 --- /dev/null +++ b/examples/pyerfa/planet_positions/config.toml @@ -0,0 +1 @@ +packages = ["pyerfa", "numpy", "matplotlib"] diff --git a/examples/pyerfa/planet_positions/setup.py b/examples/pyerfa/planet_positions/setup.py new file mode 100644 index 0000000..4b05c08 --- /dev/null +++ b/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"{text}"), append=True) + + +def note(text): + display(HTML(f"

{text}

"), append=True)