# The Orbit and Perigee of the Flamsteed comet

### 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?

-
Only slightly related: would you happen to have any ephemeris on hand that might have formulae for this comet's path? – J. M. May 15 '13 at 12:00
(In case it wasn't apparent why I was asking, see the output of AstronomicalData["CometC1683O1", "OrbitRules"].) – J. M. May 15 '13 at 12:32
@J.M. Let's say, no further information right now, besides that available here reference.wolfram.com/mathematica/note/…. My expectation is to determine it up to, say 6-hour precision or even one day. – Artes May 15 '13 at 13:10
– Dr. belisarius May 15 '13 at 13:30
If push comes to shove, and this can't be done with Mathematica, you might be interested in this. – J. M. May 15 '13 at 13:30

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, γ}} =
{{-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}]]}},
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:

-
Thanks for a good answer. I'll wait a bit with acceptance to encourage others with different approaches. – Artes May 18 '13 at 14:54
Sure @Artes, take all the time you need. I was also planning to post an approach based on integrating the underlying DEs, but that is proving to be more difficult to do... – J. M. May 18 '13 at 14:57
you Mustafa upvote from me for this stellar answer – cormullion May 18 '13 at 15:36
(As it turns out, Laskar's formula is built-in, but undocumented: PlanetaryAstronomyPrivateObliquityLaskar[].) – J. M. May 18 '13 at 17:59
@J.M. this is nice. Given the computational tools available, could you complicate it by adding green lines at perigees labeled w/ times and distances? – alancalvitti Mar 14 at 17:40

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[{(* Gaussian gravitational constant *) k = 0.01720209895, a, b, c, q,
r, s, v, W, α, β, γ, ε, ι, ω, Ω},

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

ε = PlanetaryAstronomyPrivateObliquityLaskar[(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],
{{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.

-
Thanks for your support. Why community wiki, can it be switched to the normal mode? Is it eligible to the bounty? Perhaps it hasn't been quite clear, but my expectation is to draw the trajectory from the geocentric point of view. I haven't checked yet but I belive it is not too difficult. – Artes Mar 22 at 6:14
Discrepancy between the time of "new" and "old" perigee is not enormous, I believe it is due to almost parabolic tracectory of the comet. I expected from the begining such an indeterminacy would be acceptable. I cannot test and compare the both solutions thoroughly at the moment since I have only the 10.3 version at hand. As far as the current answer is formulated I'm going to accept and reward the bounty to the old answer. However if additional improvements are realized I' ll switch the accept mark to this one, so that I hope you expectedly earn the populist badge. – Artes Mar 22 at 7:02
Hello, sorry for the late reply. I had set it CW since it is a mere translation of work I have previously done (even tho annoying contortions were involved). As you say making a geocentric view is indeed not too difficult, but it is definitely tedious. If you don't mind waiting, I'll edit this answer in a few to include that. (ConstellationData[] remains intransigent, tho. I may have to use StarData[]` manually to do what I want.) – J. M. Mar 22 at 22:06
BTW, it is not "almost"; the equations I used definitely make the parabolic assumption. ;) – J. M. Mar 22 at 22:06