Studying Ptolemy’s sun declination using visualization tools

from kanon.units import Sexagesimal
from kanon.tables import HTable
from astropy.units import arcsecond
from matplotlib import pyplot as plt
import math
import numpy as np

# This function redefines what is in the declination notebook

def build_declination(sin_table, obliquity):
    arcsin_table = sin_table.copy(set_index="Val")
    obl = sin_table.get(obliquity)
    obl_table = sin_table.apply("Val", lambda x: x * obl)
    return obl_table.apply("Val", lambda x: round(arcsin_table.get(x), 2))

OBLIQUITY = "23;51,20"

sin = lambda x: math.sin(math.radians(x))
asin = lambda x: math.degrees(math.asin(x))

# This table holds computed sine values
x = list(Sexagesimal.range(91))
y = [round(Sexagesimal.from_float(sin(n), 3)) for n in x]
sin_table_true = HTable([x, y], names=("Arg", "Val"), index="Arg")

# This table holds computed declination values
sin_obl = sin(Sexagesimal(OBLIQUITY))
x = list(Sexagesimal.range(1, 91))
y = [round(Sexagesimal.from_float(asin(sin(n) * sin_obl), 2)) for n in x]
decl_table_true = HTable(
    [x, y],
    names=("Arg", "Val"),
    index="Arg",
)

# We import Ptolemy's sun declination table from DISHAS
ptolemy_decl = HTable.read(214)

ptolemy_decl.plot2d()
plt.show()

# We rebuild a declination table from a sine table without odd arguments

sin_table_grid2 = sin_table_true[::2]
decl_table_grid2 = build_declination(sin_table_grid2, Sexagesimal(OBLIQUITY))[1:]

# Let's compare Ptolemy's declination with our reconstructed tables

residue_nogrid = ptolemy_decl.apply(
    "Entries", lambda x: x - decl_table_true["Val"], "Declination residue"
)
residue_nogrid["Declination residue"] = residue_nogrid["Declination residue"].to(
    arcsecond
)
# For the computed declination table residue, we plot in blue
residue_nogrid.plot2d("bx", linestyle="dashed", lw=0.4)

residue_grid2 = ptolemy_decl[1::2].apply(
    "Entries", lambda x: x - decl_table_grid2["Val"], "Declination residue"
)
residue_grid2["Declination residue"] = residue_grid2["Declination residue"].to(
    arcsecond
)
# For the built declination with grid 2 sine residue, we plot in red
residue_grid2.plot2d("r+", linestyle="dashed", lw=0.4)

plt.show()

arr_nogrid = residue_nogrid["Declination residue"].astype(float)
arr_grid2 = residue_grid2["Declination residue"].astype(float)

print("Residue with no interpolation")
print(f"mean : {np.mean(arr_nogrid):.4f}, std : {np.std(arr_nogrid):.4f}")
print(f"quartiles : {np.quantile(arr_nogrid, [0.25, 0.5, 0.75])} \n")
print("Residue with grid 2 interpolation")
print(f"mean : {np.mean(arr_grid2):.4f}, std : {np.std(arr_grid2):.4f}")
print(f"quartiles : {np.quantile(arr_grid2, [0.25, 0.5, 0.75])}")

Residue with no interpolation
mean : 0.3111, std : 2.2639
quartiles : Declination residue
       arcsec
-------------------
               -1.0
                0.0
                1.0

Residue with grid 2 interpolation
mean : -3.0222, std : 3.1516
quartiles : Declination residue
       arcsec
-------------------
               -5.0
               -3.0
               -1.0