Geocentric, Earth-fixed, coordinates of Sun, Moon (and others objects eventually)

Question

I know I can get "HelioCoordinates" of the Earth and the Moon with

PlanetData["Earth", "HelioCoordinates"]
PlanetaryMoonData["Moon", "HelioCoordinates"]

{Quantity[-0.897541, "AstronomicalUnit"], Quantity[-0.449087, "AstronomicalUnit"], Quantity[-1.19589*10^-6, "AstronomicalUnit"]}

{Quantity[-0.895127, "AstronomicalUnit"], Quantity[-0.44907, "AstronomicalUnit"], Quantity[0.000035853, "AstronomicalUnit"]}

What I eventually want is to convert those coordinates to a geocentric, earth-fixed, cartesian reference frame, essentially the same used within GeoPositionXYZ.

Now the translation between the two reference frame is obvious, i.e.

PlanetData["Earth", "HelioCoordinates"]

I admit I'm weak in astronomy, and so I searched but I haven't found so far a fast and easy way to get the angles between two reference frames and/or do the transformation.

What I tried to do is to follow this way and to get the longitude $\lambda$ of the observer for which the local hour angle LHA is zero, from right ascension and Greenwich sidereal time GST, using

PlanetaryMoonData["Moon", "RightAscension"] - SiderealTime[
 GeoPosition@Entity["AstronomicalObservatory", "Greenwich"]]

but I'm not sure of these ideas and I still need two other angles.

I suppose with the help of

PlanetData["Earth", "HelioCoordinates"]
PlanetaryMoonData["Moon", "HelioCoordinates"]
MoonPosition[CelestialSystem -> "Equatorial"]
SunPosition[CelestialSystem -> "Equatorial"]
SiderealTime[...]

and other similar function should be possible, if not easy.

Trying to get some insight I read many examples on the Mathematica documentation and I also encountered a problem reproducing some example. In the end, my goal is to get the geocentric, earth-fixed, cartesian coordinates at a given date and time, but apparently when I add a "Date" qualifier for relevant Entity properties I get a "Missing[NotAvailable]" response. For example with

PlanetaryMoonData["Moon", 
   EntityProperty["PlanetaryMoon", 
    "DistanceFromEarth", {"Date" -> #}]] & /@ {DateObject[{2014, 4, 9,
     22, 0, 0}]}

Answer 1

You have most of the pieces here already.

UPDATED

The proper coordinates

In this problem, we want to get the positions of astronomical objects in terms of a Cartesian system that is geocentric and rotates with the Earth. In astronomy, the positions of objects are commonly given in terms of right ascension (RA) and declination (Dec), which are similar to longitude and latitude projected onto the sky, but do not rotate with the Earth. Because of this similarity RA and Dec are natural coordinates to start with: it is simple to take the RA and Dec and then boost them into a frame that rotates with the Earth. One can then make a simple transformation from spherical into Cartesian coordinates. $$ (r, Dec, RA) \rightarrow (r, lat, long) \rightarrow (x, y, x) $$

Aside: Because the Earth's axis moves slowly with respect the distant stars (precession and nutation etc.), there is some ambiguity as to whether RA and Dec should move with the Earth or remain fixed. As detailed at length in this question, the various astronomical Data functions when applied to the Sun, the Moon, the planets, and the other planetary moons supply precessing (true geocentric) right ascension and declination. But, for stars besides the Sun, StarData provides a non-precessing (fixed) RA and Dec. All appear to be accurate up to about a degree or RA and less than a degree of Dec. I haven't checked yet what the other astronomical Data functions do, so I will restrict this answer to the Sun, the Moon, and planets. It should be sufficiently easy to generalize.

Conversion into Earth-rotating coordinates

