Open In Colab

1cf2dc2c7ed94ce7ad57d029dff4bc21

The Heliacal Rise of Sirius, Step by Step

This notebook reconstructs the first morning visibility of Sirius from Thebes around 2782 BCE, the era of the first apokatastasis of the Sothic cycle.

The goal is not to treat montu.HeliacalRise as a black box. For each of the four models included in MontuPython we will:

  1. obtain the positions of the Sun and Sirius using MontuPython’s public positional astronomy routines;

  2. write and explain the visibility criterion;

  3. evaluate it morning by morning —and, when appropriate, minute by minute—;

  4. identify the «not visible → visible» transition that defines the heliacal rise;

  5. compare the manual reconstruction with the high-level API.

Important interpretation. A calculated heliacal rise is not a purely geometric instant like a rise or a conjunction. It is a visibility prediction and depends on the model, the atmosphere, the horizon, the observer’s acuity, and the adopted time step.

If running in Google Colab, MontuPython must be installed first. In a local copy of the repository this cell can remain commented out.

[1]:

try:
    from google.colab import drive
    %pip install -Uq montu
except ImportError:
    print("Not running in Colab, skipping installation")
    import plotly.io as pio
    pio.renderers.default = "notebook_connected"
    %load_ext autoreload
    %autoreload 2
# Create folders for figures and temporal files
!mkdir -p ./gallery/ ./montu_dem/

Not running in Colab, skipping installation

[2]:

%matplotlib inline
import matplotlib.pyplot as plt
import montu
import numpy as np
import pandas as pd

pd.options.display.float_format = "{:.3f}".format

