The Orbit and Perigee of the Flamsteed comet

Question

Historical context

This year we have the 330-th anniversary of the Battle of Vienna - one of the great formative events of European history, it took place on September 12, 1683.

Kara Mustafa, Grand Vizier of the Ottoman Empire, had laid siege to the Habsburg capital and was on the verge of capturing it when a relieving Christian army under the overall command of Jan III Sobieski, King of Poland, swept into the Turkish ranks.

During the siege of Vienna by the Islamic power, before Sobieski's forces joined (on September 11) the rest of the Holy League, there had appeared a comet (later called Flamsteed) on the sky at the end of July and could be seen until September.

Newton's Principia Mathematica on the comet

In the third book of Philosophiae Naturalis Principia Mathematica Isaac Newton says:

The comet of 1683 (also according to the observations of Hevelius) at the end of July, when it was first sighted, was moving very slowly, advancing about $40'$ or $45'$ in its orbit each day. From that time its daily motion kept increasing continually until 4 September when it came to about $5^{\circ}$. Therefore in all this time the comet was approaching the earth. This is gathered also from the diameter of the head, as measured with micrometer, since Hevelius found it to be on 6 August only $6'5''$ including the coma, but on 2 September $9'7''$. Therefore the head appeared far smaller at the begining than at the end of the motion, as Hevelius also reports. Accordingly in all this time, because of receding from the sun it decreased with respect to its light, notwithstanding its approach to the earth.

Astronomical Data

With help of built-in AstronomicalData we can easily draw the orbits of the comet and the first 4 planets:

