Calculating upper culmination of the Moon

Question

I’m completely confused by the new astronomical data system in Mathematica. For every day in some interval, I want to calculate {azimuth, altitude} of the Moon in moment of its upper culmination.

I’ve already realized that I can’t get this data directly. But we can calculate the culmination as the exact midpoint between RiseTime and SetTime, and further we can use MoonPosition.

But how can I get RiseTime and SetTime for the Moon? From PlanetaryMoon data? And how can I get correct time of upper culmination, since it’s not just a half sum of RiseTime and SetTime?

UPD, PS
Thank you for great answers! Using it, I got a graph of the Moon’s positions at transition moment of a year.

I'm interesting is it show really Moon's librations, or just inaccuracy of data?

Answer 1

Mathematica provides "upper culmination" as the TransitTime property of PlanetaryMoon data. With the transit dates, we can get azimuth and altitude from MoonPosition.

The TransitTime property has two qualifiers for date and location. For example, let's use EntityValue and TransitTime qualifiers to find the Moon's transit time for October 2023 at Lowell Observatory, Flagstaff, Arizona, United States.

EntityValue[EntityProperty["PlanetaryMoon", "TransitTime"], "Qualifiers"]
(* {"Date", "Location"} *)

Let's get the rise, transit, and set times for location and dates in October 2023. First, use the rise times to find the following transit times, and then use the transit times to find the following set times. This gives us a sequence of rise time, transit time and set time for each date.

(*Lowell Observatory, Flagstaff, Arizona, United States*)
location = Entity["AstronomicalObservatory", "LowellFlagstaff"];
dates = DateRange[DateObject[{2023, 10, 1}], DateObject[{2023, 10, 30}]];