MontuPython version 0.50.0. 𓇍𓇋𓇋𓏏𓅓𓊵 𓎛𓎡𓄿𓀭𓎛𓈖𓂝𓎡 (ii-ti m Htp, HkAx Hn'-k)

1. Setup, conventions, and input data

We need to fix three elements before discussing visibility:

  • Observer: Thebes (Luxor), from the MontuPython site catalogue.

  • Object: Sirius, from the bright-star catalogue. Its catalogue visual magnitude is approximately \(V=-1.44\).

  • Interval: July 2782 BCE.

MontuPython uses astronomical year numbering: the historical year 2782 BCE is written as -2781, because astronomical year 0 exists. We use calendar='mixed', i.e. the Julian calendar before the Gregorian reform.

All altitudes are apparent and expressed in degrees. Azimuth is measured from north toward east (\(0°\) north, \(90°\) east). JED is the Julian day on the UTC scale that MontuPython uses as a time coordinate.

[3]:

site = montu.Observer(site="thebes")
sirius = montu.Stars(subset="bright", ProperName="Sirius")
sun = montu.Sun()

start_date = montu.Time("139-07-01", calendar="mixed")
end_date = montu.Time("139-07-30", calendar="mixed")
V_SIRIUS = float(sirius.data.iloc[0].Vmag)

start_date, end_date

Loading stellar catalogue montu_stellar_catalogue_v38_bright.csv

[3]:

(Time('139-07-01 02:33:37.2'/'139-07-01 00:00:00'/'[hrw 2921] IV shemu 17'/JED 1772008.5/JTD 1772008.6066806),
 Time('139-07-30 02:33:37.2'/'139-07-30 00:00:00'/'[hrw 2922] I akhet 11'/JED 1772037.5/JTD 1772037.6066806))

2. Models available in MontuPython

MontuPython provides a set of routines to estimate the heliacal rise of a star or planet using four models:

  • Geometric models:

    • Ptolemy model (``ptolemy``): This is the classical model from Ptolemy’s Almagest. The critical quantity for morning visibility is the Arcus Visionis: the difference between the elevations of the Sun and the observed object above a given horizon (see below) at the exact moment the object rises. This is the most widely used model in archaeoastronomy.

  • Physical models

    • Schaefer basic model (``schaefer1987``): This model computes heliacal-rise conditions with atmospheric extinction. Essentially, it requires the sky to be dark enough at dawn — defined by the Sun’s depression below the horizon and the limiting visual magnitude at the object’s altitude.

    • Schaefer physiological model (``schaefer1985``): Similar to the basic model by the same author, but it includes physiological factors of observation for a more realistic estimate. The model has been tested against real observations.

    • Belokrylov et al. model (``belokrylov2011``): This model uses elements of Schaefer’s basic approach, but applies corrections derived from real observations for objects of different magnitudes. It is a semi-empirical model and one of the most recent revisions of heliacal-rise theory.

3. The HeliacalRise interface in MontuPython

The details of all models are implemented in MontuPython in a class, HeliacalRise, which computes heliacal rises in a user-friendly way without going deeply into the astronomical and physical details.

3.1 Ptolemy’s method (Arcus Visionis)

The class is used simply by defining the model and its parameters:

[4]:

model = montu.HeliacalRise(
    model="ptolemy",
    #arcus_visionis_crit=14.0,  # assumed Arcus Visionis
    #ptolemy_refraction_deg=0.0/60  # astronomical refraction (if desired)
)

Once the class is defined, we compute the conditions over the chosen date interval:

[5]:

rises = model.compute(sirius, site, start_date, end_date, verbose=True)
rises

HeliacalRise verbose — model=ptolemy
  quantities: t_rise: local object rise time; AV: solar depression at object rise; AV_crit: critical Arcus Visionis; h_sun: solar altitude
  interval: 139-07-01 00:00:00 -> 139-07-30 00:00:00
  model parameters:
    arcus_visionis_crit=automatic body default
    ptolemy_refraction_deg=0.5666666666666667
  criterion: AV_calc >= AV_crit and h_sun < 0°
  day 001 | 139-07-01 00:00:00 | visible=False | t_rise=05:22:27.403 | AV=-2.348° | AV_crit=14.000° | h_sun=2.348°
  day 002 | 139-07-02 00:00:00 | visible=False | t_rise=05:18:31.479 | AV=-1.480° | AV_crit=14.000° | h_sun=1.480°
  day 003 | 139-07-03 00:00:00 | visible=False | t_rise=05:14:35.564 | AV=-0.611° | AV_crit=14.000° | h_sun=0.611°
  day 004 | 139-07-04 00:00:00 | visible=False | t_rise=05:10:39.658 | AV=0.259° | AV_crit=14.000° | h_sun=-0.259°
  day 005 | 139-07-05 00:00:00 | visible=False | t_rise=05:06:43.751 | AV=1.128° | AV_crit=14.000° | h_sun=-1.128°
  day 006 | 139-07-06 00:00:00 | visible=False | t_rise=05:02:47.844 | AV=1.998° | AV_crit=14.000° | h_sun=-1.998°
  day 007 | 139-07-07 00:00:00 | visible=False | t_rise=04:58:51.938 | AV=2.868° | AV_crit=14.000° | h_sun=-2.868°
  day 008 | 139-07-08 00:00:00 | visible=False | t_rise=04:54:56.031 | AV=3.739° | AV_crit=14.000° | h_sun=-3.739°
  day 009 | 139-07-09 00:00:00 | visible=False | t_rise=04:51:00.125 | AV=4.609° | AV_crit=14.000° | h_sun=-4.609°
  day 010 | 139-07-10 00:00:00 | visible=False | t_rise=04:47:04.210 | AV=5.480° | AV_crit=14.000° | h_sun=-5.480°
  day 011 | 139-07-11 00:00:00 | visible=False | t_rise=04:43:08.294 | AV=6.350° | AV_crit=14.000° | h_sun=-6.350°
  day 012 | 139-07-12 00:00:00 | visible=False | t_rise=04:39:12.379 | AV=7.221° | AV_crit=14.000° | h_sun=-7.221°
  day 013 | 139-07-13 00:00:00 | visible=False | t_rise=04:35:16.464 | AV=8.092° | AV_crit=14.000° | h_sun=-8.092°
  day 014 | 139-07-14 00:00:00 | visible=False | t_rise=04:31:20.549 | AV=8.962° | AV_crit=14.000° | h_sun=-8.962°
  day 015 | 139-07-15 00:00:00 | visible=False | t_rise=04:27:24.634 | AV=9.833° | AV_crit=14.000° | h_sun=-9.833°
  day 016 | 139-07-16 00:00:00 | visible=False | t_rise=04:23:28.727 | AV=10.703° | AV_crit=14.000° | h_sun=-10.703°
  day 017 | 139-07-17 00:00:00 | visible=False | t_rise=04:19:32.812 | AV=11.573° | AV_crit=14.000° | h_sun=-11.573°
  day 018 | 139-07-18 00:00:00 | visible=False | t_rise=04:15:36.905 | AV=12.443° | AV_crit=14.000° | h_sun=-12.443°
  day 019 | 139-07-19 00:00:00 | visible=False | t_rise=04:11:40.999 | AV=13.313° | AV_crit=14.000° | h_sun=-13.313°
  day 020 | 139-07-20 00:00:00 | visible=True | t_rise=04:07:45.092 | AV=14.182° | AV_crit=14.000° | h_sun=-14.182°
    -> heliacal rise detected
       criterion satisfied: AV_calc >= AV_crit and h_sun < 0°
  day 021 | 139-07-21 00:00:00 | visible=True | t_rise=04:03:49.194 | AV=15.050° | AV_crit=14.000° | h_sun=-15.050°
  day 022 | 139-07-22 00:00:00 | visible=True | t_rise=03:59:53.288 | AV=15.919° | AV_crit=14.000° | h_sun=-15.919°
  day 023 | 139-07-23 00:00:00 | visible=True | t_rise=03:55:57.381 | AV=16.786° | AV_crit=14.000° | h_sun=-16.786°
  day 024 | 139-07-24 00:00:00 | visible=True | t_rise=03:52:01.466 | AV=17.653° | AV_crit=14.000° | h_sun=-17.653°
  day 025 | 139-07-25 00:00:00 | visible=True | t_rise=03:48:05.559 | AV=18.520° | AV_crit=14.000° | h_sun=-18.520°
  day 026 | 139-07-26 00:00:00 | visible=True | t_rise=03:44:09.644 | AV=19.386° | AV_crit=14.000° | h_sun=-19.386°
  day 027 | 139-07-27 00:00:00 | visible=True | t_rise=03:40:13.729 | AV=20.251° | AV_crit=14.000° | h_sun=-20.251°
  day 028 | 139-07-28 00:00:00 | visible=True | t_rise=03:36:17.805 | AV=21.115° | AV_crit=14.000° | h_sun=-21.115°
  day 029 | 139-07-29 00:00:00 | visible=True | t_rise=03:32:21.890 | AV=21.978° | AV_crit=14.000° | h_sun=-21.978°
  day 030 | 139-07-30 00:00:00 | visible=True | t_rise=03:28:25.975 | AV=22.841° | AV_crit=14.000° | h_sun=-22.841°

[5]:
model source day_jed jed local_time body_altitude_deg body_azimuth_deg sun_altitude_deg sun_azimuth_deg sun_altitude_formula_deg vmag target_horizon_deg body_ra_hours body_dec_deg sun_ra_hours sun_dec_deg h_star_deg h_sun_deg arcus_visionis_calc_deg arcus_visionis_crit_deg
0 ptolemy Toomer, G. J. (1998). Ptolemy's Almagest. Prin... 1772027.500 1772027.581 04:07:45.092 -0.567 107.392 -14.182 57.377 -14.182 -1.440 -0.567 5.384 -15.880 7.780 21.392 -82.792 -118.719 14.182 14.000

As you can see, the algorithm returns all detected heliacal rises in the interval as a DataFrame. To inspect the date:

[6]:

montu.Time(float(rises.jed[0]), format="jd", calendar="mixed")

[6]:

Time('139-07-20 04:30:48.2'/'139-07-20 01:57:57'/'[hrw 2922] I akhet 1'/JED 1772027.5813773/JTD 1772027.6880579)

The package includes a method to present the result more clearly:

[7]:

model.print_rises(rises)

ptolemy — 1 date(s)
  [1] 139-07-20 00:00:00  04:07:45.092  0139-07-19 00:00:00.000000  [hrw 2922] I akhet 1  body -0.57°  Sun -14.18°
  source: Toomer, G. J. (1998). Ptolemy's Almagest. Princeton University Press. Book XIII, Chapter 7: "On the heliacal risings and settings of the planets".

Or, if you prefer something a bit more polished:

[8]:

model.print_rises(rises, title="Heliacal rise of Sirius — Ptolemy", body_label="Sirius")

Heliacal rise of Sirius — Ptolemy — 1 date(s)
  [1] 139-07-20 00:00:00  04:07:45.092  0139-07-19 00:00:00.000000  [hrw 2922] I akhet 1  Sirius -0.57°  Sun -14.18°
  source: Toomer, G. J. (1998). Ptolemy's Almagest. Princeton University Press. Book XIII, Chapter 7: "On the heliacal risings and settings of the planets".

Save the result:

[9]:

rise_ptolemy_jed = rises.day_jed[0]

3.2 Schaefer (1987)

The same approach as for the other models, but with a different set of input parameters.

[10]:

model = montu.HeliacalRise(
    model="schaefer1987",
    k=0.25,  # atmospheric extinction
    limiting_mag_zenith=6.0,  # limiting visual magnitude at the zenith
    sun_depression=-11.0,  # solar depression required for observation
)
rises = model.compute(sirius, site, start_date, end_date, verbose=True)
model.print_rises(rises, title="Heliacal rise of Sirius — Schaefer 1987", body_label="Sirius")

HeliacalRise verbose — model=schaefer1987
  quantities: h_star: object altitude; X: air mass; V_observed: extinguished object magnitude; V_limit: limiting magnitude at object altitude
  interval: 139-07-01 00:00:00 -> 139-07-30 00:00:00
  model parameters:
    k=0.25
    limiting_mag_zenith=6.0
    sun_depression=-11.0
  criterion: h_star > 0° and V_observed <= V_limit(local)
  day 001 | 139-07-01 00:00:00 | visible=False | h_star=-15.410° | X=inf | V_observed=inf | V_limit=-inf
  day 002 | 139-07-02 00:00:00 | visible=False | h_star=-14.455° | X=inf | V_observed=inf | V_limit=-inf
  day 003 | 139-07-03 00:00:00 | visible=False | h_star=-13.498° | X=inf | V_observed=inf | V_limit=-inf
  day 004 | 139-07-04 00:00:00 | visible=False | h_star=-12.538° | X=inf | V_observed=inf | V_limit=-inf
  day 005 | 139-07-05 00:00:00 | visible=False | h_star=-11.576° | X=inf | V_observed=inf | V_limit=-inf
  day 006 | 139-07-06 00:00:00 | visible=False | h_star=-10.611° | X=inf | V_observed=inf | V_limit=-inf
  day 007 | 139-07-07 00:00:00 | visible=False | h_star=-9.645° | X=inf | V_observed=inf | V_limit=-inf
  day 008 | 139-07-08 00:00:00 | visible=False | h_star=-8.677° | X=inf | V_observed=inf | V_limit=-inf
  day 009 | 139-07-09 00:00:00 | visible=False | h_star=-7.707° | X=inf | V_observed=inf | V_limit=-inf
  day 010 | 139-07-10 00:00:00 | visible=False | h_star=-6.737° | X=inf | V_observed=inf | V_limit=-inf
  day 011 | 139-07-11 00:00:00 | visible=False | h_star=-5.765° | X=inf | V_observed=inf | V_limit=-inf
  day 012 | 139-07-12 00:00:00 | visible=False | h_star=-4.792° | X=inf | V_observed=inf | V_limit=-inf
  day 013 | 139-07-13 00:00:00 | visible=False | h_star=-3.819° | X=inf | V_observed=inf | V_limit=-inf
  day 014 | 139-07-14 00:00:00 | visible=False | h_star=-2.845° | X=inf | V_observed=inf | V_limit=-inf
  day 015 | 139-07-15 00:00:00 | visible=False | h_star=-1.872° | X=inf | V_observed=inf | V_limit=-inf
  day 016 | 139-07-16 00:00:00 | visible=False | h_star=-0.898° | X=inf | V_observed=inf | V_limit=-inf
  day 017 | 139-07-17 00:00:00 | visible=False | h_star=0.074° | X=770.605 | V_observed=191.211 | V_limit=-186.401
  day 018 | 139-07-18 00:00:00 | visible=False | h_star=1.046° | X=54.757 | V_observed=12.249 | V_limit=-7.439
  day 019 | 139-07-19 00:00:00 | visible=False | h_star=2.017° | X=28.405 | V_observed=5.661 | V_limit=-0.851
  day 020 | 139-07-20 00:00:00 | visible=False | h_star=2.987° | X=19.188 | V_observed=3.357 | V_limit=1.453
  day 021 | 139-07-21 00:00:00 | visible=True | h_star=3.956° | X=14.496 | V_observed=2.184 | V_limit=2.626
    -> heliacal rise detected
       criterion satisfied: h_star > 0° and V_observed <= V_limit(local)
  day 022 | 139-07-22 00:00:00 | visible=True | h_star=4.922° | X=11.655 | V_observed=1.474 | V_limit=3.336
  day 023 | 139-07-23 00:00:00 | visible=True | h_star=5.887° | X=9.750 | V_observed=0.998 | V_limit=3.812
  day 024 | 139-07-24 00:00:00 | visible=True | h_star=6.849° | X=8.386 | V_observed=0.656 | V_limit=4.154
  day 025 | 139-07-25 00:00:00 | visible=True | h_star=7.809° | X=7.360 | V_observed=0.400 | V_limit=4.410
  day 026 | 139-07-26 00:00:00 | visible=True | h_star=8.766° | X=6.562 | V_observed=0.200 | V_limit=4.610
  day 027 | 139-07-27 00:00:00 | visible=True | h_star=9.720° | X=5.923 | V_observed=0.041 | V_limit=4.769
  day 028 | 139-07-28 00:00:00 | visible=True | h_star=10.670° | X=5.401 | V_observed=-0.090 | V_limit=4.900
  day 029 | 139-07-29 00:00:00 | visible=True | h_star=11.617° | X=4.966 | V_observed=-0.199 | V_limit=5.009
  day 030 | 139-07-30 00:00:00 | visible=True | h_star=12.560° | X=4.598 | V_observed=-0.290 | V_limit=5.100
Heliacal rise of Sirius — Schaefer 1987 — 1 date(s)
  [1] 139-07-21 00:00:00  04:24:56.311  0139-07-20 00:00:00.000000  [hrw 2922] I akhet 2  Sirius 3.96°  Sun -11.00°
  source: Schaefer, B. E. (1987). "Heliacal rise phenomena". Journal for the History of Astronomy, 18(11), 19-33.

Save the result:

[11]:

rise_schaefer1987_jed = rises.day_jed[0]

3.3 Schaefer (1985)

This model scans the twilight. step_minutes=2 means that visibility is checked every two minutes between astronomical dawn and sunrise.

[12]:

model = montu.HeliacalRise(
    model="schaefer1985",
    k=0.25,  # atmospheric extinction coefficient
    limiting_mag_zenith=6.0,  # limiting visual magnitude at the zenith
    step_minutes=2.0,  # time step
    twilight_sunbelow=-18.0,  # solar depression
)
rises = model.compute(sirius, site, start_date, end_date, verbose=True)
model.print_rises(rises, title="Heliacal rise of Sirius — Schaefer 1985", body_label="Sirius")

HeliacalRise verbose — model=schaefer1985
  quantities: h_star: object altitude; h_sun: solar altitude; V: object magnitude; V_limit: limiting magnitude at object altitude; B: sky brightness
  interval: 139-07-01 00:00:00 -> 139-07-30 00:00:00
  model parameters:
    k=0.25
    limiting_mag_zenith=6.0
    step_minutes=2.0
    twilight_sunbelow=-18.0
  criterion: h_star > 0° and V <= V_limit
  morning scan: h_sun=-18.0° to sunrise, every 2.0 minutes
  day 001 | 139-07-01 00:00:00 | visible=False
  day 002 | 139-07-02 00:00:00 | visible=False
  day 003 | 139-07-03 00:00:00 | visible=False
  day 004 | 139-07-04 00:00:00 | visible=False
  day 005 | 139-07-05 00:00:00 | visible=False
  day 006 | 139-07-06 00:00:00 | visible=False
  day 007 | 139-07-07 00:00:00 | visible=False
  day 008 | 139-07-08 00:00:00 | visible=False
  day 009 | 139-07-09 00:00:00 | visible=False
  day 010 | 139-07-10 00:00:00 | visible=False
  day 011 | 139-07-11 00:00:00 | visible=False
  day 012 | 139-07-12 00:00:00 | visible=False
  day 013 | 139-07-13 00:00:00 | visible=False
  day 014 | 139-07-14 00:00:00 | visible=False
  day 015 | 139-07-15 00:00:00 | visible=True | h_star=2.539° | h_sun=-6.949° | V=-1.440 | V_limit=-1.286 | B=653243.240
    -> heliacal rise detected
       criterion satisfied: h_star > 0° and V <= V_limit
  day 016 | 139-07-16 00:00:00 | visible=True | h_star=2.237° | h_sun=-8.159° | V=-1.440 | V_limit=-1.098 | B=200103.753
  day 017 | 139-07-17 00:00:00 | visible=True | h_star=1.938° | h_sun=-9.314° | V=-1.440 | V_limit=-1.136 | B=64385.214
  day 018 | 139-07-18 00:00:00 | visible=True | h_star=1.640° | h_sun=-10.463° | V=-1.440 | V_limit=-1.354 | B=20733.486
  day 019 | 139-07-19 00:00:00 | visible=True | h_star=1.773° | h_sun=-11.221° | V=-1.440 | V_limit=-1.052 | B=9641.547
  day 020 | 139-07-20 00:00:00 | visible=True | h_star=1.480° | h_sun=-12.363° | V=-1.440 | V_limit=-1.426 | B=3090.394
  day 021 | 139-07-21 00:00:00 | visible=True | h_star=1.617° | h_sun=-13.117° | V=-1.440 | V_limit=-1.033 | B=1425.541
  day 022 | 139-07-22 00:00:00 | visible=True | h_star=0.899° | h_sun=-14.632° | V=-1.440 | V_limit=-1.192 | B=288.621
  day 023 | 139-07-23 00:00:00 | visible=True | h_star=0.610° | h_sun=-15.759° | V=-1.440 | V_limit=-1.085 | B=67.119
  day 024 | 139-07-24 00:00:00 | visible=True | h_star=0.323° | h_sun=-16.882° | V=-1.440 | V_limit=-0.648 | B=0.000
  day 025 | 139-07-25 00:00:00 | visible=True | h_star=0.467° | h_sun=-17.628° | V=-1.440 | V_limit=-0.129 | B=0.000
  day 026 | 139-07-26 00:00:00 | visible=True | h_star=1.042° | h_sun=-18.000° | V=-1.440 | V_limit=1.465 | B=0.000
  day 027 | 139-07-27 00:00:00 | visible=True | h_star=2.045° | h_sun=-18.000° | V=-1.440 | V_limit=3.177 | B=0.000
  day 028 | 139-07-28 00:00:00 | visible=True | h_star=3.046° | h_sun=-18.000° | V=-1.440 | V_limit=4.206 | B=0.000
  day 029 | 139-07-29 00:00:00 | visible=True | h_star=4.044° | h_sun=-18.000° | V=-1.440 | V_limit=4.886 | B=0.000
  day 030 | 139-07-30 00:00:00 | visible=True | h_star=5.040° | h_sun=-18.000° | V=-1.440 | V_limit=5.366 | B=0.000
Heliacal rise of Sirius — Schaefer 1985 — 1 date(s)
  [1] 139-07-15 00:00:00  04:41:53.023  0139-07-14 00:00:00.000000  [hrw 2921] I mesut 1  Sirius 2.54°  Sun -6.95°
  source: Schaefer, B.E. 1985, Sky & Telescope 70, 261–263 (BASIC listing, lines 34–35 and 55–81).

Save the result:

[13]:

rise_schaefer1985_jed = rises.day_jed[0]

3.4 Belokrylov et al. (2011)

This model also scans the twilight. It uses reference_extinction as the reference extinction coefficient to correct the star’s magnitude.

[14]:

model = montu.HeliacalRise(
    model="belokrylov2011",
    k=0.25,  # atmospheric extinction coefficient
    reference_extinction=0.25,  # reference extinction coefficient
    step_minutes=2.0,  # time step
)

rises = model.compute(sirius, site, start_date, end_date, verbose=True)
model.print_rises(rises, title="Heliacal rise of Sirius — Belokrylov 2011", body_label="Sirius")

HeliacalRise verbose — model=belokrylov2011
  quantities: h_star: object altitude; h_sun: solar altitude; m': extinction-corrected magnitude; rho: Sun-object separation; h_theor: limiting solar altitude
  interval: 139-07-01 00:00:00 -> 139-07-30 00:00:00
  model parameters:
    k=0.25
    reference_extinction=0.25
    step_minutes=2.0
    twilight_sunbelow=-18.0
  criterion: h_star > 0° and h_sun <= h_theor
  morning scan: h_sun=-18.0° to sunrise, every 2.0 minutes
  day 001 | 139-07-01 00:00:00 | visible=False
  day 002 | 139-07-02 00:00:00 | visible=False
  day 003 | 139-07-03 00:00:00 | visible=False
  day 004 | 139-07-04 00:00:00 | visible=False
  day 005 | 139-07-05 00:00:00 | visible=False
  day 006 | 139-07-06 00:00:00 | visible=False
  day 007 | 139-07-07 00:00:00 | visible=False
  day 008 | 139-07-08 00:00:00 | visible=False
  day 009 | 139-07-09 00:00:00 | visible=False
  day 010 | 139-07-10 00:00:00 | visible=False
  day 011 | 139-07-11 00:00:00 | visible=False
  day 012 | 139-07-12 00:00:00 | visible=False
  day 013 | 139-07-13 00:00:00 | visible=False
  day 014 | 139-07-14 00:00:00 | visible=True | h_star=2.417° | h_sun=-6.105° | m'=2.628 | rho=48.064° | h_theor=-6.039°
    -> heliacal rise detected
       criterion satisfied: h_star > 0° and h_sun <= h_theor
  day 015 | 139-07-15 00:00:00 | visible=True | h_star=1.684° | h_sun=-7.776° | m'=3.575 | rho=48.591° | h_theor=-7.186°
  day 016 | 139-07-16 00:00:00 | visible=True | h_star=1.381° | h_sun=-8.944° | m'=4.087 | rho=49.113° | h_theor=-7.797°
  day 017 | 139-07-17 00:00:00 | visible=True | h_star=1.509° | h_sun=-9.704° | m'=3.859 | rho=49.644° | h_theor=-7.499°
  day 018 | 139-07-18 00:00:00 | visible=True | h_star=1.211° | h_sun=-10.851° | m'=4.414 | rho=50.184° | h_theor=-9.409°
  day 019 | 139-07-19 00:00:00 | visible=True | h_star=1.345° | h_sun=-11.608° | m'=4.154 | rho=50.734° | h_theor=-7.825°
  day 020 | 139-07-20 00:00:00 | visible=True | h_star=1.051° | h_sun=-12.747° | m'=4.757 | rho=51.291° | h_theor=-11.296°
  day 021 | 139-07-21 00:00:00 | visible=True | h_star=1.188° | h_sun=-13.500° | m'=4.461 | rho=51.859° | h_theor=-9.616°
  day 022 | 139-07-22 00:00:00 | visible=True | h_star=0.899° | h_sun=-14.632° | m'=5.118 | rho=52.433° | h_theor=-13.281°
  day 023 | 139-07-23 00:00:00 | visible=True | h_star=1.040° | h_sun=-15.382° | m'=4.783 | rho=53.017° | h_theor=-11.384°
  day 024 | 139-07-24 00:00:00 | visible=True | h_star=0.753° | h_sun=-16.507° | m'=5.500 | rho=53.607° | h_theor=-15.381°
  day 025 | 139-07-25 00:00:00 | visible=True | h_star=0.897° | h_sun=-17.254° | m'=5.123 | rho=54.206° | h_theor=-13.247°
  day 026 | 139-07-26 00:00:00 | visible=True | h_star=1.042° | h_sun=-18.000° | m'=4.778 | rho=54.812° | h_theor=-11.295°
  day 027 | 139-07-27 00:00:00 | visible=True | h_star=2.045° | h_sun=-18.000° | m'=3.066 | rho=55.427° | h_theor=-6.328°
  day 028 | 139-07-28 00:00:00 | visible=True | h_star=3.046° | h_sun=-18.000° | m'=2.038 | rho=56.049° | h_theor=-5.042°
  day 029 | 139-07-29 00:00:00 | visible=True | h_star=4.044° | h_sun=-18.000° | m'=1.358 | rho=56.678° | h_theor=-4.184°
  day 030 | 139-07-30 00:00:00 | visible=True | h_star=5.040° | h_sun=-18.000° | m'=0.878 | rho=57.312° | h_theor=-3.573°
Heliacal rise of Sirius — Belokrylov 2011 — 1 date(s)
  [1] 139-07-14 00:00:00  04:45:14.533  0139-07-13 00:00:00.000000  [hrw 2921] IV shemu 30  Sirius 2.42°  Sun -6.11°
  source: Belokrylov, R. O., Belokrylov, S. V., & Nickiforov, M. G. (2011). "Model of the stellar visibility during twilight". Bulgarian Astronomical Journal, 16, 50-72.

Save the result:

[15]:

rise_belokrylov2011_jed = rises.day_jed[0]

4. Opening the black box

Next we explain, in basic terms, what the HeliacalRise.compute method does internally for each model, so that it is clear how MontuPython’s own tools are used for this important calculation, and also which formulas each model applies.

To do this we use two dates in each case: one on which the heliacal rise does not occur and one on which it does, so we can see when the criterion fails and when it is satisfied.

4.1. Ptolemy model: Arcus Visionis

The source used by MontuPython for this model is: Toomer, G. J. (1998). Ptolemy’s Almagest. Princeton University Press. Book XIII, Chapter 7: “On the heliacal risings and settings of the planets”

This model asks where the Sun is when Sirius rises on the eastern horizon (\(h_\star=0\)). If \(h_\odot\) is the solar altitude at Sirius’s rising time, the calculated arc of vision is

\[AV_{\rm calc}\equiv h_\star-h_\odot=-h_\odot\]

Ptolemy’s method assumes that Sirius will be visible if \(AV_{\rm calc}\ge AV_{\rm crit}\).

MontuPython adopts \(AV_{\rm crit}=15°\) for first-magnitude stars. This is a historical threshold, not a consequence of geometry. For other bodies, the value of \(AV_{\rm crit}\) changes according to Ptolemy’s accumulated experience and his sources:

  • Venus: \(5^\circ\)

  • Jupiter: \(10^\circ\)

  • Mercury: \(10^\circ\)

  • Saturn: \(11^\circ\)

  • Mars: \(11.5^\circ\)

To calculate the Arcus Visionis we need, each day: 1) Sirius’s rising time, 2) the Sun’s altitude above the horizon, and 3) the Arcus Visionis between the Sun and Sirius at rising. In MontuPython we do it like this:

