Maple solves satellite Section F_7

This is the Title PhL 1.11.15 F.7 Numerical Solution of Satellite Equations of Motion The equations of motion are stated in (F.5.8) which we replicate here. 1 θ + sinθcosθ(3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 // 2 φ + 2 cotθ – ωcosφ (ωsinφ -2) = 0 // (F.5.7) 3 = vθ 4 = vφ (F.5.8) To keep our Maple code compact, we make these name changes: θ ↔ T Theta φ ↔ P Phi vθ ↔ vT vTheta vφ ↔ vP vPhi ω ↔ w The first step is to enter the four equations and other setup data including ω = 1 : For each solution run we must select some initial conditions. Notice that we have set ω = 1 T = 2π = 6.28 seconds Verification of in-plane libration The initial conditions are taken to be θ = 0.2 = 0 φ = π/2 = 0 . Here is a plot of the result where θ(t) is magnified by a factor of 10: The azimuth φ(t) holds at π/2 while θ(t) oscillates. From (F.5.10) for in-plane libration one predicts, Tocs1 = T/ = 6.28/1.73 = 3.63 sec which is verified by the vertical line in the figure. Verification of out-of-plane libration The initial conditions are taken to be θ = 0.2 = 0 φ = 0 = 0 . Here is a plot of the result where again θ(t) is magnified by a factor of 10: Maple has chosen to show the oscillation by keeping θ(t) positive and having φ(t) jump by π at each zero crossing. This is just another way to represent a cosine function. As suggested below (F.5.13), one cannot maintain φ(t) = 0 for the out-of-plane case and one sees φ(t) creeping up from 0 around t = 0.5 sec. From (F.5.13) for out-of-plane libration one predicts, Tocs2 = T/2 = 6.28/2 = 3.14 which is verified by the vertical line in the figure. Maple can also plot (t) and (t) and here is a plot showing all four curves θ,φ,, = red,black,green,blue: This plot hints at a technical problem with the ODE system of equations as entered above. The second equation contains a factor cot(θ) which blows up at θ = 0, φ + 2 cotθ – ωcosφ (ωsinφ -2) = 0, so whenever the red θ(t) curve approaches 0, one sees φ getting very large (the slope of the blue curve = vφ). Unless = 0, the solution cannot pass through the θ = 0 value. The problem is that θ = 0 is a singular point in spherical coordinates where azimuth φ is undefined. One solution to this problem might be to convert the differential equations from spherical coordinates θ and φ to Cartesian coordinates x,y,z where x = rsinθcosφ y = rsinθsinφ z = rcosθ (E.2.5) A more general solution so θ is a small angle. Here then is the solution run and a plot ofθ(t) in red and φ(t) in black (θ is magnified by a factor of 10), One can see that the period of the small-θ oscillation is about 3.7 seconds so ωosc = 2π/3.7 = 1.7. This is in agreement with (F.5.10) which says ωosc1 = ω = 1.7 * 1 = 1,7 for in-plane libration. If θ = 1 radian, so the polar angle is no longer small, the period gets longer : as happens in a spherical pendulum. Here we plot 30*θ(t) and φ(t) on the same graph (we magnify θ(t) in order to see both plots), The motion of this example is quite complicated: the polar angle θ(t) of the dumbbell oscillates back and forth without reaching 0 while the dumbbell rotates around the vertical axis in azimuth φ since we started it with = -1, so the φ just winds around in the negative direction. Another way to plot this motion is as a trajectory in θ-φ space, In this graph, mass m1 starts at the top right and then moves along the curve. Each descent of π takes the satellite halfway around the vertical axis. The azimuthal velocity is quite uneven as energy shifts between the two degrees of freedom. Here we add the velocity in green and in blue to our original plot ** to get As the red curve gets close to θ = 0, the cotθ in equation #2 gets large, forcing φ to be large as indicated by the large negative blue peaks. The simulation is very non-robust and dies if the point θ = 0 is hit since then cotθ = ∞. This is unfortunate, since we would like to plot the in-plane and out-of-plane oscillations which should pass through θ = 0. The problem of course is that θ = 0 is a singular point of the spherical coordinate system where the value of φ is undefined, so our simple simulation not surprisingly has problems at θ = 0. One could fix this problem by changing from θ,φ to x,y.z coordinates, but we leave that task to the reader. x = r1sinθcosφ y = r1sinθsinφ z = r1cosθ = r1cosθcosφ - r1sinθ sinφ = r1cosθsinφ + r1sinθ cosφ = -r1sinθ The motion is a bit difficult to visualize. From t = 0 to t = 1 θ(t) drops down from its initial 0.2 to some small θmin while φ moves from 0 to -π, so the dumbbell has gone half way around the vertical axis. The lower dumbbell mass moves somewhat like a spherical pendulum (while the upper mass as an upside down spherical pendulum). The angle θ(t) shows cosine-like behavior between some θmax and θmin, while the azimuthal coordinate φ just winds around the vertical axis roughly once per oscillation of θ. One can see that the period of the oscillation is about 1.6 seconds so ωosc ≈ 2π/1.6 ≈ 4. For this example the dumbbell initial conditions are: θ = 0.2 = 0 φ = 0 = -1 The dumbbell oscillates in θ while it rotates in azimuth in the - direction. Next we group the four equations into set "eqs", and put initial conditions into "init", and set ω = 1: The call to solve the system of equations is then, We then plot T = θ and P = φ for t =