mirror/astronomy
1
0 Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-08-07 05:22:16 -04:00
Code Issues Packages Projects Releases Wiki Activity

Improved agreement of precision among the 4 supported languages.

Browse Source 
Before making these changes, I had the following discrepancies
between the calculations made by the different programming
language implementations of Astronomy Engine:

    C vs C#: 5.55112e-17, worst line number = 6
    C vs JS: 2.78533e-12, worst line number = 196936
    C vs PY: 1.52767e-12, worst line number = 159834

Now the results are:

    Diffing calculations: C vs C#
    ctest(Diff): Maximum numeric difference = 5.55112e-17, worst line number = 5

    Diffing calculations: C vs JS
    ctest(Diff): Maximum numeric difference = 1.02318e-12, worst line number = 133677

    Diffing calculations: C vs PY
    ctest(Diff): Maximum numeric difference = 5.68434e-14, worst line number = 49066

    Diffing calculations: JS vs PY
    ctest(Diff): Maximum numeric difference = 1.02318e-12, worst line number = 133677

Here is how I did this:

1. Use new constants HOUR2RAD, RAD2HOUR that directly convert between radians and sidereal hours.
   This reduces tiny roundoff errors in the conversions.

2. In VSOP longitude calculations, keep clamping the angular sum to
   the range [-2pi, +2pi], to prevent it from accumulating thousands
   of radians. This reduces the accumulated error in the final result
   before it is fed into trig functions.

The remaining discrepancies are largely because of an "azimuth amplification" effect:
When converting equatorial coordinates to horizontal coordinates, an object near
the zenith (or nadir) has an azimuth that is highly sensitive to the input
equatorial coordinates. A tiny change in right ascension (RA) can cause a much
larger change in azimuth.

I tracked down the RA discrepancy, and it is due to a different behavior
of the atan2 function in C and JavaScript. There are cases where the least
significant decimal digit is off by 1, as if due to a difference of opinion
about rounding policy.

My best thought is to go back and have a more nuanced diffcalc that
applies less strict tests for azimuth values than the other calculated values.
It seems like every other computed quantity is less sensitive, because solar
system bodies tend to stay away from "poles" of other angular coordinate
systems: their ecliptic latitudes and equatorial declinations are usually
reasonably close to zero. Therefore, right ascensions and ecliptic longitudes
are usually insensitive to changes in the cartesian coordinates they
are calculated from.
This commit is contained in:
Don Cross committed 2021-04-18 21:15:17 -04:00
1 parent 6b01510b33
commit cbcacc4b57
32 files changed
+11067 -10503
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