Let’s try first with the day before the heliacal rise. The rising time is:

[16]:

t_test = montu.Time(rise_ptolemy_jed - 1, format="jd", calendar="mixed")
conditions_sirius = sirius.conditions_in_sky(t_test, site)
jed_rise_star_utc = conditions_sirius.rise_time[0]
site.get_local_time(jed_rise_star_utc)

[16]:

'04:11:49.837'

Sirius’s and the Sun’s altitudes at rising are:

[17]:

mt_rise = montu.Time(jed_rise_star_utc, format="jd", calendar="mixed")
h_sirius = sirius.where_in_sky(mt_rise, site).el[0]
sun.where_in_sky(mt_rise, site)
h_sun = sun.position.el

arcus_visionis_calc = h_sirius - h_sun

h_sirius, h_sun, arcus_visionis_calc

[17]:

(np.float64(-0.5349187387273291),
 -13.28456084132389,
 np.float64(12.74964210259656))

As we can see, the calculated Arcus Visionis for that day is smaller than the threshold \(AV_{\rm crit}=15°\). Therefore that day is not considered the heliacal-rise day.

Now let’s try the heliacal-rise day:

[18]:

t_test = montu.Time(rise_ptolemy_jed, format="jd", calendar="mixed")
conditions_sirius = sirius.conditions_in_sky(t_test, site)
jed_rise_star_utc = conditions_sirius.rise_time[0]

