[AstroPy] sun angles from astropy deviate from US Navy algorithm

[Koot, Sjaak]

Airbus Amber
Hi,

I am comparing the sun angles computed with astropy V5.3.3 with the algorithm as described
here https://aa.usno.navy.mil/faq/sun_approx.

According to that website:

Given below is a simple algorithm for computing the Sun's angular coordinates to an accuracy of about 1 arcminute within two centuries of 2000.

So my understanding is that the algorithm should not deviate (much) more than 0.0166666667 degrees for dates from 2000 - 2023.

However when I plot the delta's for the Ra and dec from 2000 - 2023 I get the following results:

[cid:image001.png at 01D9EA0F.442D2A80]


[cid:image002.png at 01D9EA0F.442D2A80]

I have no astronomical background whatsoever so from the plots I do not understand what is happening.

I have attached the python3 code (as .txt files) generating the data.

Do you have any idea what is happening?

Best regards,
Sjaak


M   +31 (0)6 53 79 80 30
E    j.koot at airbusDS.nl



-- -----------------------------------------------------------------------
Airbus Netherlands B.V. te Leiden. KvK nummer: 28086907.
Airbus Netherlands B.V., Leiden, The Netherlands.
Chamber of Commerce number 28086907.
-- -----------------------------------------------------------------------
This communication is intended for use by the addressee and may
contain confidential or privileged information. If you receive this
communication unintentionally, please notify us immediately and
delete the message from your computer without making any copies.
-- -----------------------------------------------------------------------
Please consider the environment before printing this email
-- -----------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20230918/093cecc8/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 36416 bytes
Desc: image001.png
URL: <https://mail.python.org/pipermail/astropy/attachments/20230918/093cecc8/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image002.png
Type: image/png
Size: 54722 bytes
Desc: image002.png
URL: <https://mail.python.org/pipermail/astropy/attachments/20230918/093cecc8/attachment-0003.png>
-------------- next part --------------
#!/usr/bin/env python3
# https://aa.usno.navy.mil/faq/sun_approx

from math import sin, cos, asin, atan2, degrees, radians, fmod

print("{:11s},{:7s},{:7s}".format('Time', 'Ra', 'dec'))

JD = 2451544.50
#JD = 2459945.50

JD_0 = JD
JDD  = JD - JD_0

while JDD < 20*365:

    D = JD - 2451545.0 
    
    g = 357.529 + 0.98560028 * D
    q = 280.459 + 0.98564736 * D
    L = q + 1.915 * sin(radians(g)) + 0.020 *  sin(radians(2*g))
    
    R = 1.00014 - 0.01671 * cos(radians(g)) - 0.00014 * cos(radians(2*g))
    e = 23.439 - 0.00000036 * D
    
    
    #tan RA = cos e sin L / cos L
    RA = degrees(atan2(cos(radians(e))*sin(radians(L)), cos(radians(L))))
    #sin d = sin e sin L 
    d = degrees(asin(sin(radians(e))*sin(radians(L))))

    #qm = fmod(q, 360)
    #EqT = (qm - RA)
    #if EqT > 180.0 : EqT -= 360.0
    ## From degrees to minutes
    #EqT = 60*EqT/15 

    if RA < 0 : RA += 360

    #print(f"JD = {JD:6.2f} D = {D:6.2f} qm = {qm:6.2f} EqT = {EqT:5.2f} RA = {RA:6.2f} d = {d:6.2f}")
    print(f"{JD:.2f},{RA:6.2f},{d:6.2f}")

    JD += 1
    JDD = JD - JD_0
-------------- next part --------------
#!/usr/bin/env python3
from astropy.coordinates import get_sun
import astropy.coordinates as coord
from astropy.time import Time
import astropy.units as u
import time

my_time = Time('2000-01-01T00:00:00', format='isot', scale='utc')
#my_time = Time('2459945.5', format='jd')
my_time.format='jd'

dt = 0

print("{:11s},{:7s},{:7s}".format('Time', 'Ra', 'dec'))

while dt < 20*365:

    radec = coord.GCRS()
    sun_rade  = get_sun(my_time).transform_to(radec)

    print("{:.2f},{:6.2f},{:6.2f}".format(my_time.jd, sun_rade.ra.degree, sun_rade.dec.degree))

    my_time = my_time + 1 * u.day
    dt = dt + 1