(* Content-type: application/mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 6.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 145, 7] NotebookDataLength[ 53883, 1761] NotebookOptionsPosition[ 49141, 1615] NotebookOutlinePosition[ 49533, 1632] CellTagsIndexPosition[ 49490, 1629] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["Kepler's Laws", "Title"], Cell["\<\ AP1603 M. E. Mauel\ \>", "Subsubtitle"], Cell[CellGroupData[{ Cell["Kepler", "Section"], Cell[TextData[{ "The German astronomer Johannes Kepler (b. Dec. 27, 1571, d. Nov. 15, \ 1630) was the first strong supporter of Copernicus's view that the planets \ revolved the Sun and the discoverer of the famous three laws of planetary \ motion. Kepler had planned to enter religious life, attended seminaries at \ Adelberg and Maulbronn, and studied theology, philosophy, and mathematics at \ the University of Tubingen. At Tubingen, Kepler's scientific ability \ attracted notice. In his early 20's, Kepler became a professor of mathematics \ and astronomy at Graz.\nAt the age of 24, Kepler published ", StyleBox["Mysterium cosmographicum", FontSlant->"Italic"], " (Cosmographic Mystery, 1596), in which he defended the Copernican theory \ and described his ideas on the structure of the planetary system. Influenced \ by the Pythagoreans, Kepler viewed the universe as being governed by \ geometric relationships that conform to the inscribed and circumscribed \ circles of the five regular polygons.\nAlthough he was not a Copernican \ himself, Tycho Brahe, the mathematician at the court of Emperor Rudolph II at \ Prague, was so impressed with Kepler's work that in 1600 he invited Kepler to \ come to Prague as his assistant. Confronted with the Catholic persecution of \ the Protestant minority in Graz, Kepler gladly accepted. When Brahe died the \ following year, Kepler was appointed his successor and thus inherited Brahe's \ scientific legacy.\nThis legacy included many accurate positional \ determinations of the planets, especially those of Mars. Kepler now embarked \ on an intensive study of the true orbits of the planets. Abandoning the \ ancient belief that the planets must move in perfect circles, Kepler \ concentrated on Mars. He proved that the orbit of Mars is an ellipse, with \ the Sun occupying one of its two foci. This, the first of Kepler's laws of \ planetary motion, appeared in ", StyleBox["Astronomia nova", FontSlant->"Italic"], " (New Astronomy) in 1609, with the second \"law of areas\" governing \ planetary velocity.\nAlways guided by the concept of beauty in the structure \ of the universe, and specifically by a theory of harmony in geometric \ figures, numbers, and music, Kepler, in his ", StyleBox["Harmonices mundi", FontSlant->"Italic"], " (Harmonies of the World, 1619), announced his third law--a relationship \ between the orbital periods and the distances of the planets from the Sun. \ His belief that the Sun regulates the velocity of the planets was a milestone \ in scientific thought, laying the foundation for Newton's theory of universal \ gravitation.\nAmong Kepler's numerous scientific contributions are an \ influential treatise on the theory of optics (1604), a treatise on optics as \ applied to telescope lenses (1611), a work offering physical explanations of \ the appearance of a nova in 1604, and an enthusiastic acceptance of and \ elaboration on Galileo's observations with a telescope (1610). His Epitome \ astronomiae ", StyleBox["Copernicanae", FontSlant->"Italic"], " (Introduction to Copernican Astronomy, 1618-21) became one of the most \ widely read treatises on astronomy in Europe. Kepler's last great work, \ known as the ", StyleBox["Rudolphine", FontSlant->"Italic"], " ", StyleBox["Tables", FontSlant->"Italic"], " (1627), was a widely used compilation of accurate tables of planetary \ motion.\nThe posthumous ", StyleBox["Somnium", FontSlant->"Italic"], " (Dream, 1634), on which Kepler labored until shortly before his death, is \ indicative of his fertile mind. In this work, Kepler describes a journey to \ the Moon and discusses the existence of lunar inhabitants. \n(From Steven J. \ Dick, ", StyleBox["Groiler's Encyclopedia.", FontSlant->"Italic"], ")" }], "Text"], Cell[CellGroupData[{ Cell["Kepler's Three Laws", "Subsection"], Cell[TextData[{ "Kepler's three laws of planetary motion describe the shape of planetary \ motion on the plane, the velocities of the planets, and the distances of the \ planets from the Sun. The laws are:\n1. All planets move in elliptic orbits \ with the Sun at one focus.\n2. The line joining a planet to the Sun sweeps \ out equal areas in equal times.\n3. If ", StyleBox["T", FontSlant->"Italic"], " is equal to the orbital period of a planet (", StyleBox["i.e.", FontSlant->"Italic"], " the time required to complete one orbit), and if ", StyleBox["a", FontSlant->"Italic"], " is the semimajor axis of the orbit (", StyleBox["i.e.", FontSlant->"Italic"], " one half the long distance across the ellispe), then the ratio ", Cell[BoxData[ FormBox[ SuperscriptBox["T", "2"], TraditionalForm]]], "/", Cell[BoxData[ FormBox[ SuperscriptBox["a", "3"], TraditionalForm]]], " is a constant." }], "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Planetary Motion about a Massive Star", "Section"], Cell[GraphicsData["CompressedBitmap", "\<\ eJztW3uPFEUQH25vb3fvgXc8RTh273g/TmCPQwQBEQRORCWXUzGYeJ4QHkF5 CYoPfIavwP9gDPEDmPAZ+Ap8AUn4CutUdXVNz+yvZ3vuzsREJ2Futqq6qrq7 urq6qjk1c/PCuSszNy/OzjSOX5+5euHi7I3GsS+ux6DSoiiKnsX//mxE9N2K P+V1oBU/UYl+tOSJ/uB39D/qv4raSO8+At3xEw3Tu0agPv265idfQe8KgYb1 65KffIjeZQJt069zfvIBencTaEK/ZvzkVeV5QL/O+snL9GbQcf064yfvUqJT +vU+Jqc12EXvKaWccii7VBJ/nTaofdrFs9rqHYN6RQEloNjHij1pyCeUdQ8g /0Sxk4Z8NyBCDT8FsGOGRVN51gDRrGKFfBcgQg3PAdgRw2KH8uwHROcVK+Rj ClgMyC8o9rAh3w6IUMNLAHbQsNiqPIcA0WXFCvkWQIQaXgGw/YbFJuW5DBB9 rlgh3wiIUMOrALbXsNgAUMsB7JrKlobrAdEKALsBYBOGxTqAWglgNwHMWHw0 omq9BIi+BLBx07ABUIjFLQBrGhZ1gFoNYLcBzCyfaC1ArQGwrwBsp2ExDFAI 9jWA7TAskESk2R0Ae3nBWKChqwPYNwA2ZligCUQT/S2AbS/G4jsA22ZYrAKo EQD7HsC2GhYvAtQogN0FsC0LxgKtRbRmfwCwzYYF8gjIc/wIYJsWjAVyaMjx /QRgxsdCt4rc789eUQvBYqm3j2nYL94xC2ZRFdj9+/cZTgERhzYxhmD0PHny pDU4OMh4+vvo0aMWx03t7BBs1q/lEoDanKMlwUQrCbBqqjk9IyMjrdu32RO3 nj59yr+jAqp+BmBmKcAtHoUCVQDrzsBIJaN+mQfXoZUJ+EdUWyzKyPxZcd0C zmhRUIlBgNoKYHxYiSo8PTR1R48edafV6X+QRijwHC2mUVU0ev78ORsPGQ2Z +fT0dOvevXukaY4BlbNaMj2K/nNUfQGgtnlVLbWazSarFqsMDwBGtVqbatQ1 Gnh6pKsmpIb7WGG1Otl8t0ctMkjXwcRrYqHVSnuNkus1YPAnotBRAh05/G60 pJbtE9XoIKqHQL8Wa8hOdIBAaLPLacg71RCB0Haf05C3rqW+hnXTcMDXcDmB 7gKsCU562HIfP37MQOu/YjPpmitfEzcZ90srghY7TVT8N8phyvvWynym3WrP 5FHifzn82HmuIhAKE21wZzfglhrwPJjRKNJodmDGnmk1gVAMbJmRE6GhoxU7 Z2YrZS7Io1E3SUMZQjsXKHWwLpcphxZ8ZOz1s2C9hgmEzgrMgv17ZV4sqsnm Ff1F34vmyo53ek4M5oxLwmLawZ4/z6kWGFMPqZblkOH+QLBkljRnikWdZy9Y DWH7oWBlz036gthmYCYn1cH2rACyL9pBfQLWFmFGK4CYxQZblFmyXs8IloIh 0i6erKLMEs/0Ef820RN5IBrLWDt0Ns3hl7jP94Qf9VL2Nuq1iTuMBxV/wg8N REFZfDDhTWdKsDb2arVshrPMgyNWzA8NUgNwq/sFcVTM2+JOYU1+Rk4JtmM0 aHVpTGiSazRoiUeJB8PuRfTEfipHD7RZsh59iR4V5md9Hz02GI2FSyfLOqH0 EI4CTXRayhHMYRLn/puCdTbWVNQqgmTayrpBysjTNzxQdopmMiHhBKCzapBI mgBXDSfA4n+0X89TDR0Xn42j2DpNGn9bBWncLDl90wQugIIW1r4wzBTacbL5 DtKOVib9QjYi4kNDaSveXQ9WPKlDarniCWZ3HHt2Ehr+hZxbweh+TAYcrQp7 +EDGLeEK/7IZVGoXuWE72iXnqB9aOEa/qgZ2cpZjGmvhMe2kX4vQwy33kLvU L/qgFWS2fXMEJn0tjsZPDunRCcC94Emb1wgf/6xPR6tlh1fQvITX6d2bBDnp NZJUqdDUi6jQREtD+2lFdVoPSQ0KBVIF8zxsIk7QCeyb1zO3QiMr8kKzc6Mq Lyn+ogigYGqSteD1jbIvzYUUtV6Hg0esSzdIFo82dxEQmuxFaegknmDj42N+ 3ds2WFRiyWMhbEOT5ruUbRJDoAKUJNJD0/njAMaLgUGotlSwXrA7rTdyj1hA aFlkjwpAaxFVzKRoElq62asCkOmielrB2tA+FYAmHpXKpH4VWgJ7LVARxE+q baEFu+TaBhKARkQqgqF1yUMqAIlHRiNVy9Da6WEVgLDI7KWyGloiPqICkEo5 AkLL2G8CGJpa5HqkVl0HKDRDx7QvSABymVJPR5OHbOyECkBTi/YXqfmjdYB0 fEsFoOFAW7DcSxjN7frbgcwGlYvc1gi98HFKGyK2KEgWAaFXUd5VASjCG1Cs XFQJvSRzWhuioxfKUIkA5HzR/E+pAMSsD8BeNQI2a0PE1oyyNxGO7sGhi1o1 FWKccofbT0YTnCY/qC2RoIpi5d5U6E0tM2ltSfTXlV8FNEIHdZGbf6HMmoib Yk8cJEeD3frlXjzLv9ZWk05k8glWdTcHj5bKuDKPSgYSpWtKLkwuziUBWi/g aGFu/sBq4ybx+7XHfJWRo3rjz2GohubCimrPFeAUf1UF7cmKTII3JKhPmBIs mwFg/zDpsErCtHIOK3RYZ7PphfoloRkq+Fqm6MTNi62ifrfqU8wMXLoUHvMz ttGjqR9rZ2IXqIAetypJK6ciyY+0Nhusc711QDqATunLlIr3o5IOQil/jbZX f91rJp0vbCTt7bDa9s7vVHubyTTbuxMh2gyxmweYVGyom0EzFCX3O8hrZlBt 0yNa2ybujRbbZI2jfKZlFOmu0inN8Ib2Dm3adpG6smlsaEBLuJeoyo9uA7hj n+HSyphMW67isJKwzbF+sZZyUSA1SrbPNOhysSSGxX0I05zPysj0UuqnUcnG cFK/TKhYyXaTLZGmh9YPEmOdiFRs2laCrYrZFPRc7QdMe0m9Gn+hkNJjAl3C xF2K9IgWBr2Y0eD6iUuKVnLB2eB9ly1kv36hI0puX9LXOIyrraX6Rjh0acVe Ccp6fHVJJcfPIv15r2G1xvXLhAnFrQkddnOGLql2b9GvDfoVd4pmK7vg0JEp RwYvjH71P1wBu+wn53TLEgLd8hP9S/7PzrxQibEKgI8EHBQ4wR4VPH83P3iP rrgm6xbjf3OIelxn55bWH/qI3PL4w6x+bkbZlr1dIrZYG8i45WwPUVuZ+oGh 49NJKZGI68wusRPb4VpxG3HiutsLug86ztkz9Cda9DdCtKQN\ \>"], "Graphics", Evaluatable->False, ImageSize->{235, 183}, ImageMargins->{{0, 0}, {0, 0}}, ImageRegion->{{0, 1}, {0, 1}}], Cell[TextData[{ "According to Sir Issac Newton, the force between a planet and the Sun is\n\t\ ", Cell[BoxData[ FormBox[ SubscriptBox["F", "G"], TraditionalForm]]], " = ", StyleBox["G ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["M", "s"], TraditionalForm]]], StyleBox["m / ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SuperscriptBox["r", "2"], TraditionalForm]]], "\nwhere ", StyleBox["m", FontSlant->"Italic"], " is the mass of the planet and ", Cell[BoxData[ FormBox[ SuperscriptBox["r", "2"], TraditionalForm]]], "= ", Cell[BoxData[ FormBox[ SuperscriptBox["x", "2"], TraditionalForm]]], "+ ", Cell[BoxData[ FormBox[ SuperscriptBox["y", "2"], TraditionalForm]]], " is the square of the distance between the planet and the Sun. \nThe force \ of gravity always attracts one object towards the other. Force is a vector, \ and as the planet orbits the Sun, the direction of the force also revolves. \ From the diagram above, the ", StyleBox["x", FontSlant->"Italic"], "-directed force is equal to the total force, ", Cell[BoxData[ FormBox[ SubscriptBox["F", "G"], TraditionalForm]]], ", times the ratio, ", StyleBox["x/r", FontSlant->"Italic"], ". The ", StyleBox["y", FontSlant->"Italic"], "-component of the force is ", StyleBox["y/r", FontSlant->"Italic"], " times ", Cell[BoxData[ FormBox[ SubscriptBox["F", "G"], TraditionalForm]]], ".\nNewton also explained the force on the Sun due to the planet is equal \ and opposite to the force on the planet from the Sun. However, if the Sun is \ much, much more massive than the planet, the Sun does not significantly move \ as compared to the planet. In this approximation, we need only follow the \ trajectory of the planet. (Physics student's should note: Newton also showed \ how to describe the motion of ", StyleBox["any", FontSlant->"Italic"], " two objects interacting through gravity, like binary stars, and this \ description is often called Kepler's ammended laws.)\nIn the ", StyleBox["x-y", FontSlant->"Italic"], " plane, the equations of motion of a planet about a massive star are:\n\t", Cell[BoxData[ FormBox[ RowBox[{ SuperscriptBox["x", "\[Prime]"], "[", "t", "]"}], TraditionalForm]]], " = vx[t]\n\t", Cell[BoxData[ FormBox[ RowBox[{ SuperscriptBox["y", "\[Prime]"], "[", "t", "]"}], TraditionalForm]]], " = vy[t]\n\t", Cell[BoxData[ FormBox[ RowBox[{ SuperscriptBox["vx", "\[Prime]"], "[", "t", "]"}], TraditionalForm]]], " = - ", StyleBox["G", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["M", "S"], TraditionalForm]]], " ", StyleBox["x / ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SuperscriptBox["r", "3"], TraditionalForm]]], "\n\t", Cell[BoxData[ FormBox[ RowBox[{ SuperscriptBox["vy", "\[Prime]"], "[", "t", "]"}], TraditionalForm]]], " = - ", StyleBox["G", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["M", "S"], TraditionalForm]]], " ", StyleBox["y / ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SuperscriptBox["r", "3"], TraditionalForm]]], "\nwhere the prime denotes teh first derivative with respect to time, ", StyleBox["df /dt = ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SuperscriptBox["f", "\[Prime]"], TraditionalForm]]], ".\nThe product of Newton's gravitational constant and the mass of the sun, \ ", StyleBox["G", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["M", "S"], TraditionalForm]]], ", takes on a simple form if we measure distance in ", StyleBox["astronomical units", FontSlant->"Italic"], ", or AU, and if we measure time in years. In these units, ", StyleBox["G", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["M", "S"], TraditionalForm]]], " = 4", Cell[BoxData[ FormBox[ SuperscriptBox["\[Pi]", "2"], TraditionalForm]]], " ", Cell[BoxData[ FormBox[ SuperscriptBox["AU", "3"], TraditionalForm]]], "/", Cell[BoxData[ FormBox[ SuperscriptBox["yr", "2"], TraditionalForm]]], ". One AU is equal to the average distance between the Sun and the Earth, \ about 1.5 \[Cross] ", Cell[BoxData[ FormBox[ SuperscriptBox["10", "11"], TraditionalForm]]], "m." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Equations of Motion", "Section"], Cell["\<\ The equations of motion for a planet describe the time evolution of the \ planet's phase-space location and its initial position.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"p", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", " ", RowBox[{"y", "[", "t", "]"}], ",", " ", RowBox[{"vx", "[", "t", "]"}], ",", " ", RowBox[{"vy", "[", "t", "]"}]}], "}"}]}], ";"}]], "Input", CellLabel->"In[1]:="], Cell[BoxData[ RowBox[{"dpdt", " ", "=", " ", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{"vx", "[", "t", "]"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"y", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{"vy", "[", "t", "]"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"-", "4"}], " ", SuperscriptBox["\[Pi]", "2"], " ", RowBox[{ RowBox[{"x", "[", "t", "]"}], "/", RowBox[{"r", "^", "3"}]}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"vy", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"-", "4"}], " ", SuperscriptBox["\[Pi]", "2"], " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], "/", RowBox[{"r", "^", "3"}]}]}]}]}], "}"}], " ", "/.", " ", RowBox[{"r", " ", "\[Rule]", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], "^", "2"}]}], "]"}]}]}]}]], "Input", CellLabel->"In[2]:="], Cell[BoxData[ RowBox[{"p0", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "0", "]"}], " ", "\[Equal]", " ", "x0"}], ",", " ", RowBox[{ RowBox[{"y", "[", "0", "]"}], " ", "\[Equal]", " ", "y0"}], ",", RowBox[{ RowBox[{"vx", "[", "0", "]"}], " ", "\[Equal]", " ", "vx0"}], ",", " ", RowBox[{ RowBox[{"vy", "[", "0", "]"}], " ", "\[Equal]", " ", "vy0"}]}], "}"}]}]], "Input", CellLabel->"In[3]:="], Cell[BoxData[ RowBox[{"eqs", " ", "=", " ", RowBox[{"dpdt", " ", "~", " ", "Join", " ", "~", "p0"}]}]], "Input", CellLabel->"In[4]:="], Cell[CellGroupData[{ Cell["Earth's Orbit", "Subsection"], Cell[TextData[{ "In astronomical units (AUand year), the trajectory of the Earth is \ relatively easy to calculate. The key is finding the appropriate initial \ condition. \nLet's set the position of Earth to be {", StyleBox["x, y", FontSlant->"Italic"], "} = {1, 0}. The initial velocity is approximately the circumference, \ 2\[Pi]", StyleBox["r", FontSlant->"Italic"], " = 2\[Pi]", StyleBox[", ", FontSlant->"Italic"], "dividedby 1 year, or {", StyleBox["vx, vy", FontSlant->"Italic"], "} = {0, 2\[Pi]}." }], "Text"], Cell[TextData[{ "In Cartesian coordinates, the equations of motion are two difficult for ", StyleBox["Mathematica", FontSlant->"Italic"], " to solve using ", StyleBox["DSolve", FontWeight->"Bold"], ". (If we do some algebra and transformation of variables, we can use ", StyleBox["DSolve", FontWeight->"Bold"], " for the two-body central-force problem, but we'll discuss this later.)" }], "Text"], Cell[BoxData[ RowBox[{"earth", " ", "=", " ", RowBox[{"DSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", "\[Pi]"}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", "t"}], "]"}]}]], "Input", CellLabel->"In[5]:="], Cell[TextData[{ "Instead, we must use ", StyleBox["NDSolve", FontWeight->"Bold"], "." }], "Text"], Cell[BoxData[ RowBox[{"earth", " ", "=", " ", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", "\[Pi]"}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "1.2"}], "}"}]}], "]"}]}]], "Input",\ CellLabel->"In[6]:="], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earth"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "1"}], "}"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input", CellLabel->"In[7]:="], Cell["\<\ The orbit is a circle because of our choice of the initial condition.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"Off", "[", RowBox[{ RowBox[{"General", "::", "spell"}], ",", " ", RowBox[{"General", "::", "spell1"}]}], "]"}], ";"}]], "Input", CellLabel->"In[8]:="], Cell[BoxData[{ RowBox[{ RowBox[{"earthA", " ", "=", " ", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", RowBox[{"\[Pi]", "/", "1.1"}]}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "1.2"}], "}"}]}], "]"}]}], ";"}], "\n", RowBox[{ RowBox[{"earthB", " ", "=", " ", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", RowBox[{"\[Pi]", "/", "1.2"}]}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "1.2"}], "}"}]}], "]"}]}], ";"}], "\n", RowBox[{ RowBox[{"earthC", " ", "=", " ", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", "\[Pi]", " ", "1.1"}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "1.2"}], "}"}]}], "]"}]}], ";"}], "\n", RowBox[{ RowBox[{"earthD", " ", "=", " ", RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "1"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", RowBox[{"2", " ", "\[Pi]", " ", "1.2"}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "1.2"}], "}"}]}], "]"}]}], ";"}]}], "Input", CellLabel->"In[9]:="], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earth"}], ",", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earthA"}], ",", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earthB"}], ",", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earthC"}], ",", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "earthD"}]}], "}"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "1"}], "}"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", "}"}], ",", RowBox[{"{", RowBox[{"Hue", "[", FractionBox["0", "4"], "]"}], "}"}], ",", RowBox[{"{", RowBox[{"Hue", "[", FractionBox["1", "4"], "]"}], "}"}], ",", RowBox[{"{", RowBox[{"Hue", "[", FractionBox["2", "4"], "]"}], "}"}], ",", RowBox[{"{", RowBox[{"Hue", "[", FractionBox["3", "4"], "]"}], "}"}]}], "}"}]}], ",", RowBox[{"Compiled", "\[Rule]", "False"}]}], "]"}]], "Input", CellLabel->"In[13]:="] }, Closed]], Cell[CellGroupData[{ Cell["Halley's Comet", "Subsection"], Cell["\<\ Halley's coment has anorbital period of 76 yr and a closest distance tothe \ Sun of 0.59 AU. The challenge is finding the initial condition which determines the orbit. Below, we try the \"brute force\" or \"rial-and-error\" approach--perfect for \ computers...\ \>", "Text"], Cell[BoxData[ RowBox[{"halley", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "0.59"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", "11.46"}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "76"}], "}"}], ",", RowBox[{"MaxSteps", " ", "->", " ", "2000"}]}], "]"}], " ", "//", " ", "Flatten"}]}]], "Input", CellLabel->"In[14]:="], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}]}], "]"}]], "Input", CellLabel->"In[15]:="], Cell[BoxData[ RowBox[{"Show", "[", RowBox[{"%", ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input", CellLabel->"In[16]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"vx", "[", "t", "]"}], ",", RowBox[{"vy", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"RGBColor", "[", RowBox[{"1", ",", "0", ",", "0"}], "]"}], "}"}], ",", RowBox[{"{", RowBox[{"RGBColor", "[", RowBox[{"0", ",", "1", ",", "0"}], "]"}], "}"}]}], "}"}]}]}], "]"}]], "Input", CellLabel->"In[17]:="], Cell["\<\ Notice that the orbital period is not quite 76 years. We'll try forming a \ table of possible orbits.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"halleyTable", " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "0.59"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", "v"}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "76"}], "}"}], ",", RowBox[{"MaxSteps", " ", "->", " ", "2000"}]}], "]"}], " ", "//", " ", "Flatten"}], ",", " ", RowBox[{"{", RowBox[{"v", ",", " ", "11.46", ",", " ", "11.48", ",", " ", "0.004"}], "}"}]}], "]"}]}], ";"}]], "Input", CellLabel->"In[18]:="], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "halleyTable"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "36", ",", "40"}], "}"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{"Hue", "[", FractionBox["i", "5"], "]"}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "0", ",", "4"}], "}"}]}], "]"}]}]}], "]"}]], "Input",\ CellLabel->"In[19]:="], Cell[TextData[{ "So the velocity of Halley's comet at its closest distance to the Sun is \ about ", StyleBox["v", FontSlant->"Italic"], " \[TildeTilde] 11.472 AU/yr." }], "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["The Ellipse (optional)", "Section"], Cell[TextData[{ "(This section is optional. Students interested in techniques to use ", StyleBox["Mathematica", FontSlant->"Italic"], " to manipulate algebraic equations can look at the expressions below.)" }], "Text"], Cell[TextData[{ "An ellipse is is a curve such that the sum of the distances between each \ point on the ellipse to two fixed points (", StyleBox["i.e.", FontSlant->"Italic"], " the foci) are constant." }], "Text"], Cell[CellGroupData[{ Cell["A Centered, Symmetric Ellipse", "Subsection"], Cell[TextData[{ "If the two foci are located at {", StyleBox["x, y", FontSlant->"Italic"], "}", StyleBox[" ", FontSlant->"Italic"], "= {\[PlusMinus]", StyleBox["c, ", FontSlant->"Italic"], "0}, and if the distance between the ellipse when it crosses the ", StyleBox["x", FontSlant->"Italic"], " axes is ", StyleBox["a, ", FontSlant->"Italic"], "then the equation for the ellipse is" }], "Text"], Cell[BoxData[ RowBox[{"ellipse1", " ", "=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{"x", "-", "c"}], ")"}], "2"], "+", " ", SuperscriptBox["y", "2"]}], "]"}], " ", "+", " ", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{"x", "+", "c"}], ")"}], "2"], "+", " ", SuperscriptBox["y", "2"]}], "]"}]}], " ", "==", " ", RowBox[{"2", " ", "a"}]}]}]], "Input", CellLabel->"In[20]:="], Cell[TextData[{ "This can be turned into a more familiar form of the ellipse by replacing ", Cell[BoxData[ FormBox[ SuperscriptBox["c", "2"], TraditionalForm]]], " with the difference between the major and minor axes." }], "Text"], Cell[BoxData[ RowBox[{"eSol1", " ", "=", " ", RowBox[{ RowBox[{ RowBox[{"Solve", "[", RowBox[{"ellipse1", ",", "y"}], "]"}], "//.", " ", RowBox[{ SuperscriptBox["c", "2"], "->", RowBox[{ SuperscriptBox["a", "2"], " ", "-", " ", RowBox[{"b", "^", "2"}]}]}]}], " ", "//", " ", "Simplify", StyleBox[" ", FontWeight->"Plain"]}]}]], "Input", CellLabel->"In[21]:="], Cell[BoxData[ RowBox[{"ellipse2", " ", "=", " ", RowBox[{ RowBox[{ RowBox[{ SuperscriptBox["x", "2"], "/", SuperscriptBox["a", "2"]}], " ", "+", " ", RowBox[{ SuperscriptBox["y", "2"], "/", SuperscriptBox["b", "2"]}]}], " ", "==", " ", "1"}]}]], "Input", CellLabel->"In[22]:="], Cell[BoxData[ RowBox[{"eSol2", " ", "=", " ", RowBox[{ RowBox[{"Solve", "[", RowBox[{"ellipse2", ",", " ", "y"}], "]"}], " ", "//", " ", "Simplify"}]}]], "Input", CellLabel->"In[23]:="], Cell[BoxData[ RowBox[{ RowBox[{"Expand", "[", RowBox[{ RowBox[{"y", "^", "2"}], " ", "/.", " ", RowBox[{"eSol1", "[", RowBox[{"[", "1", "]"}], "]"}]}], "]"}], " ", "==", " ", RowBox[{"Expand", "[", RowBox[{ RowBox[{"y", "^", "2"}], "/.", RowBox[{"eSol2", "[", RowBox[{"[", "1", "]"}], "]"}]}], "]"}], " ", RowBox[{"(*", " ", "QED", " ", "*)"}]}]], "Input", CellLabel->"In[24]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{"y", "/.", "\[InvisibleSpace]", "eSol1"}], "/.", "\[InvisibleSpace]", RowBox[{"b", "\[Rule]", FractionBox["1", "2"]}]}], "/.", "\[InvisibleSpace]", RowBox[{"a", "\[Rule]", "1"}]}], "]"}], ",", RowBox[{"{", RowBox[{"x", ",", RowBox[{"-", "1"}], ",", "1"}], "}"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input", CellLabel->"In[25]:="], Cell[TextData[{ "The eccentricity, ", StyleBox["e = c / a", FontSlant->"Italic"], ", is defined to be the ratio of ", StyleBox["c", FontSlant->"Italic"], " to the semimajor axes, ", StyleBox["a.", FontSlant->"Italic"], " (This assumes the ordering, ", StyleBox["a > b", FontSlant->"Italic"], " and 0", StyleBox[" \[GreaterEqual] ", FontSlant->"Italic"], "c", StyleBox[" \[GreaterEqual] a.", FontSlant->"Italic"], ") From the replacements demonstrated above, the eccentricity is also given \ by the equation ", Cell[BoxData[ FormBox[ SuperscriptBox["b", "2"], TraditionalForm]]], "= ", Cell[BoxData[ FormBox[ SuperscriptBox["a", "2"], TraditionalForm]]], " (1 - ", Cell[BoxData[ FormBox[ SuperscriptBox["e", "2"], TraditionalForm]]], "), or ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ SuperscriptBox["a", "2"], SuperscriptBox["e", "2"]}], "=", " ", RowBox[{ SuperscriptBox["a", "2"], " ", "-", " ", SuperscriptBox["b", "2"]}]}], TraditionalForm]]], ".When ", StyleBox["e", FontSlant->"Italic"], " = 1 (and ", StyleBox["c", FontSlant->"Italic"], " = ", StyleBox["a", FontSlant->"Italic"], "), the ellipse becomes a line. When ", StyleBox["e", FontSlant->"Italic"], "= ", StyleBox["c = ", FontSlant->"Italic"], "0, the ellipse becomes a circle." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["One Foci at the Origin", "Subsection"], Cell["\<\ The orbit of a planet forms an ellipse with the Sun at one of the foci.\ \>", "Text"], Cell[BoxData[ RowBox[{"ellipse2", " ", "=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox["x", "2"], "+", " ", SuperscriptBox["y", "2"]}], "]"}], " ", "+", " ", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox[ RowBox[{"(", RowBox[{"x", "+", RowBox[{"2", "c"}]}], ")"}], "2"], "+", " ", SuperscriptBox["y", "2"]}], "]"}]}], " ", "==", " ", RowBox[{"2", " ", "a"}]}]}]], "Input", CellLabel->"In[26]:="], Cell[BoxData[ RowBox[{"eSol3", " ", "=", " ", RowBox[{ RowBox[{ RowBox[{"Solve", "[", RowBox[{"ellipse2", ",", "y"}], "]"}], "//.", " ", RowBox[{"c", "->", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox["a", "2"], " ", "-", " ", RowBox[{"b", "^", "2"}]}], "]"}]}]}], " ", "//", " ", "Simplify", StyleBox[" ", FontWeight->"Plain"]}]}]], "Input", CellLabel->"In[27]:="], Cell[BoxData[ RowBox[{"c1", " ", "=", " ", RowBox[{"Sqrt", "[", RowBox[{ SuperscriptBox["1", "2"], "-", SuperscriptBox[ RowBox[{"(", RowBox[{"1", "/", "2"}], ")"}], "2"]}], "]"}], " ", RowBox[{"(*", " ", RowBox[{ RowBox[{"An", " ", "example", " ", "for", " ", "a"}], " ", "=", " ", RowBox[{ RowBox[{"1", " ", "and", " ", "b"}], " ", "=", " ", RowBox[{"1", "/", "2"}]}]}], " ", "*)"}]}]], "Input", CellLabel->"In[28]:="], Cell[BoxData[ RowBox[{ RowBox[{"N", "[", "c1", "]"}], " ", RowBox[{"(*", " ", RowBox[{"c", " ", "=", " ", RowBox[{ RowBox[{"e", " ", "when", " ", "a"}], " ", "=", " ", "1"}]}], " ", "*)"}]}]], "Input", CellLabel->"In[29]:="], Cell[BoxData[ RowBox[{ RowBox[{"{", RowBox[{"c1", ",", " ", RowBox[{"2", "c1"}], ",", " ", RowBox[{ RowBox[{"-", "1"}], " ", "-", "c1"}]}], "}"}], " ", "//", " ", "N"}]], "Input", CellLabel->"In[30]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{"y", "/.", "\[InvisibleSpace]", "eSol3"}], "/.", "\[InvisibleSpace]", RowBox[{"b", "\[Rule]", FractionBox["1", "2"]}]}], "/.", "\[InvisibleSpace]", RowBox[{"a", "\[Rule]", "1"}]}], "]"}], ",", RowBox[{"{", RowBox[{"x", ",", RowBox[{ RowBox[{"-", "1"}], "-", "c1"}], ",", RowBox[{"1", "-", "c1"}]}], "}"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input", CellLabel->"In[31]:="] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{"For", " ", "an", " ", "ellipse"}], ",", " ", RowBox[{ SuperscriptBox["y", "2"], " ", "varies", " ", "as", " ", "a", " ", "quadratic", " ", "in", " ", "x"}]}]], "Subsection"], Cell[TextData[{ "For an ellipse, the variation of ", Cell[BoxData[ FormBox[ SuperscriptBox["y", "2"], TraditionalForm]]], "is given by a quadratic polynomial in ", StyleBox["x.", FontSlant->"Italic"] }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"y", "^", "2"}], " ", "/.", " ", "eSol3"}], " ", "//", " ", RowBox[{ RowBox[{"Collect", "[", RowBox[{"#", ",", "x"}], "]"}], "&"}]}]], "Input", CellLabel->"In[32]:="] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Kepler's First Law", "Section"], Cell[TextData[{ "Since we have the numerical solution to Newton's equations of motion, a \ relatively easy method to prove the solution is represented by an ellipse is \ to perform a ", StyleBox["least-squares fit", FontSlant->"Italic"], " between our (approximate numerical) solutions and a perfect ellipse.\nFor \ an ellipse, the variation of ", Cell[BoxData[ FormBox[ SuperscriptBox["y", "2"], TraditionalForm]]], "is given by a quadratic polynomial in ", StyleBox["x.", FontSlant->"Italic"], " (See last section.)" }], "Text"], Cell["\<\ Making a list of the coordinates of Halley once every year...\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"halleyEllipse", " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{ RowBox[{"y", "[", "t", "]"}], "^", "2"}]}], "}"}], " ", "/.", " ", "halley"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", RowBox[{"76", "/", "2"}]}], "}"}]}], "]"}]}], ";"}]], "Input", CellLabel->"In[33]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{"halleyEllipse", ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"PointSize", "[", "0.02`", "]"}]}], ",", RowBox[{"AxesLabel", "\[Rule]", RowBox[{"{", RowBox[{"\"\\"", ",", "\"\\""}], "}"}]}]}], "]"}]], "Input", CellLabel->"In[34]:="], Cell["The best fit ellipse to the numerical solution is", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"yfit", "[", "x_", "]"}], " ", "=", " ", RowBox[{"Sqrt", "[", RowBox[{"Fit", "[", RowBox[{"halleyEllipse", ",", " ", RowBox[{"{", RowBox[{"1", ",", "x", ",", RowBox[{"x", "^", "2"}]}], "}"}], ",", " ", "x"}], "]"}], "]"}]}]], "Input", CellLabel->"In[35]:="], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"yfit", "[", RowBox[{"x", "[", "t", "]"}], "]"}]}], "}"}]}], "}"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", FractionBox["76", "2.`"]}], "}"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", "}"}], ",", RowBox[{"{", RowBox[{"RGBColor", "[", RowBox[{"1", ",", "0", ",", "0"}], "]"}], "}"}]}], "}"}]}]}], "]"}]], "Input", CellLabel->"In[36]:="], Cell["\<\ The error between the perfect elliptical fit and the numerical solution is \ pretty good:\ \>", "Text"], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{"y", "[", "t", "]"}], "-", RowBox[{"yfit", "[", RowBox[{"x", "[", "t", "]"}], "]"}]}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", FractionBox["76", "2"]}], "}"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[37]:="], Cell[TextData[{ "(B.T.W., you should note that ", StyleBox["Mathematica", FontSlant->"Italic"], " reports errors in the plots when the coordinates of the orbit move from ", StyleBox["y", FontSlant->"Italic"], " > 0 to ", StyleBox["y", FontSlant->"Italic"], " < 0. Do you know why?)" }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Kepler's Second Law", "Section"], Cell["\<\ Kepler's second law states: The line joining a planet to the Sun sweeps out \ equal areas in equal times.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"p1", " ", "=", " ", RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], " ", "/.", " ", "halley"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "76"}], "}"}], ",", " ", RowBox[{"DisplayFunction", "->", " ", "Identity"}]}], "]"}]}], ";"}]], "Input", CellLabel->"In[38]:="], Cell[BoxData[ RowBox[{ RowBox[{"p2", " ", "=", " ", RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], " ", "/.", " ", "halley"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", "0", ",", " ", "2"}], "}"}], ",", " ", RowBox[{"DisplayFunction", "->", " ", "Identity"}], ",", RowBox[{"PlotStyle", "->", RowBox[{"{", RowBox[{ RowBox[{"Hue", "[", "0.7", "]"}], ",", RowBox[{"Thickness", "[", "0.02", "]"}]}], "}"}]}]}], "]"}]}], ";"}]], "Input", CellLabel->"In[39]:="], Cell[BoxData[ RowBox[{ RowBox[{"p3", " ", "=", " ", RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], " ", "/.", " ", "halley"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", "35", ",", " ", "37"}], "}"}], ",", " ", RowBox[{"DisplayFunction", "->", " ", "Identity"}], ",", RowBox[{"PlotStyle", "->", RowBox[{"{", RowBox[{ RowBox[{"Hue", "[", "0.3", "]"}], ",", RowBox[{"Thickness", "[", "0.02", "]"}]}], "}"}]}]}], "]"}]}], ";"}]], "Input", CellLabel->"In[40]:="], Cell[BoxData[ RowBox[{"Show", "[", RowBox[{"p1", ",", "p2", ",", "p3", ",", RowBox[{"DisplayFunction", "\[Rule]", "$DisplayFunction"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input", CellLabel->"In[41]:="], Cell[TextData[{ "Kepler's second law is a statement of conservation of angular momentum.\n\ The angular momentum about the Sun (located at the origin) is equal to \n\t", StyleBox["L = ", FontSlant->"Italic"], StyleBox["r \[Times] v", FontWeight->"Bold"], "\n\t", StyleBox["L = x ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["v", "y"], TraditionalForm]]], " - ", StyleBox["y ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SubscriptBox["v", "x"], TraditionalForm]]], "\nExpressed in this way (without the mass of the planet), ", StyleBox["L", FontSlant->"Italic"], " is equal to the incremental increase of the area spanned by the line \ connecting the planet to the star per incremental time." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"l", "[", "t_", "]"}], " ", "=", " ", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], " ", RowBox[{"vy", "[", "t", "]"}]}], " ", "-", " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], " ", RowBox[{"vx", "[", "t", "]"}]}]}], " ", "/.", " ", "halley"}]}]], "Input",\ CellLabel->"In[42]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"l", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}]}], "]"}]], "Input", CellLabel->"In[43]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ FractionBox[ RowBox[{"100", " ", RowBox[{"(", RowBox[{ RowBox[{"l", "[", "t", "]"}], "-", RowBox[{"l", "[", "0", "]"}]}], ")"}]}], RowBox[{"l", "[", "0", "]"}]], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}], ",", RowBox[{ "PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[44]:="], Cell[TextData[{ "The total area of the ellipse is given by the integral of ", StyleBox["L", FontSlant->"Italic"], " around one complete orbit, ", StyleBox["A", FontSlant->"Italic"], " = \[Pi]", StyleBox["ab.", FontSlant->"Italic"] }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Kepler's Third Law", "Section"], Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ RowBox[{"Kepler", "'"}], "s", " ", "third", " ", "law", " ", RowBox[{"states", ":", " ", RowBox[{"If", " ", StyleBox["T", FontSlant->"Italic"], " ", "is", " ", "equal", " ", "to", " ", "the", " ", "orbital", " ", "period", " ", "of", " ", "a", " ", "planet", " ", RowBox[{"(", RowBox[{ RowBox[{ StyleBox["i", FontSlant->"Italic"], StyleBox[".", FontSlant->"Italic"], StyleBox["e", FontSlant->"Italic"], StyleBox[".", FontSlant->"Italic"], " ", "the"}], " ", "time", " ", "required", " ", "to", " ", "complete", " ", "one", " ", "orbit"}], ")"}]}]}]}], ",", " ", RowBox[{"and", " ", "if", " ", StyleBox["a", FontSlant->"Italic"], " ", "is", " ", "the", " ", "semimajor", " ", "axis", " ", "of", " ", "the", " ", "orbit", RowBox[{"(", RowBox[{ RowBox[{ StyleBox["i", FontSlant->"Italic"], StyleBox[".", FontSlant->"Italic"], StyleBox["e", FontSlant->"Italic"], StyleBox[".", FontSlant->"Italic"], " ", "one"}], " ", "half", " ", "the", " ", "long", " ", "distance", " ", "across", " ", "the", " ", "ellispe"}], ")"}]}], ",", " ", RowBox[{"then", " ", "the", " ", "ratio", " ", RowBox[{ FormBox[ SuperscriptBox["T", "2"], TraditionalForm], "/", FormBox[ SuperscriptBox["a", "3"], TraditionalForm]}], " ", "is", " ", "a", " ", RowBox[{"constant", "."}]}]}], TextForm]], "Text"], Cell["\<\ To prove Kepler's third law, we must develop a way to find the period of a \ planet's orbit and the length of its major axis.\ \>", "Text"], Cell[CellGroupData[{ Cell["Harmonic Constant for Halley", "Subsection"], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}]}], "]"}]], "Input", CellLabel->"In[45]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "60", ",", "70"}], "}"}]}], "]"}]], "Input", CellLabel->"In[46]:="], Cell[BoxData[ RowBox[{ RowBox[{"tHalley", " ", "=", " ", RowBox[{"t", " ", "/.", " ", RowBox[{"FindRoot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "halley"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "63"}], "}"}]}], "]"}]}]}], ";"}]], "Input", CellLabel->"In[47]:="], Cell[BoxData[ RowBox[{ RowBox[{"aHalley", " ", "=", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"(", " ", RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], "/.", " ", "halley"}], " ", "/.", RowBox[{"t", "->", "0"}]}], ")"}], "-", " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], "/.", " ", "halley"}], "/.", RowBox[{"t", "->", RowBox[{"tHalley", "/", "2"}]}]}], ")"}]}], ")"}], "/", "2"}]}], " ", ";"}]], "Input", CellLabel->"In[48]:="] }, Closed]], Cell[CellGroupData[{ Cell["Harmonic Constant for Earth", "Subsection"], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "earth"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "1"}], "}"}]}], "]"}]], "Input", CellLabel->"In[49]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "earth"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0.8`", ",", "1"}], "}"}]}], "]"}]], "Input", CellLabel->"In[50]:="], Cell[BoxData[ RowBox[{ RowBox[{"tEarth", " ", "=", " ", "1.0"}], ";", " ", RowBox[{"aEarth", " ", "=", " ", "1.0"}], ";"}]], "Input", CellLabel->"In[51]:="], Cell[BoxData[ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"tHalley", "^", "2"}], "/", RowBox[{"aHalley", "^", "3"}]}], ",", " ", RowBox[{ RowBox[{"tEarth", "^", "2"}], "/", RowBox[{"aEarth", "^", "3"}]}]}], "}"}], " ", "//", " ", "N"}]], "Input",\ CellLabel->"In[52]:="] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Summary", "Section"], Cell[TextData[{ "We've used ", StyleBox["NDSolve[\[Ellipsis]]", FontWeight->"Bold"], " to numerically integrate the equations of motion for an object orbiting a \ massive star. These solutions depend sensitively on the ", StyleBox["initial conditions", FontSlant->"Italic"], " of the object. We examined in detail the orbits of the planet Earth and \ the comet Halley. We showed that the solutions produced by ", StyleBox["Mathematica", FontSlant->"Italic"], " are in agreement with the three Laws of Kepler." }], "Text"] }, Closed]] }, Open ]] }, WindowSize->{596, 821}, WindowMargins->{{2, Automatic}, {Automatic, 1}}, PrintingCopies->1, PrintingPageRange->{1, Automatic}, FrontEndVersion->"6.0 for Mac OS X x86 (32-bit) (May 21, 2008)", StyleDefinitions->"TutorialBook.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[590, 23, 30, 0, 51, "Title"], Cell[623, 25, 49, 3, 41, "Subsubtitle"], Cell[CellGroupData[{ Cell[697, 32, 25, 0, 86, "Section"], Cell[725, 34, 3815, 67, 809, "Text"], Cell[CellGroupData[{ Cell[4565, 105, 41, 0, 34, "Subsection"], Cell[4609, 107, 932, 26, 167, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[5590, 139, 56, 0, 54, "Section"], Cell[5649, 141, 3453, 59, 191, "Graphics", Evaluatable->False], Cell[9105, 202, 4271, 151, 584, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[13413, 358, 38, 0, 54, "Section"], Cell[13454, 360, 153, 3, 41, "Text"], Cell[13610, 365, 310, 9, 24, "Input"], Cell[13923, 376, 1498, 46, 114, "Input"], Cell[15424, 424, 464, 13, 42, "Input"], Cell[15891, 439, 139, 3, 26, "Input"], Cell[CellGroupData[{ Cell[16055, 446, 35, 0, 34, "Subsection"], Cell[16093, 448, 536, 17, 90, "Text"], Cell[16632, 467, 409, 11, 59, "Text"], Cell[17044, 480, 491, 14, 60, "Input"], Cell[17538, 496, 102, 5, 23, "Text"], Cell[17643, 503, 563, 16, 60, "Input"], Cell[18209, 521, 442, 13, 45, "Input"], Cell[18654, 536, 93, 2, 23, "Text"], Cell[18750, 540, 199, 6, 24, "Input"], Cell[18952, 548, 2357, 67, 284, "Input"], Cell[21312, 617, 1850, 56, 197, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[23199, 678, 36, 0, 34, "Subsection"], Cell[23238, 680, 285, 6, 83, "Text"], Cell[23526, 688, 650, 17, 95, "Input"], Cell[24179, 707, 387, 12, 44, "Input"], Cell[24569, 721, 152, 4, 26, "Input"], Cell[24724, 727, 693, 22, 78, "Input"], Cell[25420, 751, 125, 3, 23, "Text"], Cell[25548, 756, 884, 23, 113, "Input"], Cell[26435, 781, 657, 21, 120, "Input"], Cell[27095, 804, 184, 6, 23, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[27328, 816, 41, 0, 54, "Section"], Cell[27372, 818, 223, 5, 41, "Text"], Cell[27598, 825, 219, 6, 41, "Text"], Cell[CellGroupData[{ Cell[27842, 835, 51, 0, 34, "Subsection"], Cell[27896, 837, 415, 17, 42, "Text"], Cell[28314, 856, 536, 17, 52, "Input"], Cell[28853, 875, 238, 6, 41, "Text"], Cell[29094, 883, 414, 13, 48, "Input"], Cell[29511, 898, 316, 10, 35, "Input"], Cell[29830, 910, 202, 6, 26, "Input"], Cell[30035, 918, 427, 13, 45, "Input"], Cell[30465, 933, 524, 15, 82, "Input"], Cell[30992, 950, 1358, 58, 78, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[32387, 1013, 44, 0, 34, "Subsection"], Cell[32434, 1015, 95, 2, 23, "Text"], Cell[32532, 1019, 504, 16, 51, "Input"], Cell[33039, 1037, 424, 13, 68, "Input"], Cell[33466, 1052, 477, 14, 51, "Input"], Cell[33946, 1068, 254, 8, 24, "Input"], Cell[34203, 1078, 228, 8, 26, "Input"], Cell[34434, 1088, 579, 17, 82, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[35050, 1110, 218, 5, 38, "Subsection"], Cell[35271, 1117, 222, 8, 23, "Text"], Cell[35496, 1127, 224, 7, 26, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[35769, 1140, 37, 0, 54, "Section"], Cell[35809, 1142, 546, 15, 89, "Text"], Cell[36358, 1159, 85, 2, 23, "Text"], Cell[36446, 1163, 480, 15, 44, "Input"], Cell[36929, 1180, 317, 8, 62, "Input"], Cell[37249, 1190, 65, 0, 23, "Text"], Cell[37317, 1192, 329, 10, 28, "Input"], Cell[37649, 1204, 870, 28, 138, "Input"], Cell[38522, 1234, 113, 3, 23, "Text"], Cell[38638, 1239, 473, 14, 83, "Input"], Cell[39114, 1255, 307, 11, 41, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[39458, 1271, 38, 0, 54, "Section"], Cell[39499, 1273, 129, 3, 41, "Text"], Cell[39631, 1278, 534, 16, 63, "Input"], Cell[40168, 1296, 711, 21, 81, "Input"], Cell[40882, 1319, 713, 21, 81, "Input"], Cell[41598, 1342, 244, 5, 45, "Input"], Cell[41845, 1349, 753, 24, 166, "Text"], Cell[42601, 1375, 363, 12, 28, "Input"], Cell[42967, 1389, 193, 6, 24, "Input"], Cell[43163, 1397, 457, 15, 87, "Input"], Cell[43623, 1414, 254, 10, 24, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[43914, 1429, 37, 0, 54, "Section"], Cell[43954, 1431, 1640, 49, 85, "Text"], Cell[45597, 1482, 149, 3, 41, "Text"], Cell[CellGroupData[{ Cell[45771, 1489, 50, 0, 34, "Subsection"], Cell[45824, 1491, 289, 9, 26, "Input"], Cell[46116, 1502, 290, 9, 26, "Input"], Cell[46409, 1513, 379, 11, 60, "Input"], Cell[46791, 1526, 558, 18, 61, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[47386, 1549, 49, 0, 34, "Subsection"], Cell[47438, 1551, 287, 9, 26, "Input"], Cell[47728, 1562, 290, 9, 26, "Input"], Cell[48021, 1573, 163, 4, 24, "Input"], Cell[48187, 1579, 314, 11, 44, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[48550, 1596, 26, 0, 54, "Section"], Cell[48579, 1598, 534, 13, 78, "Text"] }, Closed]] }, Open ]] } ] *) (* End of internal cache information *)