To convert from geocentric, precessing RA into a coordinate in an Earth-rotating reference frame, we calculate the local hour angle as $$ LHA(long,t) = LST(long,t) - RA(t) $$ Here, $LST$ is the local sidereal time. Since we are interested in the position of the Sun or Moon with respect to Earth longitude and latitude, we pick our earth position to be $(lat, long) = (0, 0)$. Conveniently, the local sidereal time, on the prime meridian is Greenwich sidereal time (GST, closely related to GMT), so a "longitude-type" angle (only "longitude-type", because we don't need to require this angle to only run between -180 degrees and 180 degrees) for the moon or sun is given as $$ LHA(0,t) = 360^\circ - \left( GST(t) - RA(t) \right) $$ (the 360 degrees converts the angle from westward-directed to eastward-directed, so that we can get a right-handed coordinate system, with the x-axis sticking out from the equator-prime meridian point, the y-axis going off through the equator at 90 degrees, and the z-axis along up the axis of rotation)

Mathematica implementation

We can put this into Mathematica now fairly easily:

Sun position

SunPositionXYZ[date_] := With[
  {date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},

  Module[{pos, GST, r, lat, long},
   pos = StarData["Sun",
     {EntityProperty["Star", "RightAscension", {"Date" -> date0}], 
      EntityProperty["Star", "Declination", {"Date" -> date0}]}];
   pos = UnitConvert[pos, "Degrees"];

   GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];

   r = PlanetData["Earth", 
     EntityProperty["Planet", "DistanceFromSun", {"Date" -> date0}]];

   long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
   lat = pos[[2]];

   <|
    x -> r Cos[lat] Cos[long],
    y -> r Cos[lat] Sin[long],
    z -> r Sin[lat]
    |>
 ]
]

Moon position

MoonPositionXYZ[date_] := With[
  {date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},

  Module[{pos, GST, r, lat, long},
   pos = PlanetaryMoonData["Moon",
     {EntityProperty["PlanetaryMoon", "RightAscension", {"Date" -> date0}], 
      EntityProperty["PlanetaryMoon", "Declination", {"Date" -> date0}]}];
   pos = UnitConvert[pos, "Degrees"];

   GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];

   r = PlanetaryMoonData["Moon", 
     EntityProperty["PlanetaryMoon", "DistanceFromSun", {"Date" -> date0}]];

   long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
   lat = pos[[2]];

   <|
    x -> r Cos[lat] Cos[long],
    y -> r Cos[lat] Sin[long],
    z -> r Sin[lat]
    |>
 ]
]

Planet positions

PlanetPositionXYZ[planet_, date_] := With[
  {date0 = TimeZoneConvert[date, 0], pos0 = GeoPosition[{0, 0}]},

  Module[{pos, GST, r, lat, long},
   pos = PlanetData[planet,
     {EntityProperty["Planet", "RightAscension", {"Date" -> date0}], 
      EntityProperty["Planet", "Declination", {"Date" -> date0}]}];
   pos = UnitConvert[pos, "Degrees"];

   GST = UnitConvert[SiderealTime[pos0, date0], "Degrees"];

   r = PlanetData[planet, 
     EntityProperty["Planet", "DistanceFromEarth", {"Date" -> date0}]];

   long = Quantity[360, "Degrees"] - (GST - pos[[1]]);
   lat = pos[[2]];

   <|
    x -> r Cos[lat] Cos[long],
    y -> r Cos[lat] Sin[long],
    z -> r Sin[lat]
    |>
 ]
]

For moons, just make the substitution PlanetData -> PlanetaryMoonData and "Planet" -> "PlanetaryMoon".

Date and Time

Run these using DateObject to specify the time and date:

MoonPositionXYZ[DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> -3]]
SunPositionXYZ[DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> 10]]
PlanetPositionXYZ["Venus", 
 DateObject[{2015, 4, 17, 10, 0, 0}, TimeZone -> 5]]

(* -> <|x -> 0.857983 au, y -> -0.512016 au, z -> 0.0685721 au|>

      <|x -> -0.9875 au, y -> -0.00628945 au, z -> 0.179209 au|>

      <|x -> -0.824832 au, y -> 0.583867 au, z -> 0.431657 au|> *)

Mathematica should convert all of the DateObjects into some kind of universal date (probably Julian date) internally when they are passed to the various astronomical position functions, so as long as the date and time are uniquely specified, it shouldn't matter what timezone is used (note that Mathematica may use your timezone as default). That said, it is probably safest to specify the time and date fully (as above), rather than just using DateObject[TimeZone -> x], as this seems to give an unexpected date sometimes (see this question).

Other objects

For objects such as stars for which Mathematica gives non-precessing RA and Dec, we must first convert it to a precessing RA. Mathematica already does this, but the functionality is buried in the innards of the astronomical Data functions. I have not yet come up with a good way of doing the conversion.