mt_rise = montu.Time(jed_rise_star_utc, format="jd", calendar="mixed")
h_sirius = sirius.where_in_sky(mt_rise, site).el[0]
sun.where_in_sky(mt_rise, site)
h_sun = sun.position.el

arcus_visionis_calc = h_sirius - h_sun

site.get_local_time(jed_rise_star_utc), h_sirius, h_sun, arcus_visionis_calc

[18]:

('04:07:53.948',
 np.float64(-0.5348629198277512),
 -14.153722050001877,
 np.float64(13.618859130174126))

On this day the Arcus Visionis is already larger, and therefore this is the heliacal rise.

The following day gives:

[19]:

t_test = montu.Time(rise_ptolemy_jed + 1, format="jd", calendar="mixed")
conditions_sirius = sirius.conditions_in_sky(t_test, site)
jed_rise_star_utc = conditions_sirius.rise_time[0]

mt_rise = montu.Time(jed_rise_star_utc, format="jd", calendar="mixed")
h_sirius = sirius.where_in_sky(mt_rise, site).el[0]
sun.where_in_sky(mt_rise, site)
h_sun = sun.position.el

arcus_visionis_calc = h_sirius - h_sun

site.get_local_time(jed_rise_star_utc), h_sirius, h_sun, arcus_visionis_calc

