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.