generate/template/astronomy.c
+27 -15
View File
@@ -1342,7 +1342,7 @@ static astro_equatorial_t vector2radec(const double pos[3], astro_time_t time)
}
else
{
equ.ra = atan2(pos[1], pos[0]) / (DEG2RAD * 15.0);
equ.ra = RAD2HOUR * atan2(pos[1], pos[0]);
if (equ.ra < 0)
equ.ra += 24.0;
@@ -1866,6 +1866,7 @@ static const vsop_model_t vsop[] =
static void VsopCoords(const vsop_model_t *model, double t, double sphere[3])
{
int k, s, i;
double incr;
for (k=0; k < 3; ++k)
{
@@ -1879,9 +1880,12 @@ static void VsopCoords(const vsop_model_t *model, double t, double sphere[3])
for (i=0; i < series->nterms; ++i)
{
const vsop_term_t *term = &series->term[i];
sum += term->amplitude * cos(term->phase + (t * term->frequency));
sum += term->amplitude * cos(term->phase + (t * term->frequency));
}
sphere[k] += tpower * sum;
incr = tpower * sum;
if (k == LON_INDEX)
incr = fmod(incr, PI2); /* improve precision for longitudes, which can be hundreds of radians */
sphere[k] += incr;
tpower *= t;
}
}
@@ -1909,8 +1913,10 @@ static terse_vector_t VsopRotate(const double ecl[3])
static void VsopSphereToRect(double lon, double lat, double radius, double pos[3])
{
double r_coslat = radius * cos(lat);
pos[0] = r_coslat * cos(lon);
pos[1] = r_coslat * sin(lon);
double coslon = cos(lon);
double sinlon = sin(lon);
pos[0] = r_coslat * coslon;
pos[1] = r_coslat * sinlon;
pos[2] = radius * sin(lat);
}
@@ -2965,20 +2971,26 @@ astro_horizon_t Astronomy_Horizon(
astro_time_t *time, astro_observer_t observer, double ra, double dec, astro_refraction_t refraction)
{
astro_horizon_t hor;
double latrad, lonrad, decrad, rarad;
double uze[3], une[3], uwe[3];
double uz[3], un[3], uw[3];
double p[3], pz, pn, pw, proj;
double az, zd;
double spin_angle;
double sinlat = sin(observer.latitude * DEG2RAD);
double coslat = cos(observer.latitude * DEG2RAD);
double sinlon = sin(observer.longitude * DEG2RAD);
double coslon = cos(observer.longitude * DEG2RAD);
double sindc = sin(dec * DEG2RAD);
double cosdc = cos(dec * DEG2RAD);
double sinra = sin(ra * 15 * DEG2RAD);
double cosra = cos(ra * 15 * DEG2RAD);
latrad = observer.latitude * DEG2RAD;
lonrad = observer.longitude * DEG2RAD;
decrad = dec * DEG2RAD;
rarad = ra * HOUR2RAD;
double sinlat = sin(latrad);
double coslat = cos(latrad);
double sinlon = sin(lonrad);
double coslon = cos(lonrad);
double sindc = sin(decrad);
double cosdc = cos(decrad);
double sinra = sin(rarad);
double cosra = cos(rarad);
/*
Calculate three mutually perpendicular unit vectors
@@ -3090,7 +3102,7 @@ astro_horizon_t Astronomy_Horizon(
proj = sqrt(pr[0]*pr[0] + pr[1]*pr[1]);
if (proj > 0)
{
hor.ra = atan2(pr[1], pr[0]) * (RAD2DEG / 15.0);
hor.ra = RAD2HOUR * atan2(pr[1], pr[0]);
if (hor.ra < 0.0)
hor.ra += 24.0;
}
@@ -3098,7 +3110,7 @@ astro_horizon_t Astronomy_Horizon(
{
hor.ra = 0.0;
}
hor.dec = atan2(pr[2], proj) * RAD2DEG;
hor.dec = RAD2DEG * atan2(pr[2], proj);
}
}
generate/template/astronomy.cs
+29 -15
View File
@@ -1943,11 +1943,22 @@ $ASTRO_ADDSOL()
/// </summary>
public const double RAD2DEG = 57.295779513082321;
/// <summary>
/// The factor to convert radians to sidereal hours = 12/pi.
/// </summary>
public const double RAD2HOUR = 3.819718634205488;
/// <summary>
/// The factor to convert degrees to radians = pi/180.
/// </summary>
public const double DEG2RAD = 0.017453292519943296;
/// <summary>
/// The factor to convert sidereal hours to radians = pi/12.
/// </summary>
public const double HOUR2RAD = 0.2617993877991494365;
// Jupiter radius data are nominal values obtained from:
// https://www.iau.org/static/resolutions/IAU2015_English.pdf
// https://nssdc.gsfc.nasa.gov/planetary/factsheet/jupiterfact.html
@@ -2278,7 +2289,7 @@ $ASTRO_CSHARP_VSOP(Neptune)
}
}
private static double VsopFormulaCalc(vsop_formula_t formula, double t)
private static double VsopFormulaCalc(vsop_formula_t formula, double t, bool clamp_angle)
{
double coord = 0.0;
double tpower = 1.0;
@@ -2287,7 +2298,10 @@ $ASTRO_CSHARP_VSOP(Neptune)
double sum = 0.0;
foreach (vsop_term_t term in series.term)
sum += term.amplitude * Math.Cos(term.phase + (t * term.frequency));
coord += tpower * sum;
double incr = tpower * sum;
if (clamp_angle)
incr %= PI2; // improve precision: longitude angles can be hundreds of radians
coord += incr;
tpower *= t;
}
return coord;
@@ -2319,9 +2333,9 @@ $ASTRO_CSHARP_VSOP(Neptune)
double t = time.tt / DAYS_PER_MILLENNIUM; /* millennia since 2000 */
/* Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. */
double lat = VsopFormulaCalc(model.lat, t);
double lon = VsopFormulaCalc(model.lon, t);
double rad = VsopFormulaCalc(model.rad, t);
double lon = VsopFormulaCalc(model.lon, t, true);
double lat = VsopFormulaCalc(model.lat, t, false);
double rad = VsopFormulaCalc(model.rad, t, false);
/* Convert ecliptic spherical coordinates to ecliptic Cartesian coordinates. */
TerseVector eclip = VsopSphereToRect(lon, lat, rad);
@@ -2401,14 +2415,14 @@ $ASTRO_CSHARP_VSOP(Neptune)
double t = tt / DAYS_PER_MILLENNIUM; /* millennia since 2000 */
/* Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. */
double lat = VsopFormulaCalc(model.lat, t);
double lon = VsopFormulaCalc(model.lon, t);
double rad = VsopFormulaCalc(model.rad, t);
double lon = VsopFormulaCalc(model.lon, t, true);
double lat = VsopFormulaCalc(model.lat, t, false);
double rad = VsopFormulaCalc(model.rad, t, false);
TerseVector eclip_pos = VsopSphereToRect(lon, lat, rad);
double dlat_dt = VsopDerivCalc(model.lat, t);
double dlon_dt = VsopDerivCalc(model.lon, t);
double dlat_dt = VsopDerivCalc(model.lat, t);
double drad_dt = VsopDerivCalc(model.rad, t);
/* Use spherical coords and spherical derivatives to calculate */
@@ -3136,7 +3150,7 @@ $ASTRO_IAU_DATA()
}
else
{
ra = Math.Atan2(pos.y, pos.x) / (DEG2RAD * 15.0);
ra = RAD2HOUR * Math.Atan2(pos.y, pos.x);
if (ra < 0)
ra += 24.0;
@@ -3325,7 +3339,7 @@ $ASTRO_IAU_DATA()
case Body.Saturn:
case Body.Uranus:
case Body.Neptune:
return VsopFormulaCalc(vsop[(int)body].rad, time.tt / DAYS_PER_MILLENNIUM);
return VsopFormulaCalc(vsop[(int)body].rad, time.tt / DAYS_PER_MILLENNIUM, false);
default:
/* For non-VSOP objects, fall back to taking the length of the heliocentric vector. */
@@ -3589,8 +3603,8 @@ $ASTRO_IAU_DATA()
double coslon = Math.Cos(observer.longitude * DEG2RAD);
double sindc = Math.Sin(dec * DEG2RAD);
double cosdc = Math.Cos(dec * DEG2RAD);
double sinra = Math.Sin(ra * 15 * DEG2RAD);
double cosra = Math.Cos(ra * 15 * DEG2RAD);
double sinra = Math.Sin(ra * HOUR2RAD);
double cosra = Math.Cos(ra * HOUR2RAD);
// Calculate three mutually perpendicular unit vectors
// in equatorial coordinates: uze, une, uwe.
@@ -3679,7 +3693,7 @@ $ASTRO_IAU_DATA()
proj = Math.Sqrt(prx*prx + pry*pry);
if (proj > 0.0)
{
hor_ra = Math.Atan2(pry, prx) * (RAD2DEG / 15.0);
hor_ra = RAD2HOUR * Math.Atan2(pry, prx);
if (hor_ra < 0.0)
hor_ra += 24.0;
}
@@ -3687,7 +3701,7 @@ $ASTRO_IAU_DATA()
{
hor_ra = 0.0;
}
hor_dec = Math.Atan2(prz, proj) * RAD2DEG;
hor_dec = RAD2DEG * Math.Atan2(prz, proj);
}
}
else if (refraction != Refraction.None)
generate/template/astronomy.py
+20 -8
View File
@@ -56,6 +56,8 @@ CALLISTO_RADIUS_KM = 2410.3 #<const> The mean radius of Jupiter's moon Calli
_CalcMoonCount = 0
_RAD2HOUR = 3.819718634205488 # 12/pi = factor to convert radians to sidereal hours
_HOUR2RAD = 0.2617993877991494365 # pi/12 = factor to convert sidereal hours to radians
_DAYS_PER_TROPICAL_YEAR = 365.24217
_PI2 = 2.0 * math.pi
_EPOCH = datetime.datetime(2000, 1, 1, 12)
@@ -929,7 +931,7 @@ def _vector2radec(pos, time):
else:
dec = +90.0
else:
ra = math.degrees(math.atan2(pos[1], pos[0])) / 15.0
ra = _RAD2HOUR * math.atan2(pos[1], pos[0])
if ra < 0:
ra += 24
dec = math.degrees(math.atan2(pos[2], math.sqrt(xyproj)))
@@ -1234,11 +1236,16 @@ _vsop = [
$ASTRO_LIST_VSOP(Neptune),
]
def _VsopFormula(formula, t):
def _VsopFormula(formula, t, clamp_angle):
tpower = 1.0
coord = 0.0
for series in formula:
coord += tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
incr = tpower * sum(A * math.cos(B + C*t) for (A, B, C) in series)
if clamp_angle:
# Longitude angles can be hundreds of radians.
# Improve precision by keeping each increment within [-2*pi, +2*pi].
incr = math.fmod(incr, _PI2)
coord += incr
tpower *= t
return coord
@@ -1285,7 +1292,9 @@ def _VsopSphereToRect(lon, lat, rad):
def _CalcVsop(model, time):
t = time.tt / _DAYS_PER_MILLENNIUM
(lon, lat, rad) = [_VsopFormula(formula, t) for formula in model]
lon = _VsopFormula(model[0], t, True)
lat = _VsopFormula(model[1], t, False)
rad = _VsopFormula(model[2], t, False)
eclip = _VsopSphereToRect(lon, lat, rad)
return _VsopRotate(eclip).ToAstroVector(time)
@@ -1298,7 +1307,10 @@ class _body_state_t:
def _CalcVsopPosVel(model, tt):
t = tt / _DAYS_PER_MILLENNIUM
(lon, lat, rad) = [_VsopFormula(formula, t) for formula in model]
lon = _VsopFormula(model[0], t, True)
lat = _VsopFormula(model[1], t, False)
rad = _VsopFormula(model[2], t, False)
(dlon_dt, dlat_dt, drad_dt) = [_VsopDeriv(formula, t) for formula in model]
# Use spherical coords and spherical derivatives to calculate
@@ -1358,7 +1370,7 @@ def _VsopHelioDistance(model, time):
# The caller only wants to know the distance between the planet and the Sun.
# So we only need to calculate the radial component of the spherical coordinates.
# There is no need to translate coordinates.
return _VsopFormula(model[2], time.tt / _DAYS_PER_MILLENNIUM)
return _VsopFormula(model[2], time.tt / _DAYS_PER_MILLENNIUM, False)
def _CalcEarth(time):
return _CalcVsop(_vsop[Body.Earth.value], time)
@@ -2284,7 +2296,7 @@ def Horizon(time, observer, ra, dec, refraction):
latrad = math.radians(observer.latitude)
lonrad = math.radians(observer.longitude)
decrad = math.radians(dec)
rarad = math.radians(ra * 15.0)
rarad = ra * _HOUR2RAD
sinlat = math.sin(latrad)
coslat = math.cos(latrad)
@@ -2377,7 +2389,7 @@ def Horizon(time, observer, ra, dec, refraction):
pr = [(((p[j] - coszd0 * uz[j]) / sinzd0)*sinzd + uz[j]*coszd) for j in range(3)]
proj = math.sqrt(pr[0]*pr[0] + pr[1]*pr[1])
if proj > 0:
hor_ra = math.degrees(math.atan2(pr[1], pr[0])) / 15
hor_ra = _RAD2HOUR * math.atan2(pr[1], pr[0])
if hor_ra < 0:
hor_ra += 24
else:
generate/template/astronomy.ts
+38 -22
View File
@@ -42,15 +42,26 @@ export type FlexibleDateTime = Date | number | AstroTime;
export const KM_PER_AU = 1.4959787069098932e+8;
/**
* @brief The factor to convert radians to degrees = pi/180.
* @brief The factor to convert degrees to radians = pi/180.
*/
export const DEG2RAD = 0.017453292519943296;
/**
* @brief The factor to convert sidereal hours to radians = pi/12.
*/
export const HOUR2RAD = 0.2617993877991494365;
/**
* @brief The factor to convert degrees to radians = 180/pi.
* @brief The factor to convert radians to degrees = 180/pi.
*/
export const RAD2DEG = 57.295779513082321;
/**
* @brief The factor to convert radians to sidereal hours = 12/pi.
*/
export const RAD2HOUR = 3.819718634205488;
// Jupiter radius data are nominal values obtained from:
// https://www.iau.org/static/resolutions/IAU2015_English.pdf
// https://nssdc.gsfc.nasa.gov/planetary/factsheet/jupiterfact.html
@@ -1314,10 +1325,10 @@ function vector2radec(pos: ArrayVector, time: AstroTime): EquatorialCoordinates
return new EquatorialCoordinates(0, +90, dist, vec);
}
let ra = Math.atan2(vec.y, vec.x) / (DEG2RAD * 15);
let ra = RAD2HOUR * Math.atan2(vec.y, vec.x);
if (ra < 0)
ra += 24;
const dec = Math.atan2(pos[2], Math.sqrt(xyproj)) / DEG2RAD;
const dec = RAD2DEG * Math.atan2(pos[2], Math.sqrt(xyproj));
return new EquatorialCoordinates(ra, dec, dist, vec);
}
@@ -1381,8 +1392,8 @@ export function Horizon(date: FlexibleDateTime, observer: Observer, ra: number,
const coslon = Math.cos(observer.longitude * DEG2RAD);
const sindc = Math.sin(dec * DEG2RAD);
const cosdc = Math.cos(dec * DEG2RAD);
const sinra = Math.sin(ra * 15 * DEG2RAD);
const cosra = Math.cos(ra * 15 * DEG2RAD);
const sinra = Math.sin(ra * HOUR2RAD);
const cosra = Math.cos(ra * HOUR2RAD);
// Calculate three mutually perpendicular unit vectors
// in equatorial coordinates: uze, une, uwe.
@@ -1438,7 +1449,7 @@ export function Horizon(date: FlexibleDateTime, observer: Observer, ra: number,
if (proj > 0) {
// If the body is not exactly straight up/down, it has an azimuth.
// Invert the angle to produce degrees eastward from north.
az = -Math.atan2(pw, pn) * RAD2DEG;
az = -RAD2DEG * Math.atan2(pw, pn);
if (az < 0) az += 360;
} else {
// The body is straight up/down, so it does not have an azimuth.
@@ -1447,7 +1458,7 @@ export function Horizon(date: FlexibleDateTime, observer: Observer, ra: number,
}
// zd = the angle of the body away from the observer's zenith, in degrees.
let zd = Math.atan2(proj, pz) * RAD2DEG;
let zd = RAD2DEG * Math.atan2(proj, pz);
let out_ra = ra;
let out_dec = dec;
@@ -1466,14 +1477,14 @@ export function Horizon(date: FlexibleDateTime, observer: Observer, ra: number,
}
proj = Math.sqrt(pr[0]*pr[0] + pr[1]*pr[1]);
if (proj > 0) {
out_ra = Math.atan2(pr[1], pr[0]) * RAD2DEG / 15;
out_ra = RAD2HOUR * Math.atan2(pr[1], pr[0]);
if (out_ra < 0) {
out_ra += 24;
}
} else {
out_ra = 0;
}
out_dec = Math.atan2(pr[2], proj) * RAD2DEG;
out_dec = RAD2DEG * Math.atan2(pr[2], proj);
}
}
@@ -1746,7 +1757,7 @@ export function GeoMoon(date: FlexibleDateTime): Vector {
return new Vector(mpos2[0], mpos2[1], mpos2[2], time);
}
function VsopFormula(formula: any, t: number): number {
function VsopFormula(formula: any, t: number, clamp_angle: boolean): number {
let tpower = 1;
let coord = 0;
for (let series of formula) {
@@ -1754,7 +1765,10 @@ function VsopFormula(formula: any, t: number): number {
for (let [ampl, phas, freq] of series) {
sum += ampl * Math.cos(phas + (t * freq));
}
coord += tpower * sum;
let incr = tpower * sum;
if (clamp_angle)
incr %= PI2; // improve precision for longitudes: they can be hundreds of radians
coord += incr;
tpower *= t;
}
return coord;
@@ -1800,18 +1814,20 @@ function VsopRotate(eclip: ArrayVector): TerseVector {
function VsopSphereToRect(lon: number, lat: number, radius: number): ArrayVector {
// Convert spherical coordinates to ecliptic cartesian coordinates.
const r_coslat = radius * Math.cos(lat);
const coslon = Math.cos(lon);
const sinlon = Math.sin(lon);
return [
r_coslat * Math.cos(lon),
r_coslat * Math.sin(lon),
r_coslat * coslon,
r_coslat * sinlon,
radius * Math.sin(lat)
];
}
function CalcVsop(model: any[], time: AstroTime): Vector {
const t = time.tt / DAYS_PER_MILLENNIUM; // millennia since 2000
const lon = VsopFormula(model[LON_INDEX], t);
const lat = VsopFormula(model[LAT_INDEX], t);
const rad = VsopFormula(model[RAD_INDEX], t);
const lon = VsopFormula(model[LON_INDEX], t, true);
const lat = VsopFormula(model[LAT_INDEX], t, false);
const rad = VsopFormula(model[RAD_INDEX], t, false);
const eclip = VsopSphereToRect(lon, lat, rad);
return VsopRotate(eclip).ToAstroVector(time);
}
@@ -1820,9 +1836,9 @@ function CalcVsopPosVel(model: any[], tt: number): body_state_t {
const t = tt / DAYS_PER_MILLENNIUM;
// Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates.
const lon = VsopFormula(model[LON_INDEX], t);
const lat = VsopFormula(model[LAT_INDEX], t);
const rad = VsopFormula(model[RAD_INDEX], t);
const lon = VsopFormula(model[LON_INDEX], t, true);
const lat = VsopFormula(model[LAT_INDEX], t, false);
const rad = VsopFormula(model[RAD_INDEX], t, false);
const dlon_dt = VsopDeriv(model[LON_INDEX], t);
const dlat_dt = VsopDeriv(model[LAT_INDEX], t);
@@ -2422,7 +2438,7 @@ export function HelioVector(body: Body, date: FlexibleDateTime): Vector {
export function HelioDistance(body: Body, date: FlexibleDateTime): number {
const time = MakeTime(date);
if (body in vsop) {
return VsopFormula(vsop[body][RAD_INDEX], time.tt / DAYS_PER_MILLENNIUM);
return VsopFormula(vsop[body][RAD_INDEX], time.tt / DAYS_PER_MILLENNIUM, false);
}
return HelioVector(body, time).Length();
}
@@ -5663,7 +5679,7 @@ function GeoidIntersect(shadow: ShadowInfo): GlobalSolarEclipseInfo {
// Adjust longitude for Earth's rotation at the given UT.
const gast = sidereal_time(peak);
longitude = (RAD2DEG * Math.atan2(py, px) - (15*gast)) % 360.0;
longitude = (RAD2DEG*Math.atan2(py, px) - (15*gast)) % 360.0;
if (longitude <= -180.0) {
longitude += 360.0;
} else if (longitude > +180.0) {