Assignment set #2 - PHYD57 - Adv. Comp. Methods for Physics =========================================================== Due on 17 Feb 2022, 23:59pm Note: this is 2 days after normal due date, to leave 10 days for the solution. Problem 2.1 - Radial velocity curve _____________________________________ There is a planet, which we "see" via the reflex motion of its host star. Center of mass of the star-planet system is at rest with respect to us, the observers. Mass of the planet is suspected to be 1000 times smaller than the star, but that needs to be checked by Doppler spectroscopy of the stellar spectrum. In preparation, we need to be able to plot the theoretical radial velocity curve, Vr(t). The mass of the star is in practice equal to the system's mass, M = 0.5 solar masses, and from the periodicity of the radial motion and Kepler's laws astronomers have found that the semi-major axis of the elliptical relative orbit of the system equals a = 1.2 AU. Finally, e = 0.64. An ellipse in polar coordinate system (r,f), where r is the distance from the center and f the so-called true anomaly (angle between the vector drawn from the focus to the object travelling on ellipse, and the line of apsides, i.e. the line joining the periastron and apoastron (the closest and the furthest points of the orbit), has mathematical equation r(f) = a(1-e*e)/[1 + e*cos(f)], where e = nondimensional eccentricity parameter describing the shape of the ellipse, 0 <= e < 1. Both f(t) and r(t) change in time; true anomaly changes between 0 and 2 Pi radians. In elliptic motion, angle f depends on time in a complicated way, called Kepler equation. We introduce two angles called mean anomaly nt. The first is a re-scale time, the product of n = mean angular speed = const., and time variable t, counted from last periastron; mean anomaly is the phase of motion on an ellipse, expressed as angle in radians, from 0 to 2*pi. The second angle, eccentric anomaly E, is a non-uniformly increasing angle, obeying Kepler equation valid at any time t: nt = E - e sin(E). Over one hundred numerical methods have been tried since the times of J. Kepler to quickly solve this transcendental equation for E(t), or E(nt). When high accuracy is required and eccentricity is not very large, the equation can be solved by iteration of the kind E = nt + e*sin(E) where E in the argument of sine function is the previously obtained approximation, and the next estimate for E is the whole r.h.s., also found on the left-hand side. After many iterations, which are usually started by initially guessing E(t) = nt, the updates become very small, and eventually less than machine accuracy. You can stop the iterations when the increase (abs value) becomes less than 5e-16, or after 120 iterations, whichever comes first. Then you use the relationship of true anomaly f(t) on eccentric anomaly E(t): tan f/2 = tan(E/2) sqrt[(1+e)/(1-e)] to get f(t) for any given t. Every t requires separate iterative solution of Kepler equation. Once f(t) is found, you can plug it into equation of ellipse or do other interesting things with it. All this is by the way described in wiki pages on eccentric anomaly and Kepler equation. It is important to know that the full time derivative of f, df/dt, is the angular speed of the system, r df/dt is the transverse linear speed at right angle to radius, and r^2 df/dt is a constant angular momentum sometimes denoted by l or L, whose value is L = sqrt(GMa(1-e*e)) = const. So far we talked about velocity components of the *planet*. The motion of the star is a reflex motion and by assumption exactly 1000 times smaller, (both the distance and the velocity components) since the star is 1000 times closer to C.M. than the planet. For instance, if the planet moves at 30 km/s, the star moves at 30 m/s. Your task is to derive and plot the simulated radial velocity of the star as seen from Earth, when the line-of-sight is inclined by an angle w (it should be written as omega, here w resembles omega, the proper symbol, also called the longitude of periastron). The viewer may be at any angle with respect to the line of apsides. By projecting the velocity onto the line-of-sight, show that the line-of-sight component of velocity in elliptic motion, which a distant observer will call radial velocity, equals Vr = sqrt{GM/[a(1-e^2)]}*[ sin(w+f) (1+ e*cos(f)) - e sin(f) cos(w+f)] Using ** for raising to a power, like in Pytho or Fortran, Vr = 30 m/s [(M/M_sun)(AU/a)/(1-e^2)]**0.5 * *[sin(w+f) (1+ e*cos(f)) -e sin(f) cos(w+f)] Plot Vr(t) for M, a, e mentioned at the beginning, and 3 choices of angle w: w=0, w=1, w = 2.5 rad. You can overplot those 3 curves in one figure. Do computations in the range nt = 0...8 radians, use a few hundred points to display a little more than one full period of motion, since 8 > 2*Pi. Is the amplitude of the velocity curve (maximum minus minimum) independent of w? Include the plots. Please explain why the shapes are different and do not resemble a sinusoid. This problem can all be programmed in Python for easy plotting. Write it in a moderately modular way: have a self-standing procedure/subroutine or function, which solves the Kepler equation, and a separate procedure doing the plot. We'd also love to see how your little Kepler equation solver looks and works in Fortran or C(++) on some example value of n*t phase angle, to know that you are learning active use of one of these languages. Can you print during the iteration of order 20 lines stating iteration number, current best value of E, and the difference of that value with respect to the previous approximation. Doing graphics is not necessary. Work on your code bottom-up, first writing a tiny driver main program that states the values of and a procedure that computes E corresponding to nt = const = 1.0 rad. Verify that iteration works well during debug stage by printing values at every iteration. The needed input parameters (e, nt) must be transferred to or otherwise known to the subprogram; output parametars should be the number of iterations done, and the resulting E. Problem 2.2 - Satellite subject to constant radiation pressure ______________________________________________________________ This exercise can be done in any language (Python is OK). It involves solving two ODE equations of Newtonian dynamics using a leapfrog method. It should have been explained in a prerequisite course PSCB54. There's also wikipedia but it may give too many confusing variants. Just do like that: for every direction, first do the half-timestep update of positions (that's called "push"). Then evaluate the accelerations, and update the velocity components (do a full timestep "kick"). Finish the timestep with the remaining half-dt push. The two half-step pushes come one after another & so are sometimes coalesced into one full-step push to savearithmetic operations, except in the first and last step of the whole simulation. The name leapfrog refers to x,y,z leapfroging over vx,vy,vz: these quantities are never evaluated at the same time, they refer to time shifted by dt/2. In this exercise, don't do that. Write a function that evolves (updates) all variables over one full time interval dt. Return the updated variables to the calling program, unrestanding that the returned quantities are not in sync: positions are at end of timestep (t+dt), whereas the velocities at the midpoint t+dt/2. The simulated body is a solar sailboat capturing momentum from the sun's rays, and from the particles of the solar wind. It may look like a thin black sphere constructed from a completely absorbing, thinnest and lightest material possible, inflated by a small amount of air, to achieve maximum area-to-mass ratio of the spacecraft. The intention is that once placed in the starting, circular orbit 500 km above the surface of the Earth by a rocket, the solar sailboat will be picked up by solar radiation strong enough to carry it beyond the confines of Earth's gravity. Gravity of course has infinite range, but having captured enough momentum and energy from the wind, the spacecraft can escape, i.e. has a chance to never return to Earth. There are 2 levels of fidelity in the description of the motion of solar sailboat. You will explore only the simplest here: Level-1 description. We treat the Earth as an isolated body at rest, placed in a constant stream of solar radiation. In this approximation, sun's gravity and the motion of Earth around the sun are neglected. It's not a very bad assumption, in our case described below, but we won't know it until we consider a generalization that includes the third body (sun) in Level-2 description, which we defer to a future day during this course. Our spacecraft is subject to two forces and their accelerations: (i) Earth's gravitational acceleration g directed toward the origin of coordinate system at (0,0,0), with radial magnitude g = -GM/r^2, and correct direction. (Directions are expressed as versors, for any vector f, its direction is the unit vector a.k.a. versor: f/|f|), and (ii) acceleration of 'wind', directed to the right along x axis, with magnitude +f G M/(r0)^2, where G is gravitational constant, M mass of the Earth, r the distance from the center of the Earth, r0 is r(t=0), and f is a non-dimensional parameter of the wind force. The point of release of the spacecraft is y=-(R+h), R=6371 km (Earth's radius), h=500 km, and x=z=0 (find the other needed constants G and M online, or better yet, renormalize the equations of motion to the unit of initial orbital distance being 1, speed being 1, and gravity acceleration also 1; this way you'll also appreciate that the problem works the same way for all planets in the universe, not just Earth). The initial velocity of the object has zero y and z components, while v_x = +sqrt(GM/r0) is the Keplerian (circular) speed for an object subject to zero perturbing force. Given the vector sum of accelerations just described, formulate the differential equations in x and y coordinates (z coordinate will stay zero forever, since there are no forces in z direction, and initial velocity had zero z component). Choose timestep dt=const short enough that in the test case f=0 the position of the spacecraft after T=1.25 orbital periods is accurate to better than 1e-5 times the initial Earth-satellite distance R+h. Compute the theoretical position precisely. We use 1.25 or some other non-integer number of periods, because it may happen that a given numerical scheme produces much smaller positional error near the starting point than at subsequent points along the orbit. (Precise value of T is not important for error estimate, but obviously it must be the same in theoretical and numerical solution). Describe how you reached a good value of dt that is short enough but not unnecessarily short. Perform these numerical experiments: Try f=0.02 and observe what happens with the orbit. Monitor the current distance from the surface of the Earth and print out t=t_R when it becomes negative, but continue the calculation until the distance is less than 1/100 of the initial radius. At what time does that happen? Then increase the influence of solar radiation to 13.5 percent of Earth gravity force on the starting orbit (f=0.134). Plot the trajectory. You should observe a qualitative difference. What happens with the orbit? Finally overplot (or plot separately, but in Python overplotting should be very easy) the trajectory with f=0.2. The above calculations and plotting should be restricted to a square box with side equal to 6*r0 (coordinates in range +-3 r0). Finally, find as accurately as you can (however, consider more than 3 accurate significant digits an overkill) the critical value of f, above which the spacecraft escapes from the vicinity of the planet, given it's initial position. For this exploration, you may need to allow the computational box to be a few times wider, to allow for possible trajectories, which reach a distance larger than x,y = +-3*r0 and yet turn back and hit the surface of the Earth after that. Just to be safe in your conclusions! Include a thoroughly commented listing of the program, and the plot or precise hand-drawn sketch of 3 trajectories; provide answers to all questions. Problem 2.3 - Protons and positrons ___________________________________ This is how physicists managed to arrange at time t=0 two positrons (x) and two protons (o), at rest in a perfect square of side length d = 100 fm (femtometers). They release them at time t=0. o x x o Particles have electric charge opposite to electron charge. Gravity and other forces can be neglected. Show that because of symmetry, particles trace semi-infinite straight lines coinciding with diagonal directions of a square. Find particle trajectories numerically as distance from square's center vs. time at t > 0. Use symmetry to reduce the amount of work. At t--> oo (oo is infinity symbol in Python), find the ratio of positron to proton speeds, and the ratio of their kinetic energies. Use Fortran or C. Comment the codes thoroughly. Monitor the accuracy by tracking the total energy E (kinetic and potential energy of the system) and printing from time to time during the calculation the relative error in energy dE/E0 = (E-E0)/E0, where E0 = E(t=0). Either vary timestep of your integration scheme as a function of distance between particles, or keep it constant for simplicity, but assure (and describe) a reasonably good accuracy of total energy conservation. Integration scheme is up to you; motivate your choice. Compare your answer with the simplified analytical solution done with paper and pencil as follows. Assume that protons, being much more massive (find the physical laws and constants on wikipedia), move so slowly compared with the lighter positrons, that the latter actually 'think' that protons are always at rest, while protons 'think' that electrons have escaped instantaneously, before they had a chance to move much from rest. Assume energy conservation law for e+ and p+ particles separately in the analytical calculation. Hint: The given value of d is unimportant for the shape of trajectories and the asymptotic ratio of v_e+ / v_p+. Can you convince yourself of that? How can you simplify the equations to reflect this: do you really need all the physical constants in correct SI units in the equations? Do you need x and y coordinates or just one along the trajectory? Maybe it is enough to solve for the motion of two types of masses with the correct mass ratio, starting at (square side) distance = 1, under force law F = +1/r^2 between the particles, and of course using the diagonal symmetry. If you're not convinced, solve the problem the hard way first, with dimensional constants, precise laws of electrostatics etc. NOTE ON PLOTTING In problem 2.1 plots are required but can be done in Python. In problem 2.3 plotting of the trajectory is not required, but some plots may be helpful. Future assignments will ask for plots to accompany calculations done in C/F. 1. Start learning the usage of DISLIN graphics package for Fortran, since it is installed and working on the art-2 server (example will be provided in T6), or some other graphics (instructions how to install DISLIN on your computer and manuals are available online). 2. Alternatively, and easy for you, since you are familiar with PyPlot, write results to a file, which you'll read with a Python script and display with matplotlib.pyplot. NOTE ON COMPUTATIONS This set requires some limited programming in C or Fortran, in addition to Python. Install the Windows compilers for one or more these languages, or use Linux on your own machine (familiarity with Linux is a required outcome in this course), or the art-2 server, account phyd57@art-2.utsc.utoronto.ca with the password communicated in lectures and tutorials. Work in your own subdirectory. Compilers: gcc, g++, gfortran, ifort, icc, icpc, pgf95, pgcc, pgc++, are available. Remember to use sftp (safe file transfer protocol) or other secure ftp program from your computer, in order to transfer files both ways. After you have finished working on your file on art-2, save it on your local machine & log out (CTRL-D or 'exit' are normal ways), after delete assignment files from subaccount. First, we do not provide backups of your work, and secondly your work would remain visible to other students. We do pay attention to possible plagiarism in submissions but leaving source files visible is unnecessarily enticing and you don't want to get into a dispute on who created the implementation first etc. Finally, the basic algorithmic ideas and workflow of all your codes must be clearly explained in your own words through extensive blocks of comments inside. C/F aourse files must compile on art-2. Python must run on your machine, you may submit to Qurcus a jupyter notebook with results if you wish, or a .py file plus description.