How do I calculate the closest point to the sun on Earth's surface?

Question

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.

Answer 1

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.