NavList:

   

Reply

I'm starting a new thread since the Julian date formula formerly in
question now seems accurate. And, after the bug in the Copilot GHA
formula is fixed (it subtracted ecliptical longitude instead of RA from
GMST) the Sun coordinates are accurate too. In a Monte Carlo test I get
about 1.1′ and 0.3′ root mean square error for GHA and dec respectively
at random times in the 21st century.

Some simplification is possible. In step 7, true longitude is corrected
to apparent longitude:

λ = θ - 0.00569 - 0.00478*SIN(RADIANS(125.04 - 1934.136*T))

Well, you can omit that. Just do λ = θ. I find RMS errors in the 21st
century practically the same: 1.3′ and 0.3′.

Here's the corrected and simplified code in the C# language. Like all
languages in the C family, the keyword "double" declares a double
precision (64 bit) floating point variable, e.g., "double L0". The
ATAN2(y, x) function considers its arguments (note the order) as
Cartesian coordinates and returns the corresponding angle. C# has no ^
exponentiation operator, so expressions have been factored.

In variable names C# allows Unicode characters, such as Greek letters.
Normally I don't do that, but in this case I've retained the same names
as the code posted by Jim.

All angles are radians. If the formula gives its result in degrees I
immediately convert to radians. Also, I don't normalize angles to a
conventional range such as 0 - 360.

The time argument is UT1. To practical accuracy we can assume UT1 = UTC.
To be strictly correct, the GMST expression should use UT1, and
everything else TT. However, the error from ignoring this 69 second
difference is only a few 0.01′


// Julian century
double T = (UT1 - 2451545.0) / 36525.0;

// Geometric Mean Longitude
double L0 = RADIANS(280.46646 + (36000.76983 + T * 0.0003032) * T);

// Mean Anomaly
double M = RADIANS(357.52911 + (35999.05029 - T * 0.0001537) * T);

// Equation of Center
double C = RADIANS(SIN(M) * 1.914602 - SIN(2 * M) * 0.004817 + SIN(3 *
M) * 0.000014);

// True Longitude
double θ = L0 + C;

// Apparent Longitude
double λ = θ;

// Obliquity of the Ecliptic
double ε = RADIANS(23.439291 - 0.0130042 * T);

// RA
double α = ATAN2(COS(ε) * SIN(λ), COS(λ));

// Declination
DEC = ASIN(SIN(ε) * SIN(λ));

// Greenwich Mean Sidereal Time
double GMST = RADIANS(280.46061837 + 360.98564736629 * (UT1 - 2451545) +
(0.000387933 -
T / 38710000) * T * T);

// GHA
GHA = GMST - α;


More accuracy (0.3′ GHA, 0.05′ dec) is possible with the formulas from
Astronomical Algorithms (chapter 24) by Meeus.


// centuries since J2000
double T = (UT1 - 2451545.0) / 36525.0;

// geometric mean longitude wrt mean equinox
double L0 = RADIANS(280.46645 + (36000.76983 + .0003032 * T) * T);

// mean anomaly
double M = RADIANS(357.52910 + (35999.05030 + (-.0001559 - 4.8e-7 * T) *
T) * T);

// equation of center
double C = RADIANS(
(1.914600 + (-.004817 - 1.4e-5 * T) * T) * SIN(M)
+ (.019993 - 1.01e-4 * T) * SIN(2 * M)
+ 2.90e-4 * SIN(3 * M));

// true longitude
double θ = L0 + C;

// longitude of Moon mean orbit ascending node
double Ω = RADIANS(125.04 - 1934.136 * T);

// apparent longitude
double λ = θ + RADIANS(-0.00569 - 0.00478 * SIN(Ω));

// obliquity of the ecliptic
double ε = RADIANS(23.439291 - 0.0130042 * T + .00256 * COS(Ω));

// RA
double α = ATAN2(COS(ε) * SIN(λ), COS(λ));

// Greenwich Mean Sidereal Time
double GMST = RADIANS(280.46061837 + 360.98564736629 * (UT1 - 2451545) +
(0.000387933 -
T / 38710000) * T * T);

GHA = GMST - α;
DEC = ASIN(SIN(ε) * SIN(λ));

--
Paul Hirose
sofajpl.com

   

Reply