Need help on my labs.planet.dvi Physics 336/536 – Computational Physics – Fall 2022 Code plan for...

Question

Need help on my labs.planet.dvi Physics 336/536 – Computational Physics – Fall 2022 Code plan for planet orbiting fixed Sun in circular orbit: • Goal: Calculate circular orbit of planet around fixed Sun and create 3D visualization. Use to test ideas of planetary motion: Kepler’s laws. • Equations modeled: ∆~v = ~a∆t ∆~r = ~v∆t ~a = ~Fg M⊕ = − GM⊙ r3 ~r • Strategy: – Use astronomic units: GM⊙ = 4π 2 – Use initial conditions to guarantee circular orbit. – Use Euler-Cromer or RK algorithm with implicit Cromer correc- tion – Measure period by sign change in one position coordinate. – Borrow planet visualization from visual Demo, planet.py – Plot 2D plots as needed • Pseudocode outline: 1. Set values for givens: – Orbital radius R in AU (1 for Earth) – Orbit period T in years – Circular orbit speed √ GM⊙/R⊕ – OR circular speed = 2πR/T (2π for Earth) – Time step of simulation ∆t ∼ T/100 2. Set initial conditions, for circular orbit. – position (x, y, z) = (R, 0, 0) – velocity (vx, vy, vz) = (0, v, 0) 1 3. Start 3D visual: – x, y, z axes drawn – Scene viewed on (111) axis – create sphere for Sun, for planet – create curve to draw trail for planet. 4. Do Euler-Cromer calculation – Calculate period on the fly 5. Finish up • Pseudo-code for Euler-Cromer calculation INPUT: initial position, initial velocity, ∆t, planet and trail. while (1): Calculate acceleration on planet. Update velocity update position 2nd (Cromer correction) tnew = told +∆t move planet to new position every 10 moves, add point to trail. repeat RETURN: latest value of position and velocity (allows restart) • Checks: 1. Look at quality of orbit as function of ∆t to obtain a reasonable value of ∆t. Do this by means suggested in problem 4.1 of the textbook. • Questions: – Modify the planetary program to check Kepler’s 2nd law – that the orbit sweeps out equal areas in equal times. This can be done geometrically or by noting that this is equivalent to conservation of angular momentum. (Why?) Try for an elliptical orbit of the sort generated in your check exercise above. P.S. – recall basic definition of angular momentum: ~L = ~r × ~p and remember that vpython has a cross product function. 2 – Check numerically Kepler’s 3rd law (T 2yr = R 3 AU) for circular or- bits. You may wish to use the radii of the planets (most are close to circular in orbit) and compare your calculated orbits to pub- lished values, say on wikipedia. Or pick several arbitrary radii and calculate the period and obtain a relationship graphically. Create a plot showing your results in an effective way. (Hint – it helps here to use v = √ GM⊙/R, obtained by setting F = mv 2/r and solving for v. This allows you to choose R and v without knowing T beforehand.) 3 • Followup for next week: Looking ahead to next week – you will have one of two options to add onto this lab. 1. Repeat the above lab, putting in the proper boundary conditions (semi-major axis a and eccentricity e) for at least the five planets of table 4.2. How are your results improved? One can also show the first 5 planets in a solar system demo. An issue here would be picking a time step that will make Mercury accurate and Jupiter not take forever to move. 2. We will discuss the 3-body (section 4.3) in class next week, and look at the effect of Jupiter on Mars. This is an example of the first significant computational problem in physics. One technique is to plot two trajectories for Mars, one being Mars in the Sun’s gravitational field only, and the second being the solution using the three-body problem result, starting from the same initial con- ditions. Also the trick of making Jupiter even more gigantic or making the Sun finite in mass may be explored. • Graduate level problem: 1. Any problem from sections 4.2 through 4.5. 4    SHOassignment20.dvi Physics 336/536 – Computational Physics Spring 2022 SHO and timesteps: Read the code description for the simulation of a , attached. Please take the partially completed code SHO simple.py, which implements it. After finishing the code by filling in question marks, and answering questions, your goal is to explore the consequence of choice of integrating scheme and timestep size on the dynamics of a simple oscillator. 1. Pick two or more time steps, both much less than one period of motion. 2. Set up the Euler equation for the SHO, following the code plan, page 5. Make sure to only use original values of force, θ and ω. So, update accelaration, then θ, then ω. 3. Plot the energy versus time: E(t) = 1 2 ml2ω2(t) + 1 2 mglθ2(t) = 1 2 mgl ( θ2 + (ωT/2π)2 ) for each time step where T is the period. Do so for at least ten periods of SHO motion. Also make a phase-space plot, plotting ω(t)T/2π versus θ(t). 4. Repeat this for Euler-Cromer as described in the code plan. This uses information from the “future” (already updated data) to help deter- mine the update of current values of other data. For example, calcu- late acceleration, then update velocity, and use new value of velocity to update angle. BONUS Make a function that calculates the motion of a spring given initial conditions, as done in class for the trajectory of a projectile. Do two functions, one for Euler method, one for Euler-Cromer. 5. Analysis and discussion – answer the following: • What curve should the phase-space plot produce if energy is con- served? Does it do this? 1 • What is the long-time behavior of E(t) for each method, and which is better? How do energies compare to the theoretical energy? • What is the effect of using a smaller time step – can you make problems of a given method go away or at least get better? 6. For a report, make it short. I will just take plots with good captions and labels and two paragraphs – one discussing the input parameters you picked and why, and one of discussion of your results. E.g., a plot comparing several time-steps for one method, and repeating for the second. 7. GRADUATE: redo this experiment using the leapfrog method – used in Astronomy. It obeys time reversal symmetry, meaning you can use it backwards in time and forwards in time and get the same trajectory. (Not possible with Euler!). The rule is: vi+1/2 = vi + ai ∆t 2 , xi+1 = xi + vi+1/2∆t, vi+1 = vi+1/2 + ai+1 ∆t 2 . Here, vi+1/2 is a temporary throw-away point which is leap-frogged over in time when updating x. You only keep xi+1 and vi+1. An equivalent rule that hides the leapfrogging is xi+1 = xi + vi ∆t+ 1 2 ai ∆t 2, vi+1 = vi + 1 2 (ai + ai+1)∆t. For more information, see Leapfrog integration in Wikipedia and Hut and Makino, The Art of Computational Science, Chapter 18. Web link is http://www.artcompsci.org/kali/pub/msa/title.html. BONUS Repeat for Runge-Kutta method (second order), described in appendix A2 of textbook and in the E-textbook by Landau et al. This does not conserve energy but is much more accurate than Euler’s method. 2 Physics 336/536 – Computational Physics – Spring 2022 Code plan for small angle pendulum: • Goal: Implement model for SHO of pendulum for small angles. Include plot of θ versus t and ω versus t. Include phase-space plot (ω versus θ) on a separate plot. • Equations modeled: ∆θ = ω∆t ∆ω = − gθ l ∆t • Strategy: – use Euler or Euler-Cromer algorithm for ω(t) and θ(t) – use Vpython 2D graphics windows for graphing OR matplotlib. – Use natural units for plots t in units of T , ω in units of 2π/T . • Pseudocode outline: 1. Import visual and visual.graph 2. Set values for givens: – value for g, gravitational constant. – length of pendulum l – initial angle θ0 – initial angular speed ω0 – time step of simulation ∆t – time of simulation tmax 3. Calculate theoretical information for checking numerical accuracy. – predicted period T – a theoretical value for energy – amplitude of oscillation? 4. Start 2D vpython graph and start the θ versus t curve. Repeat graph and curve commands for each curve. 3 5. Do actual Euler calculation (listed here as separate function eu- ler calculate) Stop after reasonable number of oscillation periods. Collect points in lists or add plot points during simulation. 6. Finish up 4 • Pseudo-code for euler function (Vpython) INPUT: initial value θ, initial value ω, g/l, tmax, ∆t, names of each curve. set t = 0 set θ = initial value set ω = initial value while (t

planet.dviPhysics 336/536 – Computational Physics – Fall 2022Code plan for planet orbiting fixed Sun in circular orbit:• Goal: Calculate circular orbit of planet around fixed Sun and create3D...

Get Answer To This Question

Related Questions & Answers

Submit New Assignment