<p>There are two things I can see in this code:</p>

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

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

<p><em>edit</em>
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.</p>

<p>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 <code>LHA = GMST * 360/24 + alpha</code> should be <code>longitude = GMST*360/24 - alpha</code>. (and return <code>longitude</code>) 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.</p>

<p>You seem to be attempting a very accurate ephemerides for the sun. You can check you accuracy against <a href="http://rhodesmill.org/pyephem/index.html" rel="nofollow">PyEphem</a>, or the <a href="http://ssd.jpl.nasa.gov/horizons.cgi" rel="nofollow">Nasa Horizon ephemerides</a>. </p>

CC BY-SA 3.0
2015-11-05T15:50:18.650Z
Thanks, I made these corrections, very helpful. However, the results still seem wrong.
Tim Dierks
CC BY-SA 3.0
2015-11-05T18:35:50.863Z
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
James K
CC BY-SA 3.0
2015-11-06T02:26:39.267Z
Thanks, I'll look at those. Shouldn't I expect GMST to be zero on the spring equinox?
Tim Dierks
CC BY-SA 3.0
2015-11-04T19:47:53.567Z
2015-11-05T18:51:06.643Z
2015-11-05T18:51:06.643Z
James K
<p>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.)</p>

<p>Here is my python code:</p>

<pre><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)
</code></pre>

<p>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.</p>

<p>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?</p>

<p>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.</p>

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

<pre><code>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
</code></pre>

<p>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.</p>

<p>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.</p>

<p>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.</p>

<p>Geocentric:</p>

<pre><code> 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
</code></pre>

<p>Greenwich:</p>

<pre><code> 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
</code></pre>

<p>Thanks again for any help you can provide.</p>

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

CC BY-SA 3.0
2015-11-04T19:12:59.883Z
Could you show us the output too?
Eubie Drew
CC BY-SA 3.0
2015-11-04T19:21:03.727Z
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)`
James K
CC BY-SA 3.0
2015-11-04T19:25:47.300Z
@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.
Eubie Drew
CC BY-SA 3.0
2015-11-04T19:38:03.027Z
Next thing that seems weird is that "alpha" is being calcuated in radians, instead of the customary Hours Min Sec. and  lambda is in degrees.
James K
CC BY-SA 3.0
2015-11-05T14:59:48.410Z
Thanks for all the comments. @barrycarter -- It prints the values right near the end.
Tim Dierks
CC BY-SA 3.0
2015-11-05T15:04:19.133Z
@james-kilfiger: thanks; I was making that correction for the sidereal time calculation (see n - 0.5), but for some reason I was under the impression that the Julian days used for the other calculations started at midnight. And, of course you're right about alpha!
Tim Dierks
<the-sun><earth><coordinate><software><orbital-elements>
Math for calculating the terrestrial longitude directly under the sun with time
Tim Dierks
# Math for calculating the terrestrial longitude directly under the sun with time

(Asked by: username_0 on 2015-11-04)

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]

**Tags:** the-sun, earth, coordinate, software, orbital-elements

> **Comments:**
> **username_1**: Could you show us the output too?  
> **username_2**: 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)`  
> **username_1**: @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.  
> **username_2**: Next thing that seems weird is that "alpha" is being calcuated in radians, instead of the customary Hours Min Sec. and  lambda is in degrees.  
> **username_0**: Thanks for all the comments. @barrycarter -- It prints the values right near the end.  
> **username_0**: @james-kilfiger: thanks; I was making that correction for the sidereal time calculation (see n - 0.5), but for some reason I was under the impression that the Julian days used for the other calculations started at midnight. And, of course you're right about alpha!  

---

## 1 answers

### Answer by username_2 on 2015-11-04. Score: 2

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](http://rhodesmill.org/pyephem/index.html), or the [Nasa Horizon ephemerides](http://ssd.jpl.nasa.gov/horizons.cgi).

> **Comments:**
> **username_0**: Thanks, I made these corrections, very helpful. However, the results still seem wrong.  
> **username_2**: 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  
> **username_0**: Thanks, I'll look at those. Shouldn't I expect GMST to be zero on the spring equinox?  

---

