Location of the point on Earth with the Sun at the zenith

Question

This should be an easy question, but I don't know the correct syntax. I'm trying to create an animation based on NightHemisphere during a year. For each day, the instant correspond to the time of the Sunset at Paris Looking on previous answers, I have constructed this function that works:

f := Function[x, DateObject[{2022, 1, 1}] + Quantity[(x - 1), "Days"]];
tg = Table[Show[
    GeoGraphics[
     NightHemisphere[
      Sunset[GeoPosition[{48.858042`, 2.2910492`}], DateObject[f[n]], 
       TimeZone -> +1]], GeoRange -> "World", 
     GeoCenter -> GeoPosition[{0, 0}]],
    ImageSize -> Large], {n, 1, 365}];
Export["sunset_paris.gif", tg]

Now I want to add a disk representing the position of the Earth place just below the Sun at that time, as it appears in, for instance, in https://www.timeanddate.com/worldclock/sunearth.html but I don't know how to transform SunPosition[] or other function to the position of this place. SunPosition gives the astronomical position of the Sun as seen from a location, while I'm looking for a point on Earth. Any help is welcome.

Answer 1

First provide a latitude at which you want your disk (representing the sun) to locate on your plots - let say at $x^\circ\;\mathrm{N}$. And I'm going to define a function:

zenithPole[n_] := QuantityMagnitude[
    First@SunPosition[GeoPosition[{90, 0}], DateObject[{2022, 1, n}, TimeZone -> 1]]
]

(*at 2022-01-01 00:00*)
zenithPole[1]

(* 344.18 *)

which outputs the azimuth angle of the sun, seen from the north pole GeoPosition[{90, 0}], at time {DateObject[{2022, 1, n}, TimeZone -> 1}. And I will subtract this from 360 (*degree*), then that's the longitude of every geoposition in which the zenith is. To see for yourself if this is true,

zenith[x_,n_] := SunPosition[
    GeoPosition[{x, 360 - zenithPole[n]}]
, DateObject[{2022, 1, n}, TimeZone -> 1]]

(*at 2022-01-01 00:00*)
zenith[0,1] (*at the equator*)
zenith[2.2910492`,1] (*at E 2.2910492`*)

(*
  {180.00°, -66.98°}
  {180.00°, -69.27°}
*)

As you can see, the sun is right at the zenith. So the location you want is GeoPosition[{x, 360 - zenithPole[n]}], where the x is the latitude where you want your disk to be at.

Answer 2

We can compute those locations as follows:

subSolarPoint[date_] := GeoPosition[Reverse[SunPosition[date, CelestialSystem -> "Equatorial"] - {SiderealTime[GeoPosition[{0, 0}], date], 0}]]

subLunarPoint[date_] := GeoPosition[Reverse[MoonPosition[date, CelestialSystem -> "Equatorial"] - {SiderealTime[GeoPosition[{0, 0}], date], 0}]]

Then evaluate:

paris = GeoPosition[Entity["City", {"Paris", "IleDeFrance", "France"}]];

sunsets = Table[
   Sunset[paris, DateObject[{2022, 1, n}, TimeZone -> 1]],
   {n, 1, 365}
];

tg = Table[
    GeoGraphics[{
        NightHemisphere[s],
        GeoStyling[Yellow], GeoDisk[subSolarPoint[s], Quantity[100, "Miles"]], 
        GeoStyling[Black],  GeoDisk[subLunarPoint[s], Quantity[100, "Miles"]]
    }, GeoRange -> "World", ImageSize -> Large],
    {s, sunsets}
];

ListAnimate[tg]