3
$\begingroup$

I'm trying to calculate the longitude on the earth where it's noon at some time. (That is, the longitude which is coplanar with the plane defined by the sun and the earth's axis.)

Here is my python code:

from math import sin, cos, tan, atan, pi
def sun_longitude(when):
    """Given time in ms since 1/1/1970, return the longitude the sun is over at that moment"""
    # https://en.wikipedia.org/wiki/Position_of_the_Sun
    jdn = 2440587.5 + when / (1000.0 * 3600 * 24)
    n = jdn - 2451545.0 # 1/1/2000
    L = (280.460 + 0.9856474 * n) % 360.0
    g = (357.528 + 0.9856004 * n) % 360.0
    degtorad = 2.0 * pi / 360.0
    lambda_ = L + 1.915 * sin(g * degtorad) + 0.020 * sin(2 * g * degtorad)

    # https://en.wikipedia.org/wiki/Axial_tilt#Short_term
    T = n / (365 * 100.0) # julian centuries from 2000
    # ε = 23° 26′ 21.406″ − 46.836769″ T − 0.0001831″ T2
    epsilon = 23.4392794 - 0.780612817 * T - 5.0861e-8 * T * T;
    alpha = atan(cos(epsilon * degtorad) * tan(lambda_ * degtorad)) / degtorad

    # https://en.wikipedia.org/wiki/Sidereal_time
    # https://en.wikipedia.org/wiki/Hour_angle
    GMST = (18.697374558 + 24.06570982441908 * n) % 24.0
    print "jdn = {jdn}, n = {n}, L = {L}, g = {g}, lambda_ = {lambda_}, T = {T}, epsilon = {epsilon}, alpha = {alpha}, GMST = {GMST}".format(**locals())
    LHA = GMST * 360/24 + alpha;

    return LHA

from datetime import datetime
def sun_long_from_str(whenstr):
    when = datetime.strptime(whenstr, "%Y-%m-%d %H:%M")
    secs = (when - datetime(1970, 1, 1)).total_seconds()
    return sun_longitude(secs * 1000.0)

I'm aware that I'm not correcting for the correct quandrant in figuring alpha, but I have errors elsewhere; for example, for noon GMT Jan 1, 2000, I'm getting 279° rather than close to 0, which is what I'd expect.

This isn't giving me correct answers, and I'm at a bit of a loss for how to debug it. Can anyone find my mistakes or point me to some reasonable sample code or worked-through example for this?

Because, once I get the algorithm right, I'll be reimplementing for an embedded device, I can't just use an implementing package, and I haven't found a library which is straightforward enough for me to understand how I can implement this specific question.

Thanks for a couple of error corrections. Now, when I run it, I get the following for these examples:

2000-01-01 12:00:
  jdn = 2451545.0, n = 0.0, L = 280.46, g = 357.528, lambda_ = 280.375680197, T = 0.0, epsilon = 23.4392794, alpha = -78.7141369122, GMST = 18.697374558
  result: 201.74648145775257

2000-03-20 07:35:
  jdn = 2451623.81597, n = 78.815972222, L = 358.144758099, g = 75.2090537484, lambda_ = 360.006175505, T = 0.00215934170471, epsilon = 23.4375937902, alpha = 0.00566598789588, GMST = 19.4596915825
  291.9010397253072

As you can see, the first query (noon on Jan 1, 2000) now has a correct Julian day number. At noon GMT, I would expect the sun to be above a longitude near 0, so 201 is incorrect.

The second query is the time of the spring equinox in 2000, which I expected to result in zeroing out the sidereal part of the calculation, but it does not.

I attempted to consult the NASA HORIZONS web interface for ephemerides for the sun on 1/1/2000 at noon GMT, for both "Geocentric" and Greenwich locations, but I don't know how to evaluate the result and compare it to my work above.

Geocentric:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00     18 45 09.36 -23 01 59.7 -26.78 -10.59 0.98332762653520  -0.0127281   0.0000 /?   0.0000

Greenwich:

 Date__(UT)__HR:MN     R.A._(ICRF/J2000.0)_DEC  APmag  S-brt            delta      deldot    S-O-T /r    S-T-O
**************************************************************************************************************
 2000-Jan-01 12:00 *m  18 45 09.36 -23 02 08.3 -26.78 -10.59 0.98331613178086  -0.0166655   0.0000 /?   0.0000

Thanks again for any help you can provide.

[Edited to correct Julian day calculation and use of degrees vs. radians for alpha and to add examples]

$\endgroup$
  • $\begingroup$ General debugging suggestion: print out the variables after you calculate them to see where you're going wrong? Also, I think you can do this just using the Equation of Time and the fact the Earth rotates every 24 hours (with respect to the Sun) $\endgroup$ – barrycarter Nov 4 '15 at 18:14
  • $\begingroup$ Could you show us the output too? $\endgroup$ – Aabaakawad Nov 4 '15 at 19:12
  • $\begingroup$ First thing I noticed is that Unix time is measured in seconds from midnight, whereas Julian dates are measured from midday. Try jdn = 2440587.5 + when / (1000.0 * 3600 * 24) $\endgroup$ – James K Nov 4 '15 at 19:21
  • $\begingroup$ @barrycarter it depends on how accurate you want to be. Because the Earth is in an elliptical orbit (changing speed), the literal solar day varies slightly from 24 hours according to where you are in the orbit. The solar day does average precisely 24 hours. $\endgroup$ – Aabaakawad Nov 4 '15 at 19:25
  • 1
    $\begingroup$ @Aabaakawad Yes, that's why I mentioned the Equation of Time (as a correction to the 24-hour rotation cycle w/r/t the sun) $\endgroup$ – barrycarter Nov 5 '15 at 2:48
2
$\begingroup$

There are two things I can see in this code:

First, the Julian date is measured from Midday, whereas the unix epoch is measured from midnight. jdn = 2440587.5 + when / (1000.0 * 3600 * 24) should be the correct expression.

Secondly, alpha = atan(cos(epsilon * degtorad) * tan(lambda_ * degtorad)) calculates a right ascension in radians, you should convert this to degrees for the last part of the calculations.

edit You do need to fix the quadrant, for example at when n==0, you should get a right ascension of 281.285. However the last error is in the calcualtion of the longitude.

The relevant equation is $LHA = GMST -\lambda - \alpha$. When the sun is on the meridian, LHA = 0, and (rearranging) $\lambda = GMST -\alpha$. Thus the line LHA = GMST * 360/24 + alpha should be longitude = GMST*360/24 - alpha. (and return longitude) If you do this with n=0, you get a longitude of -1.5 degrees. The sun is at meridian at noon on Jan 1 2000, if you are 1.5 degrees west of Greenwich.

You seem to be attempting a very accurate ephemerides for the sun. You can check you accuracy against PyEphem, or the Nasa Horizon ephemerides.

$\endgroup$
  • $\begingroup$ Thanks, I made these corrections, very helpful. However, the results still seem wrong. $\endgroup$ – Tim Dierks Nov 5 '15 at 15:50
  • $\begingroup$ I think the error is in the calcualtion of LHA, which should be LHA = GMST - alpha. Also the value of alpha should be in the range 0--360, whereas atan returns in the range -90--90, you will need to think about the correcting the quadrant $\endgroup$ – James K Nov 5 '15 at 18:35
  • $\begingroup$ Thanks, I'll look at those. Shouldn't I expect GMST to be zero on the spring equinox? $\endgroup$ – Tim Dierks Nov 6 '15 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.