[19]:

('04:03:58.050',
 np.float64(-0.5348388207339038),
 -15.02248540016071,
 np.float64(14.487646579426805))

The Arcus Visionis is already greater than the critical value.

It is worth noting that Sirius’s altitude at rising is negative because MontuPython’s rise and set models include atmospheric refraction, which in a standard atmosphere is close to 34 arcminutes. MontuPython actually computes this refraction from the environmental conditions of the observing site (see the Observer class).

4.2. Schaefer (1987): fixed solar depression

The source used by MontuPython is Schaefer, B. E. (1987), “Heliacal rise phenomena”, Journal for the History of Astronomy 18(11), 19–33.

This model does not wait for Sirius’s rising time. It evaluates the sky once each morning: when the Sun reaches the chosen depression \(h_\odot=-11°\). MontuPython obtains that instant with Sun.when_is_twilight().

For Sirius at altitude \(h_\star>0\), the air mass is approximated by the approximation of a plain-parallel atmosphere:

\[X=\csc h_\star=\frac{1}{\sin h_\star}.\]

With catalogue magnitude \(V\), extinction coefficient \(k\), and zenith limiting magnitude \(V_{\rm lim,z}\), the observed magnitude of the star \(V_{\rm obs}\) and the limiting magnitude at the height at which it is observed \(V_{\rm lim}(h_\star)\) is computed using:

\[V_{\rm obs}=V+kX,\]
\[V_{\rm lim}(h_\star)=V_{\rm lim,z}-k(X-1).\]

Since smaller magnitudes correspond to brighter objects, a margin \(V_{\rm lim}-V_{\rm obs}<0\) means that Sirius is still not visible. The criterion for observability is:

\[h_\star>0°\quad\text{and}\quad V_{\rm obs}\le V_{\rm lim}(h_\star).\]

This criterion is fullfilled at many nigths during the year. For finding the heliacal rise, the algorithm sweep several twilights, finding the first where the condition is fulfilled.

Each morning we need: 1) the instant when the Sun is at \(-11°\), 2) Sirius’s altitude at that instant, and 3) the air mass and the two magnitudes. We start with the day before the heliacal rise:

[20]:

t_test = montu.Time(rise_schaefer1987_jed - 1, format="jd", calendar="mixed")

# Morning instant with the Sun at -11°
jed_twilight = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-11.0))
mt_twilight = montu.Time(jed_twilight, format="jd", calendar="mixed")

site.get_local_time(jed_twilight)

[20]:

'04:24:19.401'

Now we get the elevation of the star:

[21]:

h_sirius = sirius.where_in_sky(mt_twilight, site).el[0]

k = 0.25
V_lim_zenith = 6.0
X = 1 / np.sin(np.deg2rad(h_sirius))
V_obs = V_SIRIUS + k * X
V_lim = V_lim_zenith - k * (X - 1)

print(f"Tested date: {t_test.readable.datemix}")
print(f"Elevation of Sirius at twilight ({site.get_local_time(jed_twilight)}): {h_sirius:.2f}°")
print(f"Air mass: {X:.2f}")
print(f"Observed magnitude: {V_obs:.2f}")
print(f"Limiting magnitude: {V_lim:.2f}")

Tested date: 139-07-20 00:00:00
Elevation of Sirius at twilight (04:24:19.401): 2.99°
Air mass: 19.19
Observed magnitude: 3.36
Limiting magnitude: 1.45

The previous day Sirius is already above the horizon when the Sun reaches \(-11°\), but \(V_{\rm obs}>V_{\rm lim}\): its extinguished magnitude is still too faint compared with the local threshold. That is why the morning does not count as a heliacal rise.

The heliacal-rise day:

[22]:

t_test = montu.Time(rise_schaefer1987_jed, format="jd", calendar="mixed")

jed_twilight = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-11.0))
mt_twilight = montu.Time(jed_twilight, format="jd", calendar="mixed")

h_sirius = sirius.where_in_sky(mt_twilight, site).el[0]
sun.where_in_sky(mt_twilight, site)
h_sun = sun.position.el

