5
$\begingroup$

Couldn't find an answer so here I come:

On 21st of June and 21st of December these points would lie on the tropics.

But what about the other days of the year?

Is there a simple formula to calculate the latitude, longitude of the closest point to the Sun on Earth's surface at any given date-time? And would the shadow circle of the earth be centered at the opposite point?

Edit: One part of the answer seems to be here

Edit2: After more searching, i finally stumbled upon a fully cooked up solution

However, when i compare the result to, for example https://www.timeanddate.com/worldclock/sunearth.php, they do differ quite significantly. E.g., the website reports

On Thursday, 25 August 2022, 06:58:00 UTC the Sun is at its zenith at Latitude: 10° 44' North, Longitude: 76° 03' East

while the function generates

Lat: 19.63448840599144 Lon: 77.26722364174202.

So longitude is pretty close, latitude pretty far off.

For completeness, here's the function as i translated:

function toInt(d) {
    return d | 0
}
function subsolar(ye, mo, da, ho, mi, se) {
    var ta = Math.PI * 2.0
    var ut = ho + mi / 60.0 + se / 3600.0
    var t = 367.0 * ye - toInt(7.0 * (ye + toInt((mo + 9.0) / 12.0)) / 4.0)
    var dn = t + toInt((275.0 * mo) / 9.0) + da - 730531.5 + ut / 24.0
    var sl = dn * 0.01720279239 + 4.894967873
    var sa = dn * 0.01720197034 + 6.240040768
    t = sl + 0.03342305518 * Math.sin(sa)
    var ec = t + 0.0003490658504 * Math.sin(2.0 * sa)
    var ob = 0.4090877234 - 0.000000006981317008 * dn
    var st = 4.894961213 + 6.300388099 * dn
    var ra = Math.atan2(Math.cos(ob) * Math.sin(ec), Math.cos(ec))
    var de = Math.asin(Math.sin(ob) * Math.sin(ec))
    var la = degrees(de)
    var lo = degrees(ra - st) % 360.0
    lo = (lo > 180.) ? lo - 360. : lo
    lo = (lo < -180.) ? lo + 360. : lo
    return [la, lo]
}

