Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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 Hapsburg 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 approachin 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:

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

enter image description here


I'd like to find the exact date and time of the perigee of the Flamsteed comet and to inset points of locations (at that time) of the first $4$ planets on the graphics.

share|improve this question
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…. My expectation is to determine it up to, say 6-hour precision or even one day. – Artes May 15 '13 at 13:10
As I noted, Mathematica does have some of the orbital elements missing, which complicates things. I'll see if there's a workaround in the meantime. (As a warning for other people: AstronomicalData["CometC1683O1", {"Position", {1683, 9, 12}}] returns Missing["Variable"]; similarly for other dates within Artes's period of interest.) – J. M. May 15 '13 at 13:24

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

Flamsteed comet perigee

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"],
        (* 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] :=
                     AstronomicalData["Earth", {"Position", DateList[s]}, TimeZone -> 0.]]

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

orbits of planets and a comet

share|improve this answer
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: PlanetaryAstronomy`Private`ObliquityLaskar[].) – J. M. May 18 '13 at 17:59
@cormullion: Your puns were so powerful they blew the accepted-answer checkmark away. – DumpsterDoofus Oct 5 '14 at 21:01

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.