X = 1 / np.sin(np.deg2rad(h_sirius))
V_obs = V_SIRIUS + k * X
V_lim = V_lim_zenith - k * (X - 1)

print(f"Tested date: {t_test.readable.datemix} ({t_test.readable.datesot})")
print(f"Elevation of Sirius at twilight ({site.get_local_time(jed_twilight)}): {h_sirius:.2f}°")
print(f"Air mass: {X:.2f}")
print(f"Observed magnitude: {V_obs:.2f}")
print(f"Limiting magnitude: {V_lim:.2f}")

Tested date: 139-07-21 00:00:00 ([hrw 2922] I akhet 2)
Elevation of Sirius at twilight (04:24:56.311): 3.96°
Air mass: 14.50
Observed magnitude: 2.18
Limiting magnitude: 2.63

Now Sirius is higher (\(h_\star\approx 4.5°\)), the air mass is lower, and the margin becomes positive. That is the first visible morning according to Schaefer (1987).

4.3. Schaefer (1985): twilight scan

The source used by MontuPython is Schaefer, B. E. (1985), “Predicting Heliacal Risings and Settings”, Sky & Telescope 70, $$ dawn (\(h_\odot=-18°\)) to sunrise in steps of 2 minutes. At each step it recalculates the positions of the Sun and Sirius. If Sirius is below the horizon, it moves to the next step.

The air mass uses Rozenberg’s formula, which is more stable near the horizon:

\[X=\frac{1}{\cos z+0.025\exp(-11\cos z)},\qquad z=90°-h_\star.\]

At each step, with \(z\), \(h_\odot\), and the azimuthal separation \(\Delta A\) in radians it compute the following parameter:

\[x=-0.2(V_{\rm lim,z}-7.93+k),\]
\[B_0=79.4(10^x-1)^2-589k,\]
\[L_5=4.75-\frac{\Delta A\,z}{3}+h_\odot(8.2z+12)+2.86z.\]

This parameters allows us to calculate the sky-brightnesss:

\[\begin{split}B = \begin{cases} B_0 + \dfrac{k}{0.20} \cdot 10^{L_5} & \text{si } L_5 \geq -2.07 \\ B_0 + 589k & \text{si } L_5 < -2.07 \end{cases}\end{split}\]

The physiological threshold is

\[E_{\rm th}=C_5(1+\sqrt{K_5B})^2,\]

with \((C_5,K_5)=(1.58\times10^{-10},0.0126)\) if \(B<1649\), and \((4.4668\times10^{-9},1.258\times10^{-6})\) otherwise.

Finally the limiting magnitude is given by,

\[V_{\rm lim}=-16.57-kX-2.5\log_{10}E_{\rm th}.\]

Sirius is visible when \(V\le V_{\rm lim}\). The first step that satisfies the criterion is the reported instant.

Let’s first define Rozenberg’s air mass and the model parameters. Then we evaluate the day before the heliacal rise:

[23]:

def rozenberg_airmass(h_star):
    z = np.deg2rad(90.0 - h_star)
    cos_z = np.cos(z)
    return 1 / (cos_z + 0.025 * np.exp(-11 * cos_z))

k = 0.25
V_lim_zenith = 6.0
step_minutes = 2.0

t_test = montu.Time(rise_schaefer1985_jed - 1, format="jd", calendar="mixed")
dawn_jed = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-18.0))
sun.conditions_in_sky(t_test, site)
sunrise_jed = sun.condition.rise_time

print("Scan window:", site.get_local_time(dawn_jed), "→", site.get_local_time(sunrise_jed))

Scan window: 03:43:14.538 → 05:12:00.485

Let’s scan:

[24]:

found = False
for at_jed in np.arange(dawn_jed, sunrise_jed, step_minutes / 1440):
    mt = montu.Time(at_jed, format="jd", calendar="mixed")
    star = sirius.where_in_sky(mt, site).iloc[0]
    if star.el <= 0:
        continue

    # Compute the air mass
    z = np.deg2rad(90.0 - star.el)
    h_sun = np.deg2rad(sun.position.el)
    X = rozenberg_airmass(star.el)

    # Compute the angular separation
    sun.where_in_sky(mt, site)
    delta_az = np.deg2rad(abs((star.az - sun.position.az + 180) % 360 - 180))

    # Compute the parameters for the sky brightness
    x = -0.2 * (V_lim_zenith - 7.93 + k)
    B0 = 79.4 * (10**x - 1)**2 - 589 * k
    L5 = 4.75 - delta_az * z / 3 + h_sun * (8.2 * z + 12) + 2.86 * z
    B = B0 + (k / 0.20) * 10**L5 if L5 >= -2.07 else B0 + 589 * k
    C5, K5 = (1.58e-10, 0.0126) if B < 1649 else (4.4668e-9, 1.258e-6)
    E_th = C5 * (1 + np.sqrt(K5 * B))**2
    V_lim = -16.57 - k * X - 2.5 * np.log10(E_th)

    # Check if the star is visible
    visible = V_SIRIUS <= V_lim

    print(
        site.get_local_time(at_jed),
        f"| h_star={star.el:5.2f}°",
        f"| h_sun={sun.position.el:6.2f}°",
        f"| V_lim={V_lim:6.2f}",
        f"| V_star={V_SIRIUS:6.2f}",
        f"| margin={V_lim - V_SIRIUS:6.2f}",
        f"| visible={visible}",
    )
    if visible:
        found = True
        break

print()
print("Was any visible instant found?", found)

04:35:14.537 | h_star= 0.27° | h_sun= -8.20° | V_lim= -2.52 | V_star= -1.44 | margin= -1.08 | visible=False
04:37:14.538 | h_star= 0.70° | h_sun= -7.80° | V_lim= -3.93 | V_star= -1.44 | margin= -2.49 | visible=False
04:39:14.539 | h_star= 1.13° | h_sun= -7.39° | V_lim= -3.02 | V_star= -1.44 | margin= -1.58 | visible=False
04:41:14.540 | h_star= 1.56° | h_sun= -6.97° | V_lim= -2.38 | V_star= -1.44 | margin= -0.94 | visible=False
04:43:14.532 | h_star= 1.99° | h_sun= -6.55° | V_lim= -1.93 | V_star= -1.44 | margin= -0.49 | visible=False
04:45:14.533 | h_star= 2.42° | h_sun= -6.11° | V_lim= -1.64 | V_star= -1.44 | margin= -0.20 | visible=False
04:47:14.534 | h_star= 2.84° | h_sun= -5.64° | V_lim= -1.46 | V_star= -1.44 | margin= -0.02 | visible=False
04:49:14.535 | h_star= 3.27° | h_sun= -5.12° | V_lim= -1.40 | V_star= -1.44 | margin=  0.04 | visible=True

Was any visible instant found? True

The previous day Sirius does rise during twilight, but in every step \(V>V_{\rm lim}\): the sky is still too bright or the star too low. There is no heliacal rise.

Let’s repeat the same scan on the heliacal-rise day:

[25]:

t_test = montu.Time(rise_schaefer1985_jed, format="jd", calendar="mixed")
dawn_jed = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-18.0))
sun.conditions_in_sky(t_test, site)
sunrise_jed = sun.condition.rise_time

print("Scan window:", site.get_local_time(dawn_jed), "→", site.get_local_time(sunrise_jed))
print()