To be fair, this algorithm does not look the same as what previously published, like here, which is instead what wikipedia references. However, 3 is unclear at best, not specifying

  1. what Nleap is and how to derive
  2. How to derive Fraction of day (i would assume as ((h + (m/60.0) + (s/3600.0)) / 24.0) , but i can't be sure)
  3. what Tgmt is and in what units it is expected to be given (hours? minutes? fraction of the day?)

Edit 3: the python algorithm (translated above into javascript) actually works reasonably well. The issue i was having is that javascript getMonth returns values in the [0-11] range, while this method expects values in the [1-12] range.

$\endgroup$
11
  • 3
    $\begingroup$ It depends on your definition of simple and your math skills. There is a "simple" equation to calculate the Right Ascension of the Sun, a simple equation to calculate the Declination of the Sun, and a simple formula to calculate the Greenwich Mean Sidereal Time (GMST). Knowing those three values, the subsolar point can be calculated. $\endgroup$
    – JohnHoltz
    Commented Aug 21, 2022 at 23:35
  • 1
    $\begingroup$ "And would the shadow circle of the earth be centered at the opposite point?" Draw or imagine a circle for Earth with a dot at Earth's center. Then draw/imagine a long straight line passing through the center. The two places where the line intersects the surface are indeed the closest point to the Sun and the center of Earths' shadow circle (actually a cylinder). Of course this model assumes a point Sun and spherical Earth but I think that's okay assuming we're looking for a "simple formula". In order to keep it simple, we also need Earth's orbit around the Sun to be circular as well. $\endgroup$
    – uhoh
    Commented Aug 22, 2022 at 0:25
  • 3
    $\begingroup$ BTW, this is called the subsolar point. Do you specifically want to calculate it yourself? You can get it from various websites, eg Horizons $\endgroup$
    – PM 2Ring
    Commented Aug 22, 2022 at 4:08
  • 2
    $\begingroup$ @PM2Ring thanks for the excellent pointer! one question, as layman: what exactly is jdn argument, j2000, (julian centuries, i would assume?) and how to derive these values? Probably they are what we call normally something else, but i just wanted to double check. Also, i get name 'poly' is not defined. What would anyway be its javascript equivalent? $\endgroup$
    – Pa_
    Commented Aug 22, 2022 at 6:40
  • 2
    $\begingroup$ @PM2Ring found this resource equation-of-time.info/calculating-the-equation-of-time , would that be about right? (seems easier to translate into JS or other languages than porting poly()) $\endgroup$
    – Pa_
    Commented Aug 22, 2022 at 7:06

1 Answer 1

3
+100
$\begingroup$

In celestial navigation this point is called the Geographic Position (GP) of the Sun. Since it is also useful to compute the GP of other bodies, I will give a general solution that applies to all RA/Dec coordinates.

Given the Right Ascension and Declination of an object, the GP of the body is:

$$ \begin{align*} Latitude &= Declination \\ Longitude &= 360^\circ - (360^\circ - RA + GST) \end{align*} $$

Where GST in the Greenwich Sidereal Time. For high accuracy, it is necessary to use Greenwich Apparent Sidereal Time, but since you asked for a simple formula, I will use Greenwich Mean Sidereal Time, which will still be accurate for recent times.

$$ \begin{align*} \theta &= 0.7790572732640 + 0.00273781191135448D_u + frac(JD(UT1)) \\ T &=\frac{JD_{UT1} - 2451545.0}{36525} \\ GMST &= \theta + 0.014506 + 4612.15739966T + 1.39667721T^2 + -0.00009344T^3 + 0.00001882T^4 \\ \end{align*} $$

$ D_u $ = Number of WHOLE days of UT1 since J2000 (JD(UT1) – 2451545.0)

$ \theta $ = Earth Rotation Angle (in arcseconds)

T = Centuries of UT1 since J2000

GMST = Greenwich Mean Sidereal Time (in arcseconds)

The last piece of information needed is the RA/Dec of the Sun, the Astronomical Almanac provides probably the shortest useful algorithm for the position of the Sun, valid from 1950-2050. If you want higher accuracy, or over a longer period of time, you'll need to use an ephemeris such as VSOP87. Here is the algorithm implemented in JavaScript, and you can test it here

//Low precision sun position from Astronomical Almanac page C5 (2017 ed).
//Accuracy 1deg from 1950-2050
function sunPosition(jd)    {
    const torad=Math.PI/180.0;
    n=jd-2451545.0;
    L=(280.460+0.9856474*n)%360;
    g=((375.528+.9856003*n)%360)*torad;
    if(L<0){L+=360;}
    if(g<0){g+=Math.PI*2.0;}

    lamba=(L+1.915*Math.sin(g)+0.020*Math.sin(2*g))*torad;
    beta=0.0;
    eps=(23.439-0.0000004*n)*torad;
    ra=Math.atan2(Math.cos(eps)*Math.sin(lamba),Math.cos(lamba));
    dec=Math.asin(Math.sin(eps)*Math.sin(lamba));
    if(ra<0){ra+=Math.PI*2;}
    return [ra/torad/15.0,dec/torad];
}

A higher accuracy version would use the full reduction process explained in the Explanatory Supplement for the Astronomical Almanac which corrects for precession, nutation, aberration, parallax, polar motion, time scales, and more.

The Nautical Almanac enumerates the values with all of the corrections for several stars, the sun and planets with an accuracy of .1 arcminutes.

I have implemented this algorithm here, though it is not complete (I intend to add corrections for some of the effects producing an error greater than .1 arcminutes). But it currently implements all of the parts mentioned above.

$\endgroup$
7
  • 1
    $\begingroup$ The equation for the latitude can obviously be simplified, but I chose to show it as if converted directly from celestial navigation terminology where hour angles are used which are opposite from RA and longitude. If you look at texts on celestial navigation, the way I've written it makes it more obvious as to what is going on. $\endgroup$ Commented Aug 25, 2022 at 14:10
  • 1
    $\begingroup$ Thanks for the nice answer! however, JD derivation technically is missing :-) $\endgroup$
    – Pa_
    Commented Aug 25, 2022 at 20:10
  • 1
    $\begingroup$ I found JD derivation method in your website. However, your solution also seem to produce different results compared to timeanddate.com/worldclock/sunearth.php , can you double check? $\endgroup$
    – Pa_
    Commented Aug 25, 2022 at 20:24
  • 1
    $\begingroup$ Right now, timeanddate reports On Friday, 26 August 2022, 15:02:00 UTC the Sun is at its zenith at Latitude: 10° 16' North, Longitude: 45° 03' West, your method [8.360038700401809,19.460248533502863], but in the evening the latitude error is greater even $\endgroup$
    – Pa_
    Commented Aug 26, 2022 at 15:05
  • 1
    $\begingroup$ thanks, i had to double check, and i was getting the julian date calculation (not included in the snipped here) wrong. However, the longitude straight out of the function does not seem right.. I now finally found a solution (edited the question), but If you can please fix the answer including all that is necessary to get the right values out of the method starting with gregorian date, i will still accept the answer. $\endgroup$
    – Pa_
    Commented Aug 26, 2022 at 21:22

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .