(************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 37187, 878]*) (*NotebookOutlinePosition[ 37850, 901]*) (* CellTagsIndexPosition[ 37806, 897]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Convert mean Anomaly to eccentric Anomaly", "Section 1"], Cell[BoxData[ \(meanAnom2eccAnom[e_, \ M_, debug_: False] := \ Module[\[IndentingNewLine]\[IndentingNewLine]{EAnext, \ EA, \ error, \ \[Epsilon] = 10\^\(-15\), i = 0}, \[IndentingNewLine]If[ e > 1, \ Return[0.0]]; \[IndentingNewLine]If[e < 0, \ Print["\"]; \ Abort[]]; \[IndentingNewLine]\[IndentingNewLine]EAnext = M; \[IndentingNewLine]error\ = \ 2*\[Epsilon]; \[IndentingNewLine]While[ error\ > \ \[Epsilon]\ \[And] \ \(i++\) < 15, \ \[IndentingNewLine]EA\ = \ EAnext; \ \[IndentingNewLine]EAnext = EA - \(EA - e*Sin[EA] - M\)\/\(1 - e*Cos[EA]\); \ \[IndentingNewLine]error\ = \ Abs[EAnext - EA]; \[IndentingNewLine]\[IndentingNewLine]If[ debug, \ Print["\", EA, "\< Next guess = \>", \ EAnext, "\< Difference =\>", error]\[IndentingNewLine]];\[IndentingNewLine]]; \ \[IndentingNewLine]Return[EA]\[IndentingNewLine]]\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Convert mean anomaly to true anomaly", "Section 1"], Cell[BoxData[ \(meanAnom2TrueAnom[e_, \ M_, \ debug_: False] := \ Module[\[IndentingNewLine]{EA, \ \[Theta]}, \[IndentingNewLine]EA\ = \ meanAnom2eccAnom[e, \ M]; \[IndentingNewLine]If[debug, \ Print["\", \ M, "\< E = \>", \ EA]]; \[IndentingNewLine]\[Theta] = \ ArcTan[\(Cos[EA] - e\)\/\(1 - e*Cos[EA]\), \ \(\(\@\(1 - e\^2\)\) \ Sin[EA]\)\/\(1 - e*Cos[EA]\)]; \[IndentingNewLine]Return[\[Theta]];\ \[IndentingNewLine]]\)], "Input"] }, Open ]], Cell[CellGroupData[{ Cell[" Convert Cartesian elements to Kepler Elements", "Section 1"], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{\(kepler2Cartesian[kepler_, \ debug_: False]\), ":=", " ", RowBox[{"Module", "[", RowBox[{\({a, \ e, \ i, \ \[CapitalOmega], \ \[Omega], \ M, \ IHAT, JHAT, \ R, \ \[Theta], \ r, x, y, z, vx, vy, vz, \[Mu] = 398600.436}\), ",", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{\(If[Length[kepler] \[NotEqual] \ 6, \ Print["\"]; \ \[IndentingNewLine]Abort[];\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \({a, \ e, \ i, \ \[CapitalOmega], \ \[Omega], \ M}\ = \ kepler\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["(*", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox[ RowBox[{ RowBox[{\(IH AT\), " ", "and", " ", "JHAT", " ", "are", " ", "unit", " ", "vectors", " ", "in", " ", "the", " ", "plane", " ", "of", " ", "the", " ", RowBox[{ StyleBox["orbit", FontWeight->"Plain"], StyleBox[":", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], RowBox[{ StyleBox["They", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["are", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["not", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["the", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["ihat", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["and", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["jhat", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["vectors", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["in", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["the", FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox[\((x, y, z)\), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox[\(frame!\), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], "IHAT", " ", "points", " ", "toward", " ", "perigee"}]}]}], ";", " ", \(JHAT\ is\ normal\ and\ in\ the\ plane\ of\ the\ \ orbit\)}], FontWeight->"Plain"], "\[IndentingNewLine]", StyleBox["*)", FontWeight->"Plain"]}], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(IHAT = \ {\[IndentingNewLine]Cos[\ \[CapitalOmega]]*Cos[\[Omega]] - Sin[\[CapitalOmega]]*Sin[\[Omega]]*Cos[i], Sin[\[CapitalOmega]]*Cos[\[Omega]] + Cos[\[CapitalOmega]]*Sin[\[Omega]]*Cos[i], Sin[\[Omega]]*Sin[i]\[IndentingNewLine]}\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(JHAT = {\[IndentingNewLine]\(-Cos[\ \[CapitalOmega]]\)*Sin[\[Omega]] - Sin[\[CapitalOmega]]*Cos[\[Omega]]* Cos[i], \[IndentingNewLine]\(-Sin[\[CapitalOmega]]\)* Sin[\[Omega]] + Cos[\[CapitalOmega]]*Cos[\[Omega]]*Cos[i], \ Cos[\[Omega]]*Sin[i]\[IndentingNewLine]}\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(\[Theta] = meanAnom2TrueAnom[e, \ M]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(r\ = \(a*\((1 - e\^2)\)\)\/\(1 + e\ \ Cos[\[Theta]]\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \({x, y, z} = r*Cos[\[Theta]]*IHAT + r*Sin[\[Theta]]*JHAT\), ";", "\[IndentingNewLine]", \({vx, vy, vz} = \(\@\(\[Mu]\/\(a \((1 - e\^2)\)\)\)\) \((\(-Sin[\[Theta]]\)* IHAT + \ \((e + Cos[\[Theta]])\)*JHAT)\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(If[ debug, \[IndentingNewLine]\(Print["\", IHAT, "\<\nJHAT=\>", JHAT, "\<\nW=\>", W, "\<\nM=\>", M, "\< \[Theta]=\>", \[Theta], "\< r=\>", r, "\<\n{x,y,z}=\>", {x, y, z}, "\<\n{vx,vy,vz}=\>", {vx, vy, vz}];\)\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(Return[{x, y, z, vx, vy, vz}]\), ";"}]}], "\[IndentingNewLine]", "\n", "]"}]}], ";"}]], "Input"], Cell[BoxData[ \(General::"spell1" \(\(:\)\(\ \)\) "Possible spelling error: new symbol name \"\!\(JHAT\)\" is similar to \ existing symbol \"\!\(IHAT\)\"."\)], "Message"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(Table[ kepler2Cartesian[{7700, 0.01, Pi/3, \ \((200/180)\)*Pi, \ 3*Pi/2, \ M*Pi/180}], {M, 0, 180, 10}]\)], "Input"], Cell[BoxData[ \({{\(-1303.609776285786`\), 3581.63842412549`, \(-6601.711653048776`\), \(-6.828926338636988`\), \ \(-2.485525919260343`\), 0}, {\(-2552.087806256307`\), 3063.731752628119`, \(-6498.362878662377`\), \(-6.504110930007682`\), \ \(-3.0453692304831437`\), 1.1036178122174236`}, {\(-3720.7039830179037`\), 2449.94260479205`, \(-6191.646366829996`\), \(-5.976176391284324`\), \ \(-3.509446535489579`\), 2.1716943309224956`}, {\(-4773.052367542169`\), 1759.5822158981725`, \(-5691.426252222087`\), \ \(-5.262637690294909`\), \(-3.863116460731575`\), 3.170023522894014`}, {\(-5676.584749564745`\), 1014.3458239304405`, \(-5013.732165366227`\), \ \(-4.387143373618118`\), \(-4.095564731824589`\), 4.066994814616169`}, {\(-6403.66905093455`\), 237.5638039951957`, \(-4180.167376692139`\), \ \(-3.3785047642245365`\), \(-4.200137483684298`\), 4.834710686003183`}, {\(-6932.448100639624`\), \ \(-546.5919638189616`\), \(-3217.1258196589097`\), \(-2.2695408729056097`\), \ \(-4.174486882572278`\), 5.4499079543213265`}, {\(-7247.471055675766`\), \ \(-1313.9116417854955`\), \(-2154.858292188162`\), \(-1.0958114624974253`\), \ \(-4.020533517343553`\), 5.894647090380231`}, {\(-7340.083447666443`\), \ \(-2040.928067216727`\), \(-1026.4329647003283`\), 0.10568867235577295`, \(-3.7442623201353387`\), 6.156753563784076`}, {\(-7208.574991673947`\), \ \(-2705.645685863801`\), 133.36346711912387`, 1.2978074303455387`, \(-3.3553779783778337`\), 6.230013912721875`}, {\(-6858.09575929358`\), \ \(-3288.1944190199465`\), 1289.1396721307258`, 2.4444298380589378`, \(-2.8668512363231127`\), 6.114144885012337`}, {\(-6300.360321176668`\), \(-3771.3911718189124`\ \), 2405.9888569081354`, 3.511548403661685`, \(-2.294389125646091`\), 5.814565261411292`}, {\(-5553.165649726921`\), \ \(-4141.197752285639`\), 3450.5196653037597`, 4.4681839102471645`, \(-1.6558605076155248`\), 5.342006345254369`}, {\(-4639.751974737261`\), \ \(-4387.069478898511`\), 4391.802197234311`, 5.287145386985574`, \(-0.9707041937421442`\), 4.711998826439215`}, {\(-3588.0367444255953`\), \(-4502.193340555801`\ \), 5202.2094353593575`, 5.945629751209067`, \(-0.2593413202041743`\), 3.9442715904973222`}, {\(-2429.7509029933362`\), \ \(-4483.618051725812`\), 5858.140588189331`, 6.425670134471138`, 0.4573924637446428`, 3.062093137763323`}, {\(-1199.5044864346926`\), \(-4332.280697439013`\ \), 6340.617994904145`, 6.714446966074374`, 1.1587852678142594`, 2.0915797790201753`}, {66.19430657078055`, \(-4052.935999155976`\), 6635.753097547799`, 6.804477694374749`, 1.8247529253878103`, 1.060987775067187`}, {1329.9453273218628`, \(-3653.9947559260054`\), 6735.07956523158`, 6.6937000745055615`, 2.436307584225485`, 0}}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(TableForm[%]\)], "Input"], Cell[BoxData[ TagBox[GridBox[{ {\(-1303.609776285786`\), "3581.63842412549`", \(-6601.711653048776`\), \ \(-6.828926338636988`\), \(-2.485525919260343`\), "0"}, {\(-2552.087806256307`\), "3063.731752628119`", \(-6498.362878662377`\), \ \(-6.504110930007682`\), \(-3.0453692304831437`\), "1.1036178122174236`"}, {\(-3720.7039830179037`\), "2449.94260479205`", \(-6191.646366829996`\), \ \(-5.976176391284324`\), \(-3.509446535489579`\), "2.1716943309224956`"}, {\(-4773.052367542169`\), "1759.5822158981725`", \(-5691.426252222087`\), \ \(-5.262637690294909`\), \(-3.863116460731575`\), "3.170023522894014`"}, {\(-5676.584749564745`\), "1014.3458239304405`", \(-5013.732165366227`\), \ \(-4.387143373618118`\), \(-4.095564731824589`\), "4.066994814616169`"}, {\(-6403.66905093455`\), "237.5638039951957`", \(-4180.167376692139`\), \ \(-3.3785047642245365`\), \(-4.200137483684298`\), "4.834710686003183`"}, {\(-6932.448100639624`\), \(-546.5919638189616`\), \ \(-3217.1258196589097`\), \(-2.2695408729056097`\), \(-4.174486882572278`\), "5.4499079543213265`"}, {\(-7247.471055675766`\), \(-1313.9116417854955`\), \ \(-2154.858292188162`\), \(-1.0958114624974253`\), \(-4.020533517343553`\), "5.894647090380231`"}, {\(-7340.083447666443`\), \(-2040.928067216727`\), \ \(-1026.4329647003283`\), "0.10568867235577295`", \(-3.7442623201353387`\), "6.156753563784076`"}, {\(-7208.574991673947`\), \(-2705.645685863801`\), "133.36346711912387`", "1.2978074303455387`", \(-3.3553779783778337`\), "6.230013912721875`"}, {\(-6858.09575929358`\), \(-3288.1944190199465`\), "1289.1396721307258`", "2.4444298380589378`", \(-2.8668512363231127`\), "6.114144885012337`"}, {\(-6300.360321176668`\), \(-3771.3911718189124`\), "2405.9888569081354`", "3.511548403661685`", \(-2.294389125646091`\), "5.814565261411292`"}, {\(-5553.165649726921`\), \(-4141.197752285639`\), "3450.5196653037597`", "4.4681839102471645`", \(-1.6558605076155248`\), "5.342006345254369`"}, {\(-4639.751974737261`\), \(-4387.069478898511`\), "4391.802197234311`", "5.287145386985574`", \(-0.9707041937421442`\), "4.711998826439215`"}, {\(-3588.0367444255953`\), \(-4502.193340555801`\), "5202.2094353593575`", "5.945629751209067`", \(-0.2593413202041743`\), "3.9442715904973222`"}, {\(-2429.7509029933362`\), \(-4483.618051725812`\), "5858.140588189331`", "6.425670134471138`", "0.4573924637446428`", "3.062093137763323`"}, {\(-1199.5044864346926`\), \(-4332.280697439013`\), "6340.617994904145`", "6.714446966074374`", "1.1587852678142594`", "2.0915797790201753`"}, {"66.19430657078055`", \(-4052.935999155976`\), "6635.753097547799`", "6.804477694374749`", "1.8247529253878103`", "1.060987775067187`"}, {"1329.9453273218628`", \(-3653.9947559260054`\), "6735.07956523158`", "6.6937000745055615`", "2.436307584225485`", "0"} }, RowSpacings->1, ColumnSpacings->3, RowAlignments->Baseline, ColumnAlignments->{Left}], Function[ BoxForm`e$, TableForm[ BoxForm`e$]]]], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Convert Cartesian Elements to Kepler Element", "Section 1"], Cell[BoxData[ RowBox[{\(Cartesian2Kepler[stateVector_, \ debug_: False]\), ":=", " ", RowBox[{"Module", "[", RowBox[{\({\[IndentingNewLine]x, y, z, vx, vy, vz, r, v, R, V, \[IndentingNewLine]a, e, i, \[CapitalOmega], \[Omega], M, \[Theta], \ EA, \ eccVector, \[IndentingNewLine]\[Mu] = 398600.436, \ p, \ hVector, \ h, \ n, argLat\[IndentingNewLine]}\), ",", "\[IndentingNewLine]", RowBox[{\(If[Length[stateVector] \[NotEqual] \ 6, \ Print["\"]; \ \[IndentingNewLine]Abort[];\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", \({x, y, z, vx, vy, vz} = stateVector\), ";", "\[IndentingNewLine]", \(R = {x, y, z}\), ";", "\[IndentingNewLine]", \(V = {vx, vy, vz}\), ";", "\[IndentingNewLine]", \(r\ = \ \@\(R . R\)\), ";", "\[IndentingNewLine]", \(v = \@\(V . V\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ calculate\ the\ eccentricity\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight-> "Plain"], \(eccVector\ = \ \(\(V . V - \[Mu]/r\)\/\[Mu]\) R - \(R . V\/\[Mu]\) V\), ";", "\[IndentingNewLine]", \(e = \@\(eccVector . eccVector\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ normalize\ the\ eccentricity\ vector\ *) \), FontWeight->"Plain"], "\[IndentingNewLine]", \(eccVector\ = \ \(1\/e\) eccVector\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ angular\ momentum\ per\ unit\ mass\ *) \), FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], \(hVector\ = \ Cross[R, \ V]\), ";", "\[IndentingNewLine]", \(h\ = \ \@\(hVector . hVector\)\), ";", "\[IndentingNewLine]", \(hVector\ = \ \(1\/h\) hVector\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ semimajor\ axis\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], "\[IndentingNewLine]", \(a\ = \ h\^2\/\(\[Mu] \((1 - e\^2)\)\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ inclination\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(i\ = \ ArcCos[hVector[\([3]\)]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ right\ ascension\ of\ ascending\ node\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(n\ = \ Cross[{0, 0, 1}, hVector]\), ";", "\[IndentingNewLine]", \(n\ = \ \(1\/\@\(n . n\)\) n\), ";", "\[IndentingNewLine]", \(\[CapitalOmega] = ArcCos[n[\([1]\)]]\), ";", "\[IndentingNewLine]", \(If[ n[\([2]\)] < 0, \ \[CapitalOmega] = 2*Pi - \[CapitalOmega]]\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ argument\ of\ perigee\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(\[Omega]\ = \ ArcCos[n . eccVector]\), ";", "\[IndentingNewLine]", \(If[ eccVector[\([3]\)] < 0, \ \[Omega] = 2*Pi\ - \ \[Omega]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ argument\ of\ Latitude\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(argLat\ = \ ArcCos[n . \((\(1\/r\) R)\)]\), ";", "\[IndentingNewLine]", \(If[z < 0, \ argLat\ = \ 2*Pi\ - \ argLat]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*compute\ the\ true\ anomaly\ *) \), FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(\[Theta] = argLat\ - \ \[Omega]\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ eccentric\ anomaly\ *) \), FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(EA\ = \ ArcTan[\(e + Cos[\[Theta]]\)\/\(1 + e*Cos[\[Theta]]\), \ \(\(\@\ \(1 - e\^2\)\) Sin[\[Theta]]\)\/\(1 + e*Cos[\[Theta]]\)]\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ compute\ the\ mean\ anomaly\ using\ Kepler' s\ Equation\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(M\ = \ EA\ - \ e*Sin[EA]\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", \(If[ debug, \[IndentingNewLine]\(Print[\[IndentingNewLine]"\", R, "\< r=\>", r, \[IndentingNewLine]"\<\nV=\>", V, "\< v=\>", v, \[IndentingNewLine]"\<\neccVector=\>", eccVector, "\< e=\>", e, \[IndentingNewLine]"\<\nhVector=\>", hVector, "\< h=\>", h, \[IndentingNewLine]"\<\na=\>", a, "\<\ni=\>", i, "\< rad = \>", \ i*180.0/Pi, "\< deg\>", \[IndentingNewLine]"\<\nn=\>", n, "\< magnitude = \>", \@\(n . n\), \ \[IndentingNewLine]"\ \<\n\[CapitalOmega]=\>", \[CapitalOmega], "\< rad=\>", \[CapitalOmega]*180.0/ Pi, "\< deg\>", \[IndentingNewLine]"\<\n\[Omega]=\>", \ \[Omega], "\< rad=\>", \ \[Omega]*180.0/ Pi, "\< deg\>", \[IndentingNewLine]"\<\n\[Theta]=\>", \ \[Theta], "\< rad=\>", \[Theta]*180.0/ Pi, "\< deg\>", \[IndentingNewLine]"\<\nE=\>", EA, "\< rad=\>", EA*180.0/Pi, "\< deg\>", \[IndentingNewLine]"\<\nM=\>", M, "\< rad=\>", M*180.0/Pi, "\< deg\>"];\)\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(Return[{a, e, i, \[CapitalOmega], \[Omega], M}]\), ";"}]}], "\[IndentingNewLine]", "\[IndentingNewLine]", "]"}]}]], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(Cartesian2Kepler[{7235.312243818173`, \(-1278.2214242749383`\), 7235.312243818173`, \(-4.315205778480563`\), \(-6.78421570401441`\), \ \(-6.78421570401441`\)}]\)], "Input", CellOpen->False], Cell[BoxData[ \({\(-11947.510861638562`\), 1.547105966594061`, 2.3053472670686834`, 1.8707571899222442`, 3.0648791199030234`, \(\(0.`\)\(\[InvisibleSpace]\)\) + 0.4093405000464323`\ \[ImaginaryI]}\)], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Calculate Right Ascension of Grenwich", "Section 1"], Cell[BoxData[ RowBox[{ RowBox[{\(rag[date_, debug_: False]\), ":=", " ", RowBox[{"Module", "[", "\[IndentingNewLine]", StyleBox[\( (*\ date\ is\ days\ since\ \(1/1\)/2000\ at\ 00 : 00\ GMT\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", RowBox[{\({\ TU, \[Alpha], \[Alpha]hrs, minutes, seconds}\), ",", "\[IndentingNewLine]", RowBox[{\(TU\ = \ date/36525.0\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ Compute\ angle\ in\ seconds\ of\ time\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], "\[IndentingNewLine]", StyleBox["\[IndentingNewLine]", FontWeight-> "Plain"], \(\[Alpha] = \((24110.54841 + 8640184.812866*TU + 0.093104*TU\^2 - 6.2*10\^\(-6\)*TU\^3)\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ convert\ seconds\ of\ time\ to\ hours\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], \(\[Alpha]hrs\ = \ \[Alpha]/3600.0\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ convert\ hours\ to\ degrees\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(\[Alpha] = \[Alpha]hrs*15\), ";", "\[IndentingNewLine]", \(\[Alpha] = Mod[\[Alpha], 360]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(If[ debug, \ \[IndentingNewLine]minutes\ = \ 60*\((\[Alpha]hrs - Floor[\[Alpha]hrs]\ )\); \[IndentingNewLine]seconds\ \ = 60*\((minutes - Floor[minutes])\); \[IndentingNewLine]Print["\", date, "\< \[Alpha]=\>", \[Alpha], "\< deg =\>", \[Alpha]* Pi/180, "\< rads = \>", Floor[\[Alpha]hrs], "\<:\>", Floor[minutes], "\<:\>", seconds\[IndentingNewLine]];\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ return\ angle\ in\ radians\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(Return[\[Alpha]*Pi/180.0]\), ";"}]}], "\[IndentingNewLine]", "]"}]}], ";"}]], "Input"] }, Open ]], Cell[CellGroupData[{ Cell["Calculate Latitude and Longitude from Kepler Vector", "Section 1"], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{\(latLong[kepler_, \ time_, \ debug_: False]\), ":=", " ", "\[IndentingNewLine]", StyleBox[\( (*\ time\ is\ in\ days\ since\ \(1/1\)/2000\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], RowBox[{"Module", "[", RowBox[{\({a, \ e, \ i, \ \[CapitalOmega], \[Omega], \ M, cartesian, \ x, y, z, vx, vy, vz, R, V, r, v, \[Alpha]g, longitude, \[Delta], \[Phi], \[Phi]prime, \[Delta]prime, \ height, \ new\[Phi]prime, \ f = 1/298.3, \[Epsilon] = 10\^\(-10\), error, counter = 0, rc}\), ",", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{\({a, \ e, \ i, \ \[CapitalOmega], \[Omega], \ M} = kepler\), ";", "\[IndentingNewLine]", \(cartesian\ = \ kepler2Cartesian[kepler]\), ";", "\[IndentingNewLine]", \({x, y, z, vx, vy, vz} = cartesian\), ";", "\[IndentingNewLine]", \(R = {x, y, z}\), ";", "\[IndentingNewLine]", \(r = \@\(R . R\)\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ longitude\ is\ easy\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(\[Alpha]g = rag[time]\), ";", "\[IndentingNewLine]", \(longitude = ArcTan[x, y] - \[Alpha]g\), ";", " ", "\[IndentingNewLine]", "\[IndentingNewLine]", StyleBox[\( (*\ declination\ *) \), FontWeight->"Plain"], StyleBox[" ", FontWeight->"Plain"], StyleBox["\[IndentingNewLine]", FontWeight->"Plain"], "\[IndentingNewLine]", \(\[Delta] = ArcTan[\@\(x\^2 + y\^2\), z]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \( (*\ iterate\ for\ geodetic\ latitude\ *) \), " ", "\[IndentingNewLine]", \(new\[Phi]prime = \[Delta]\), ";", "\[IndentingNewLine]", \(error\ = \ 10\), ";", "\[IndentingNewLine]", \(While[ error > \[Epsilon]\ \[And] \ \(counter++\) < 15, \[IndentingNewLine]\[Phi]prime = new\[Phi]prime; \[IndentingNewLine]rc = 6378.136*\@\(\(1 - 2*\ f + f\^2\)\/\(1 - \((2*f - f\^2)\) \ \((Cos[\[Phi]prime])\)\^2\)\); \[IndentingNewLine]\[Phi] = ArcTan[Tan[\[Phi]prime]\/\((1 - f)\)\^2]; \ \[IndentingNewLine]height\ = \ \@\(r\^2 - \(rc\^2\) \((Sin[\[Phi] - \ \[Phi]prime])\)\^2\) - rc*Cos[\[Phi] - \[Phi]prime]; \ \[IndentingNewLine]\[Delta]prime = ArcSin[\(height\/r\) Sin[\[Phi] - \[Phi]prime]]; \ \[IndentingNewLine]new\[Phi]prime = \[Delta] - \[Delta]prime; \ \[IndentingNewLine]error\ = \ Abs[new\[Phi]prime - \[Phi]prime]; \[IndentingNewLine]\ \[IndentingNewLine]If[debug, Print[\[IndentingNewLine]"\<\[Phi]=\>", \[Phi], "\< \ \[Phi]'=\>", \[Phi]prime, "\< new \[Phi]'=\>", new\[Phi]prime, "\< \[Delta]=\>", \[Delta]prime, "\< \ error=\>", \ error\[IndentingNewLine]]];\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", "\[IndentingNewLine]", "\[IndentingNewLine]", \(If[ debug, \[IndentingNewLine]\(Print[\[IndentingNewLine]"\<\ cartesian state=\>", \ cartesian, \[IndentingNewLine]"\<\nr=\>", r, "\<\n\[Alpha]g=\>", \[Alpha]g, \ "\< rad = \>", \ \ \[Alpha]g*180.0/Pi, "\< deg\>", \[IndentingNewLine]"\<\nlongitude=\>", longitude, \ "\< rad = \>", \ longitude*180.0/ Pi, "\< deg\>", \[IndentingNewLine]"\<\n\[Delta]=\>", \ \[Delta], \ "\< rad = \>", \ \[Delta]*180.0/ Pi, "\< deg\>", \[IndentingNewLine]"\<\nlatitude: \ \>", \[Phi], "\< rad = \>", \ \[Phi]*180.0/ Pi, "\< \ deg\>"\[IndentingNewLine]\[IndentingNewLine]];\)\[IndentingNewLine]]\), ";", "\[IndentingNewLine]", "\[IndentingNewLine]", \(Return[{\[Phi], \ longitude}]\)}]}], "\[IndentingNewLine]", "\[IndentingNewLine]", "]"}]}], ";"}]], "Input"], Cell[BoxData[ \(General::"spell1" \(\(:\)\(\ \)\) "Possible spelling error: new symbol name \"\!\(\[Alpha]g\)\" is \ similar to existing symbol \"\!\(\[Alpha]\)\"."\)], "Message"], Cell[BoxData[ \(General::"spell1" \(\(:\)\(\ \)\) "Possible spelling error: new symbol name \"\!\(\[Delta]prime\)\" is \ similar to existing symbol \"\!\(\[Phi]prime\)\"."\)], "Message"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(latLong[{7700, 0.01, Pi/3, \ \((200/180)\)*Pi, \ 3*Pi/2, \ 10*Pi/180}, 5]*180.0/Pi\)], "Input"], Cell[BoxData[ \({\(-58.609494617593455`\), 24.405434526428408`}\)], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Examples", "Section 1"], Cell[CellGroupData[{ Cell[BoxData[{ \(\(a = 7000;\)\), "\[IndentingNewLine]", \(\(e = 0.05;\)\), "\[IndentingNewLine]", \(\(i = 83.0*Pi/180.0;\)\), "\[IndentingNewLine]", \(\(\[CapitalOmega] = 120.0*Pi/180.0;\)\), "\[IndentingNewLine]", \(\(\[Omega] = 17.0*Pi/180.0;\)\), "\[IndentingNewLine]", \(\(M = 43.0*Pi/180.0;\)\[IndentingNewLine]\), "\[IndentingNewLine]", \(elements\ = \ {a, \ e, \ i, \ \[CapitalOmega], \ \[Omega], \ M}\)}], "Input"], Cell[BoxData[ \({7000, 0.05`, 1.4486232791552935`, 2.0943951023931953`, 0.29670597283903605`, 0.7504915783575618`}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell["convert kepler elements to cartesian", "Subsection"], Cell[CellGroupData[{ Cell[BoxData[ \(cartesianSet\ = \ kepler2Cartesian[elements]\)], "Input"], Cell[BoxData[ \({\(-2116.3168830559275`\), 2185.122037555789`, 6028.633884937644`, 3.0667088175921258`, \(-6.204415405066861`\), 3.635310149798732`}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \({x, y, z, xvelocity, \ yvelocity, \ zvelocity} = cartesianSet\)], "Input"], Cell[BoxData[ \({\(-2116.3168830559275`\), 2185.122037555789`, 6028.633884937644`, 3.0667088175921258`, \(-6.204415405066861`\), 3.635310149798732`}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(x\)], "Input"], Cell[BoxData[ \(\(-2116.3168830559275`\)\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(y\)], "Input"], Cell[BoxData[ \(2185.122037555789`\)], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["convert back to kepler and see if it really works", "Subsection"], Cell[CellGroupData[{ Cell[BoxData[ \(newElements\ = \ Cartesian2Kepler[cartesianSet]\)], "Input"], Cell[BoxData[ \({7000.000000000004`, 0.050000000000000315`, 1.4486232791552935`, 2.0943951023931953`, 0.2967059728390388`, 0.7504915783575589`}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(newElements - elements\)], "Input"], Cell[BoxData[ \({3.637978807091713`*^-12, 3.122502256758253`*^-16, 0.`, 0.`, 2.7755575615628914`*^-15, \(-2.886579864025407`*^-15\)}\)], "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Calculate the latititude & longitude of this satellite now", \ "Subsection"], Cell[CellGroupData[{ Cell[BoxData[ \(dateNow = Date[]\)], "Input"], Cell[BoxData[ \({2001, 6, 18, 7, 45, 53}\)], "Output"] }, Open ]], Cell["\<\ To use a specific date enter dateNow = {year, month, day, hours, \ minutes, seconds} The next line gives the seconds since 1/1/1900\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(secsNowSince1900 = FromDate[dateNow]\)], "Input"], Cell[BoxData[ \(3201839153\)], "Output"] }, Open ]], Cell["\<\ Calculate the seconds since 1/1/2000, and then calculate the days \ since 1/1/2000\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(secsSinceJan12000 = secsNowSince1900 - FromDate[{2000, 1, 1, 0, 0, 0}]\)], "Input"], Cell[BoxData[ \(46165553\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(daysSinceJan12000 = secsSinceJan12000/86400.0\)], "Input"], Cell[BoxData[ \(534.3235300925926`\)], "Output"] }, Open ]], Cell["\<\ Calculate the latitude and longitude now \ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \({latitude, \ longitude} = \((180.0/Pi)\)* latLong[elements, \ daysSinceJan12000]\)], "Input"], Cell[BoxData[ \({63.37064468018199`, \(-133.0316149485151`\)}\)], "Output"] }, Open ]] }, Open ]] }, Open ]] }, FrontEndVersion->"4.1 for X", ScreenRectangle->{{0, 1024}, {0, 768}}, WindowSize->{593, 599}, WindowMargins->{{Automatic, 203}, {Automatic, 38}}, StyleDefinitions -> "Report.nb" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1727, 52, 62, 0, 70, "Section 1"], Cell[1792, 54, 1052, 18, 362, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[2881, 77, 57, 0, 70, "Section 1"], Cell[2941, 79, 499, 8, 163, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[3477, 92, 67, 0, 70, "Section 1"], Cell[CellGroupData[{ Cell[3569, 96, 6258, 124, 748, "Input"], Cell[9830, 222, 181, 3, 45, "Message"] }, Open ]], Cell[CellGroupData[{ Cell[10048, 230, 153, 3, 64, "Input"], Cell[10204, 235, 2973, 50, 335, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[13214, 290, 45, 1, 32, "Input"], Cell[13262, 293, 3637, 71, 385, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[16948, 370, 65, 0, 70, "Section 1"], Cell[17016, 372, 7393, 138, 1571, "Input"], Cell[CellGroupData[{ Cell[24434, 514, 220, 4, 19, "Input", CellOpen->False], Cell[24657, 520, 237, 4, 47, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[24943, 530, 58, 0, 70, "Section 1"], Cell[25004, 532, 3050, 61, 543, "Input"] }, Open ]], Cell[CellGroupData[{ Cell[28091, 598, 72, 0, 70, "Section 1"], Cell[CellGroupData[{ Cell[28188, 602, 4635, 87, 1062, "Input"], Cell[32826, 691, 190, 3, 45, "Message"], Cell[33019, 696, 197, 3, 45, "Message"] }, Open ]], Cell[CellGroupData[{ Cell[33253, 704, 124, 2, 48, "Input"], Cell[33380, 708, 81, 1, 47, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[33510, 715, 29, 0, 70, "Section 1"], Cell[CellGroupData[{ Cell[33564, 719, 467, 8, 144, "Input"], Cell[34034, 729, 139, 2, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[34210, 736, 58, 0, 39, "Subsection"], Cell[CellGroupData[{ Cell[34293, 740, 78, 1, 32, "Input"], Cell[34374, 743, 180, 3, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[34591, 751, 101, 2, 32, "Input"], Cell[34695, 755, 180, 3, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[34912, 763, 34, 1, 32, "Input"], Cell[34949, 766, 58, 1, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[35044, 772, 34, 1, 32, "Input"], Cell[35081, 775, 52, 1, 47, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[35182, 782, 71, 0, 39, "Subsection"], Cell[CellGroupData[{ Cell[35278, 786, 81, 1, 32, "Input"], Cell[35362, 789, 175, 3, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[35574, 797, 55, 1, 32, "Input"], Cell[35632, 800, 157, 2, 67, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[35838, 808, 82, 1, 39, "Subsection"], Cell[CellGroupData[{ Cell[35945, 813, 49, 1, 32, "Input"], Cell[35997, 816, 58, 1, 47, "Output"] }, Open ]], Cell[36070, 820, 156, 5, 60, "Text"], Cell[CellGroupData[{ Cell[36251, 829, 69, 1, 32, "Input"], Cell[36323, 832, 44, 1, 47, "Output"] }, Open ]], Cell[36382, 836, 106, 3, 28, "Text"], Cell[CellGroupData[{ Cell[36513, 843, 110, 2, 48, "Input"], Cell[36626, 847, 42, 1, 47, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[36705, 853, 78, 1, 32, "Input"], Cell[36786, 856, 52, 1, 47, "Output"] }, Open ]], Cell[36853, 860, 66, 4, 60, "Text"], Cell[CellGroupData[{ Cell[36944, 868, 121, 2, 48, "Input"], Cell[37068, 872, 79, 1, 47, "Output"] }, Open ]] }, Open ]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)