for at_jed in np.arange(dawn_jed, sunrise_jed, step_minutes / 1440):
    mt = montu.Time(at_jed, format="jd", calendar="mixed")
    star = sirius.where_in_sky(mt, site).iloc[0]
    if star.el <= 0:
        continue

    sun.where_in_sky(mt, site)
    z = np.deg2rad(90.0 - star.el)
    h_sun = np.deg2rad(sun.position.el)
    delta_az = np.deg2rad(abs((star.az - sun.position.az + 180) % 360 - 180))
    X = rozenberg_airmass(star.el)

    x = -0.2 * (V_lim_zenith - 7.93 + k)
    B0 = 79.4 * (10**x - 1)**2 - 589 * k
    L5 = 4.75 - delta_az * z / 3 + h_sun * (8.2 * z + 12) + 2.86 * z
    B = B0 + (k / 0.20) * 10**L5 if L5 >= -2.07 else B0 + 589 * k
    C5, K5 = (1.58e-10, 0.0126) if B < 1649 else (4.4668e-9, 1.258e-6)
    E_th = C5 * (1 + np.sqrt(K5 * B))**2
    V_lim = -16.57 - k * X - 2.5 * np.log10(E_th)
    visible = V_SIRIUS <= V_lim

    print(
        site.get_local_time(at_jed),
        f"| h_star={star.el:5.2f}°",
        f"| h_sun={sun.position.el:6.2f}°",
        f"| V_lim={V_lim:6.2f}",
        f"| V_star={V_SIRIUS:6.2f}",
        f"| margin={V_lim - V_SIRIUS:6.2f}",
        f"| visible={visible}",
    )
    if visible:
        print()
        print("First visible instant → this is the heliacal rise according to Schaefer (1985).")
        break

Scan window: 03:43:53.029 → 05:12:29.558

04:31:53.027 | h_star= 0.40° | h_sun= -8.96° | V_lim= -4.66 | V_star= -1.44 | margin= -3.22 | visible=False
04:33:53.028 | h_star= 0.83° | h_sun= -8.57° | V_lim= -3.46 | V_star= -1.44 | margin= -2.02 | visible=False
04:35:53.028 | h_star= 1.26° | h_sun= -8.18° | V_lim= -2.61 | V_star= -1.44 | margin= -1.17 | visible=False
04:37:53.029 | h_star= 1.68° | h_sun= -7.78° | V_lim= -2.01 | V_star= -1.44 | margin= -0.57 | visible=False
04:39:53.030 | h_star= 2.11° | h_sun= -7.37° | V_lim= -1.58 | V_star= -1.44 | margin= -0.14 | visible=False
04:41:53.023 | h_star= 2.54° | h_sun= -6.95° | V_lim= -1.29 | V_star= -1.44 | margin=  0.15 | visible=True

First visible instant → this is the heliacal rise according to Schaefer (1985).

The transition is clear: for several minutes the margin is negative, and at the first step with a positive margin the model declares the heliacal rise.

4.4. Belokrylov et al. (2011): limiting solar altitude

The source used by MontuPython is Belokrylov, R. O., Belokrylov, S. V., and Nickiforov, M. G. (2011), “Model of the stellar visibility during twilight”, Bulgarian Astronomical Journal 16, 50–72, equations (5)–(8).

This model also scans the twilight from \(h_\odot=-18°\) to sunrise, but translates visibility into a limiting solar altitude. First it corrects the magnitude for extinction (with \(k_0=0.25\)):

\[m'=V+k(X-1)+(k-k_0).\]

Then it obtains the limiting solar altitude (this comes from empirical fits):

\[\begin{split}h_{\rm lim}=\begin{cases} -2.47-1.23m', & m'<4.2,\\ 15.62-5.61m', & m'\ge4.2. \end{cases}\end{split}\]

The Sun–Sirius angular separation is

\[\cos\rho=\sin h_\star\sin h_\odot+\cos h_\star\cos h_\odot\cos(A_\star-A_\odot),\]

and the correction for proximity to the Sun,

\[\Delta h=-0.0338\max(0,58°-\rho),\qquad h_{\rm theor}=h_{\rm lim}+\Delta h.\]

The star is visible if \(h_\odot\le h_{\rm theor}\).

Let’s evaluate the day before the heliacal rise with the same kind of scan:

[26]:

k = 0.25
k0 = 0.25
step_minutes = 2.0

t_test = montu.Time(rise_belokrylov2011_jed - 1, format="jd", calendar="mixed")
dawn_jed = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-18.0))
sun.conditions_in_sky(t_test, site)
sunrise_jed = sun.condition.rise_time

print("Scan window:", site.get_local_time(dawn_jed), "→", site.get_local_time(sunrise_jed))

Scan window: 03:42:36.885 → 05:11:31.895

Let’s scan:

[27]:

found = False
for at_jed in np.arange(dawn_jed, sunrise_jed, step_minutes / 1440):

    mt = montu.Time(at_jed, format="jd", calendar="mixed")
    star = sirius.where_in_sky(mt, site).iloc[0]
    if star.el <= 0:
        continue

    # Compute the air mass
    X = rozenberg_airmass(star.el)
    m_prime = V_SIRIUS + k * (X - 1) + (k - k0)
    h_lim = -2.47 - 1.23 * m_prime if m_prime < 4.2 else 15.62 - 5.61 * m_prime

    # Compute the angular separation
    sun.where_in_sky(mt, site)
    h1, A1 = np.deg2rad(star.el), np.deg2rad(star.az)
    h2, A2 = np.deg2rad(sun.position.el), np.deg2rad(sun.position.az)
    cos_rho = np.sin(h1) * np.sin(h2) + np.cos(h1) * np.cos(h2) * np.cos(A1 - A2)
    rho = np.rad2deg(np.arccos(np.clip(cos_rho, -1, 1)))

    # Compute the correction for proximity to the Sun
    delta_h = -0.0338 * max(0, 58 - rho)
    h_theor = h_lim + delta_h
    visible = sun.position.el <= h_theor

    print(
        site.get_local_time(at_jed),
        f"| h_star={star.el:5.2f}°",
        f"| h_sun={sun.position.el:6.2f}°",
        f"| m'={m_prime:5.2f}",
        f"| h_theor={h_theor:6.2f}°",
        f"| visible={visible}",
    )
    if visible:
        found = True
        break

print()
print("Was any visible instant found?", found)

04:38:36.886 | h_star= 0.15° | h_sun= -7.41° | m'= 7.59 | h_theor=-27.31° | visible=False
04:40:36.887 | h_star= 0.58° | h_sun= -7.00° | m'= 6.00 | h_theor=-18.37° | visible=False
04:42:36.888 | h_star= 1.01° | h_sun= -6.57° | m'= 4.85 | h_theor=-11.92° | visible=False
04:44:36.889 | h_star= 1.44° | h_sun= -6.13° | m'= 3.98 | h_theor= -7.72° | visible=False
04:46:36.890 | h_star= 1.87° | h_sun= -5.67° | m'= 3.30 | h_theor= -6.88° | visible=False
04:48:36.882 | h_star= 2.30° | h_sun= -5.15° | m'= 2.76 | h_theor= -6.22° | visible=False
04:50:36.883 | h_star= 2.72° | h_sun= -4.53° | m'= 2.32 | h_theor= -5.68° | visible=False
04:52:36.884 | h_star= 3.15° | h_sun= -3.56° | m'= 1.95 | h_theor= -5.23° | visible=False
04:54:36.885 | h_star= 3.58° | h_sun= -2.86° | m'= 1.64 | h_theor= -4.85° | visible=False
04:56:36.886 | h_star= 4.00° | h_sun= -2.47° | m'= 1.38 | h_theor= -4.53° | visible=False
04:58:36.887 | h_star= 4.43° | h_sun= -2.16° | m'= 1.15 | h_theor= -4.25° | visible=False
05:00:36.888 | h_star= 4.85° | h_sun= -1.88° | m'= 0.96 | h_theor= -4.01° | visible=False
05:02:36.889 | h_star= 5.28° | h_sun= -1.60° | m'= 0.78 | h_theor= -3.79° | visible=False
05:04:36.890 | h_star= 5.70° | h_sun= -1.32° | m'= 0.63 | h_theor= -3.60° | visible=False
05:06:36.882 | h_star= 6.12° | h_sun= -1.03° | m'= 0.50 | h_theor= -3.44° | visible=False
05:08:36.883 | h_star= 6.54° | h_sun= -0.73° | m'= 0.37 | h_theor= -3.29° | visible=False
05:10:36.884 | h_star= 6.96° | h_sun= -0.41° | m'= 0.27 | h_theor= -3.15° | visible=False