Graphics3D[
  {{#1, AstronomicalData[#2, "OrbitPath"]} & @@@ 
    Transpose[{ {Orange, Green, Blue, Red}, Take[ AstronomicalData["Planet"], 4]} ], 
   {Magenta, Line[ AstronomicalData[ 
                     AstronomicalData["CometC1683O1"], "OrbitPath"][[1, 28 ;; 195]]]}
  }, Boxed -> False]

Problem How can we find the exact date and time of the perigee of the Flamsteed comet and to inset points of locations (at that time) of the terrestial planets on the graphics?

Edit

To broaden the historical context it would be reasonable to mention that Scutum was one of few constellations separated in modern times in 1684 by Johannes Hevelius (whose research was supported by Sobieski) in commemoration of the victory of Christian forces led by Sobieski in the Battle of Vienna.

Originally it was named Scutum Sobiescianum (Sobieski's shield) and later abbreviated to Scutum.

ConstellationData[ Entity["Constellation", "Scutum"], 
                   EntityProperty["Constellation", "ConstellationGraphic"]]

It would be interesing to demonstrate the trajectory of the comet on the sky from the geocentric reference frame (Earth-centered inertial) against given constellations in August 1683 before the battle. Can we go further with new Mathematica functionality like PlanetData and CometData with respect to the former capabilities of AstronomicalData?

Answer 1

It took me quite a while, but finally, here's a visualization of the perigee of Flamsteed's comet:

---

I should first note two things: first, some of the needed data for computing the orbit of comet C/1683 O1 was missing in AstronomicalData["CometC1683O1", "Properties"], and I had to pull information from external sources to supplement the information available; second, after using the combined data, the orbit path I obtained didn't quite match the one from AstronomicalData["CometC1683O1", "OrbitPath"], and since I couldn't seem to access the appropriate ephemerides for a proper comparison, I'm not sure about the correctness. Nevertheless, what I have might be of some use.

As always, most of the formulae are adapted from Jean Meeus's Astronomical Algorithms (and the related book Astronomical Formulæ for Calculators, also by Meeus); pointers to formulae not in Meeus's work will be indicated in comments.

First, a few auxiliary routines. Here's a routine for the Julian Day Number (the same routine in this answer):

Options[jd] = {"Calendar" -> "Gregorian"};

jd[{yr_Integer, mo_Integer, da_?NumericQ, rest__}, opts : OptionsPattern[]] := 
  Module[{y = yr, m = mo, h}, If[m < 3, y--; m += 12];
       h = Switch[OptionValue["Calendar"],
                  "Gregorian", (Quotient[#, 4] - # + 2) &[Quotient[y, 100]],
                  "Julian", 0,
                  _, Return[$Failed]];
       Floor[365.25 y] + Floor[30.6001 (m + 1)] + da + FromDMS[{rest}]/24 + 1720994.5 + h
   ]

jd[{yr_Integer, mo_Integer, da_?NumericQ}, opts : OptionsPattern[]] :=
   jd[{yr, mo, da, 0, 0, 0}, opts]

jd[opts : OptionsPattern[]] := jd[DateList[], opts]

Here is a routine for the mean obliquity of the ecliptic $\varepsilon$. Since the period of interest is a rather long time ago, I decided to use the formula in Laskar's article that has a larger domain of validity, instead of the conventional formula used by the USNO (which was used in this answer):

MeanEclipticObliquity[args___] := Module[{T}, T = (jd[args] - 2451545)/3652500;
  (84381.406 + T (-4680.93 + T (-1.55 + T (1999.25 + T (-51.38 + T (-249.67 + T (-39.05 +
   T (7.12 + T (27.87 + T (5.79 + 2.45 T))))))))))/3600]

Here, now, is the main routine for reckoning the position (in heliocentric rectangular equatorial coordinates) of Flamsteed's comet from its orbital elements. The formulae for bodies with parabolic orbits was taken from chapter 33 of Astronomical Algorithms; the perihelion distance (one of the orbital elements missing in AstronomicalData) of C/1683 O1 was taken from here, with the data attributed to Halley.

flamsteedCometPosition[date_] := 
   Block[{(* astronomical unit *) AU = 1.495978707*^11,
          (* Gaussian gravitational constant *) k = 0.01720209895,
          a, b, c, q, r, s, v, W, α, β, γ, ε, ι, ω, Ω},
         Ω = AstronomicalData["CometC1683O1", "AscendingNodeLongitude"] °;
         ι = AstronomicalData["CometC1683O1", "Inclination"] °;
         ω = AstronomicalData["CometC1683O1", "PeriapsisArgument"] °;
         ε = MeanEclipticObliquity[date] °;
         {{a, α}, {b, β}, {c, γ}} =
         MapThread[{Norm[{##}], ArcTan[##]} &,
                   {{-Sin[Ω] Cos[ι], Cos[Ω] Cos[ι] Cos[ε] - Sin[ι] Sin[ε],
                     Cos[Ω] Cos[ι] Sin[ε] + Sin[ι] Cos[ε]},
                    {Cos[Ω], Sin[Ω] Cos[ε], Sin[Ω] Sin[ε]}}];
         (* perihelion distance of C/1683 O1 *) q = 0.5602;
         W = (3 k/Sqrt[2]) q^(-3/2)
             DateDifference[AstronomicalData["CometC1683O1", "PerihelionTime", "Epoch"],
                            date];
        (* solution of Barker's equation *) s = Root[#^3 + 3 # - W &, 1];
        (* radius vector *) r = q (1 + s^2);
        (* true anomaly *) v = 2 ArcTan[s];
        r {a, b, c} Sin[{α, β, γ} + ω + v] AU]

To reckon the date of the comet's perigee, we can now do this (note the explicit setting of the TimeZone option so that the reckoning is done in Greenwich time):

dist[s_?NumericQ] :=
   EuclideanDistance[flamsteedCometPosition[DateList[s]],
                     AstronomicalData["Earth", {"Position", DateList[s]}, TimeZone -> 0.]]

perigee =
  DateList[First[FindArgMin[dist[s],
                            {s, AbsoluteTime[{1683, 7, 1}], AbsoluteTime[{1683, 9, 30}]}]]]
   {1683, 9, 3, 3, 47, 13.4369}

Finally, here's how to generate the picture at the beginning of this answer:

With[{AU = 1.495978707*^11}, 
     Graphics3D[{{Yellow, AbsolutePointSize[30], Point[{0, 0, 0}]},
                 {LightYellow,
                  {AbsolutePointSize[4], Point[flamsteedCometPosition[perigee]/AU]},
                  {Directive[AbsoluteDashing[{5, 5}], AbsoluteThickness[1]], 
                   Line[Table[flamsteedCometPosition[DatePlus[perigee, k]]/AU,
                              {k, -30, 0}]]}},
                 {AbsoluteThickness[2], MapThread[
                  Function[{planet, color, size},
                           {{color, AbsolutePointSize[size],
                             Point[AstronomicalData[planet, {"Position", perigee},
                                                    TimeZone -> 0.]/AU]},
                            {Lighter[color, 1/5], AstronomicalData[planet, "OrbitPath"]}}],
                  {Take[AstronomicalData["Planet"], 4],
                   {Gray, Orange, Blue, Red}, {6, 12, 12, 8}}]}},
                Background -> Black, Boxed -> False, ViewPoint -> {1.3, -2.4, 1.5}]]

---

As a bonus, here's an animation of the orbits of the terrestrial planets and Flamsteed's comet, from August 1 to September 15, 1683:

Answer 2

I've decided to write the "modernization" of my old code as a separate answer, to keep the previous answer (relatively) uncluttered. Overall, I've found the need for gymnastics related to QuantityMagnitude[] and the various Entity[] functions to be quite annoying. Maybe somebody else can make the following code a bit nicer:

(* force loading of internal PlanetaryAstronomy` context *)
AstronomicalData["Earth", "Position"];

flamsteedCometPosition[date_DateObject] := 
    Block[{k, a, b, c, q, r, s, v, W, α, β, γ, ε, ι, ω, Ω},

          k = QuantityMagnitude[UnitConvert[Quantity[1, "GaussianGravitationalConstant"],
                                ("AstronomicalUnit")^(3/2)/("Days" Sqrt["SolarMass"])]];

          {Ω, ι, ω} = QuantityMagnitude[UnitConvert[CometData["CometC1683O1",
          {"AscendingNodeLongitude", "Inclination", "PeriapsisArgument"}],
          "Radians"]];

          ε = PlanetaryAstronomy`Private`ObliquityLaskar[(JulianDate[date] -
                                                          2451545)/3652500];

  {{a, α}, {b, β}, {c, γ}} = MapThread[ToPolarCoordinates[{##}] &,
  {{-Sin[Ω] Cos[ι], Cos[Ω] Cos[ι] Cos[ε] - Sin[ι] Sin[ε],
    Cos[Ω] Cos[ι] Sin[ε] + Sin[ι] Cos[ε]},
   {Cos[Ω], Sin[Ω] Cos[ε], Sin[Ω] Sin[ε]}}];

  (* perihelion distance of C/1683 O1; still missing after all these years *)
  q = 0.5602;

  W = (3 k/Sqrt[2]) q^(-3/2) QuantityMagnitude[
       DateDifference[CometData["CometC1683O1", "PeriapsisTimeLast"], 
                      date]];
  (* solution of Barker's equation *) s = Root[#^3 + 3 # - W &, 1];
  (* radius vector *) r = q (1 + s^2);
  (* true anomaly *) v = 2 ArcTan[s];
  Quantity[r {a, b, c} Sin[{α, β, γ} + ω + v], "AstronomicalUnit"]]

dist[s_?NumericQ] := EuclideanDistance[
flamsteedCometPosition[DateObject[s, TimeZone -> 0]], 
PlanetData["Earth", EntityProperty["Planet", 
"HelioCoordinates", {"Date" -> DateObject[s, TimeZone -> 0]}]]]

For some reason, however, I am now getting a different date for the perigee:

perigee = DateObject[First[FindArgMin[QuantityMagnitude[dist[s]],
                                      {s, AbsoluteTime[{1683, 7, 1}], 
                                          AbsoluteTime[{1683, 9, 30}]}]],
                     TimeZone -> 0]
   DateObject[{1683, 9, 2}, TimeObject[{8, 43, 54.864279}, TimeZone -> 0.]]

About an offset of a day; huh.

Anyway, the associated image can now be done like so:

Graphics3D[{{Yellow, AbsolutePointSize[30], Point[{0, 0, 0}]},
            {LightYellow, {AbsolutePointSize[4], 
             Point[QuantityMagnitude[flamsteedCometPosition[perigee]]]},
            {Directive[AbsoluteDashing[{5, 5}], AbsoluteThickness[1]], 
             Line[Table[QuantityMagnitude[flamsteedCometPosition[
                        DatePlus[perigee, k]]], {k, -30, 0}]]}},
            {AbsoluteThickness[2],
             MapThread[Function[{planet, size},
             {{FromEntity[PlanetData[planet, "Color"]], AbsolutePointSize[size], 
              Point[QuantityMagnitude[PlanetData[planet,
                    EntityProperty["Planet", "HelioCoordinates",
                                   {"Date" -> perigee}]]]]},
             {Lighter[FromEntity[PlanetData[planet, "Color"]], 1/10], 
              PlanetData[planet, "OrbitPath"]}}],
             {PlanetData[EntityClass["Planet", "InnerPlanet"]],
              {6, 12, 12, 8}}]}},
           Background -> Black, Boxed -> False, ViewPoint -> {1.3, -2.4, 1.5}]

As mentioned previously, adding the starry background looks to be a messy affair; I couldn't get ConstellationData[] to do what I wanted, so I omitted it for now. I'll edit this soon as I figure that one out.