(*find next rise time for each date*)
r=EntityValue[Entity["PlanetaryMoon",  "Moon"],
 EntityProperty["PlanetaryMoon",
  "RiseTime", {"Date" -> #, "Location" -> location}]] &/@ dates;

(*find transit time following each rise time*)
t=EntityValue[Entity["PlanetaryMoon", "Moon"],
 EntityProperty["PlanetaryMoon",
  "TransitTime", {"Date" -> #, "Location" -> location}]] &/@ r;

(*find set time following each transit time*)
s=EntityValue[Entity["PlanetaryMoon",  "Moon"],
 EntityProperty["PlanetaryMoon",
  "SetTime", {"Date" -> #, "Location" -> location}]] &/@ t;

Demonstate that the transit times are not the midpoint of rise time and set times. The times are accurate to the nearest minute.

r + (s - r)/2 - t
(* {5 min, 4 min, 2 min, 0 min, 1 min, -2 min, -3 min, -4 min, -5 min, -5 min, -5 min, -5 min, -4 min, -4 min, -4 min, -4 min, -3 min, -3 min, -1 min, 0 min, 2 min, 3 min, 4 min, 5 min, 5 min, 5 min, 5 min, 5 min, 4 min, 3 min} *)

Use MoonPosition to get azimuths and altitudes for each transit time. The azimuths are close to 180 degrees but not equal because of the limited accuracy of the times. Use DateListPlot to graph the Moon's altitude at transit.

{azimuths,altitudes} = MoonPosition[location,t]["PathComponents"];
Through[{Mean, StandardDeviation}[azimuths]]
(* {180.04°, 0.195364°} *)

DateListPlot[altitudes,
  PlotLabel -> Style["Moon Altitude at Transit", Bold, Black]]

Display a table of rise, transit, and set times.

fmt = {"Day", " ", "MonthNameShort", " ", "Year", " ", "Time"};
TableForm[Transpose[{DateString[#, fmt[[;;-2]]] &/@ dates,
  DateString[#, fmt] &/@ r,
  DateString[#, fmt] &/@ t,
  DateString[#, fmt] &/@ s}],
TableHeadings -> {None, {"Date", "Next Rise", "Next Transit", "Next Set"}}]

Answer 2

Just make it a maximization problem.

For example, set some given date and location:

pos = GeoPosition@
  Entity["City", {"Chicago", "Illinois", "UnitedStates"}]
t0 = DateObject[{2023, 10, 1, 0, 0}, "Minute", "Gregorian", -4.`]

Now create an objective function that returns the true altitude at each time of the day from 0 to 24 hrs:

tUnit = Quantity[1, "hour"]
obj[i_?NumericQ] := 
 QuantityMagnitude@
  Last@MoonPosition[pos, t0 + i tUnit, 
    AltitudeMethod -> "TrueAltitude"]

and maximize the altitude:

maxTime = NArgMax[{obj[x], 0 < x < 24}, x]
maxPos = 
 MoonPosition[pos, t0 + maxTime* tUnit, 
  AltitudeMethod -> "TrueAltitude"]

(*3.32927*)
(* {Quantity[181.565, "AngularDegrees"], 
 Quantity[61.1416, "AngularDegrees"]}*)

The same can be done for lower culmination as well:

minTime = NArgMin[{obj[x], 0 < x < 24}, x]
minPos = 
 MoonPosition[pos, t0 + minTime* tUnit, 
  AltitudeMethod -> "TrueAltitude"]
(*15.6206*)
(*{Quantity[358.755, "AngularDegrees"], 
 Quantity[-32.5243, "AngularDegrees"]}*)
```

Answer 3

You can get the RiseTime for a specific location and date:

PlanetaryMoonData["Moon", 
  EntityProperty["PlanetaryMoon", "RiseTime", 
  {"Location" -> GeoPosition[{30, -95}], 
   "Date" -> DateObject[{2023, 12, 1}]}] ]

For example, here are RiseTime and SetTime for five months:

risetimes = 
  PlanetaryMoonData["Moon", EntityProperty["PlanetaryMoon", "RiseTime", 
  {"Location" -> GeoPosition[{30, -95}], 
   "Date" -> #}] & /@ 
    DateRange[DateObject[{2023, 11, 1}], DateObject[{2024, 3, 31}]]];

settimes = 
  PlanetaryMoonData["Moon", EntityProperty["PlanetaryMoon", "SetTime", 
  {"Location" -> GeoPosition[{30, -95}], 
   "Date" -> #}] & /@ 
    DateRange[DateObject[{2023, 11, 1}], DateObject[{2024, 3, 31}]]];

To convert them to minutes of the day

risetimemins = #[[4]] 60 + #[[5]] & /@ DateList /@ risetimes
settimemins = #[[4]] 60 + #[[5]] & /@ DateList /@ settimes

ListPlot[{risetimemins, settimemins}, Joined -> True]

Updated

Azimuth:

azimuths = 
  PlanetaryMoonData["Moon", 
   EntityProperty["PlanetaryMoon", "Azimuth", 
   {"Location" -> GeoPosition[{30, -95}], 
    "Date" -> #}] & /@ 
    DateRange[DateObject[{2023, 11, 1}], DateObject[{2024, 3, 31}]]];

This is gives results that look like

azimuths[[1 ;; 3]]
  (* {Quantity[MixedMagnitude[{67, 10, 42.}], 
  MixedUnit[{"AngularDegrees", "Arcminutes", "Arcseconds"}]], 
 Quantity[MixedMagnitude[{60, 30, 57.}], 
  MixedUnit[{"AngularDegrees", "Arcminutes", "Arcseconds"}]], 
 Quantity[MixedMagnitude[{54, 22, 28.}], 
  MixedUnit[{"AngularDegrees", "Arcminutes", "Arcseconds"}]]} *)

ListPlot[Normal@azimuths, Joined -> True]

Transit altitude (apologies - I'm not sure if this is the same as the upper culmination):

alt = PlanetaryMoonData["Moon", 
   EntityProperty["PlanetaryMoon", 
      "AltitudeTransit", {"Location" -> GeoPosition[{30, -95}], 
       "Date" -> #}] & /@ 
    DateRange[DateObject[{2023, 11, 1}], DateObject[{2024, 3, 31}]]];

(* example results *)
alt[[1 ;; 3]]
  (* {Quantity[MixedMagnitude[{87, 43.55}], 
  MixedUnit[{"AngularDegrees", "Arcminutes"}]], 
 Quantity[MixedMagnitude[{88, 14.54}], 
  MixedUnit[{"AngularDegrees", "Arcminutes"}]], 
 Quantity[MixedMagnitude[{87, 8.70}], 
  MixedUnit[{"AngularDegrees", "Arcminutes"}]]} *)

ListPlot[Normal@alt, Joined -> True]

Answer 4

In version 14.1, there is a new function AstroRiseSet, which can also calculate culminations:

AstroRiseSet[Entity["PlanetaryMoon", "Moon"], "UpperCulmination"]
(* DateObject[{2024, 8, 1, 3, 33, 7.35101}, "Instant", "Gregorian", -5.]  *)

Plot time of upper culminations throughout this year:

upperCulminations = AstroRiseSet[Entity["PlanetaryMoon", "Moon"], "UpperCulmination", 
  Here, DateInterval[CurrentDate["Year"]]];

DateListPlot[Comap[{DateObject, TimeObject}] /@ upperCulminations, 
 ScalingFunctions -> "Date"]