The Orbit and Perigee of the Flamsteed comet
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:
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.