Was any visible instant found? False

The previous day the Sun never becomes low enough relative to \(h_{\rm theor}\). There is no heliacal rise.

The heliacal-rise day:

[28]:

t_test = montu.Time(rise_belokrylov2011_jed, format="jd", calendar="mixed")
dawn_jed = min(montu.Sun.when_is_twilight(t_test, site, sunbelow=-18.0))
sun.conditions_in_sky(t_test, site)
sunrise_jed = sun.condition.rise_time

print("Scan window:", site.get_local_time(dawn_jed), "→", site.get_local_time(sunrise_jed))
print()

for at_jed in np.arange(dawn_jed, sunrise_jed, step_minutes / 1440):
    mt = montu.Time(at_jed, format="jd", calendar="mixed")
    star = sirius.where_in_sky(mt, site).iloc[0]
    if star.el <= 0:
        continue

    sun.where_in_sky(mt, site)
    X = rozenberg_airmass(star.el)
    m_prime = V_SIRIUS + k * (X - 1) + (k - k0)
    h_lim = -2.47 - 1.23 * m_prime if m_prime < 4.2 else 15.62 - 5.61 * m_prime

    h1, A1 = np.deg2rad(star.el), np.deg2rad(star.az)
    h2, A2 = np.deg2rad(sun.position.el), np.deg2rad(sun.position.az)
    cos_rho = np.sin(h1) * np.sin(h2) + np.cos(h1) * np.cos(h2) * np.cos(A1 - A2)
    rho = np.rad2deg(np.arccos(np.clip(cos_rho, -1, 1)))
    delta_h = -0.0338 * max(0, 58 - rho)
    h_theor = h_lim + delta_h
    visible = sun.position.el <= h_theor

    print(
        site.get_local_time(at_jed),
        f"| h_star={star.el:5.2f}°",
        f"| h_sun={sun.position.el:6.2f}°",
        f"| m'={m_prime:5.2f}",
        f"| rho={rho:5.1f}°",
        f"| h_theor={h_theor:6.2f}°",
        f"| visible={visible}",
    )
    if visible:
        print()
        print("First visible instant → this is the heliacal rise according to Belokrylov (2011).")
        break

Scan window: 03:43:14.538 → 05:12:00.485

04:35:14.537 | h_star= 0.27° | h_sun= -8.20° | m'= 7.09 | rho= 48.1° | h_theor=-24.47° | visible=False
04:37:14.538 | h_star= 0.70° | h_sun= -7.80° | m'= 5.64 | rho= 48.1° | h_theor=-16.35° | visible=False
04:39:14.539 | h_star= 1.13° | h_sun= -7.39° | m'= 4.58 | rho= 48.1° | h_theor=-10.41° | visible=False
04:41:14.540 | h_star= 1.56° | h_sun= -6.97° | m'= 3.77 | rho= 48.1° | h_theor= -7.45° | visible=False
04:43:14.532 | h_star= 1.99° | h_sun= -6.55° | m'= 3.14 | rho= 48.1° | h_theor= -6.67° | visible=False
04:45:14.533 | h_star= 2.42° | h_sun= -6.11° | m'= 2.63 | rho= 48.1° | h_theor= -6.04° | visible=True

First visible instant → this is the heliacal rise according to Belokrylov (2011).

Here the criterion is read differently: instead of comparing magnitudes, we check whether the Sun’s actual altitude is already below the theoretical altitude \(h_{\rm theor}\) allowed for that position of Sirius.

5. Comparison of the four models

Let’s summarize the dates obtained with the HeliacalRise interface in section 3. All four models use the same site, the same interval, and the same star; only the visibility criterion changes.

[29]:

print("Ptolemy")
print("  Date:", montu.Time(rise_ptolemy_jed, format="jd").readable.datemix)
print("  Sothic:", montu.Time(rise_ptolemy_jed, format="jd").readable.datesot)
print()

print("Schaefer 1987")
print("  Date:", montu.Time(rise_schaefer1987_jed, format="jd").readable.datemix)
print("  Sothic:", montu.Time(rise_schaefer1987_jed, format="jd").readable.datesot)
print()

print("Schaefer 1985")
print("  Date:", montu.Time(rise_schaefer1985_jed, format="jd").readable.datemix)
print("  Sothic:", montu.Time(rise_schaefer1985_jed, format="jd").readable.datesot)
print()

print("Belokrylov 2011")
print("  Date:", montu.Time(rise_belokrylov2011_jed, format="jd").readable.datemix)
print("  Sothic:", montu.Time(rise_belokrylov2011_jed, format="jd").readable.datesot)

Ptolemy
  Date: 139-07-20 00:00:00
  Sothic: [hrw 2922] I akhet 1

Schaefer 1987
  Date: 139-07-21 00:00:00
  Sothic: [hrw 2922] I akhet 2

Schaefer 1985
  Date: 139-07-15 00:00:00
  Sothic: [hrw 2921] I mesut 1

Belokrylov 2011
  Date: 139-07-14 00:00:00
  Sothic: [hrw 2921] IV shemu 30

With the adopted parameters, Ptolemy and Schaefer (1987) place the transition on 19 July 2782 BCE; Schaefer (1985) and Belokrylov et al. (2011) place it on 11 July. The eight-day difference is not a numerical error: each algorithm defines «visible» differently.

In addition, \(k\) changes with aerosols and humidity; \(V_{\rm lim,z}\) depends on the observer and the night-sky background; an elevated horizon delays the appearance; and clouds or haze can dominate the result. That is why a calculated heliacal-rise date should always be reported together with the model and its parameters.

6. Conclusions

We have deliberately separated two layers:

  1. geometry, computed by MontuPython with Time, Observer, Stars.where_in_sky, Stars.conditions_in_sky, Sun.where_in_sky, Sun.when_is_twilight, and Sun.conditions_in_sky;

  2. visibility, reconstructed in open cells from the equations of each model.

The heliacal rise is the first morning that goes from invisible to visible, not simply the first visible instant found within an interval. The step-by-step reconstructions in section 4 show why one day fails the criterion and the next one satisfies it.

This disagreement between models is scientifically informative: it quantifies how «first visibility» depends on the observational model. For a historical study it is advisable always to report site, horizon, atmospheric parameters, model, time step, and calendar system, rather than presenting a calculated date as exact.

Implemented references

  • Schaefer, B. E. (1985), “Predicting Heliacal Risings and Settings”, Sky & Telescope 70, 261–263; BASIC listing, lines 34–35 and 55–81.

  • Schaefer, B. E. (1987), “Heliacal rise phenomena”, Journal for the History of Astronomy 18(11), 19–33.

  • Belokrylov, R. O., Belokrylov, S. V., and Nickiforov, M. G. (2011), “Model of the stellar visibility during twilight”, Bulgarian Astronomical Journal 16, 50–72, equations (5)–(8).

  • Toomer, G. J. (1998), Ptolemy’s Almagest, Princeton University Press, Book XIII, Chapter 7.

Powered by MontuPython. For more examples see MontuPython GitHub repo.

Jorge I. Zuluaga © 2023-present