MS2 UE4: Astrophysics with the computer
Hubert Baty, Joachim Köppen, Strasbourg, 2012/13
Please note
 choose your subject before the session on 9 Nov. 2012, so that you
can start working on it right away ... time runs quickly ... and the
number of afternoons is limited ...
 do not hesitate to contact Hubert Baty or me (joachim.koppen at astro....)
for more information about the subjects, and if you have any questions!
Ne hesitez pas de nous contacter!
 in this TP no question is stupid
 on parle le français ou l'anglais, comme vous le preferez!
 the programming language is your choice: you can use FORTRAN,
C, or JAVA
 for those who like to use JAVA, I shall give you two graphics libraries
which permit you to create xy plots in the usual scientific form.
One library, MorganPlot, is very good for looking at results after they
have been computed; the other, JKZoomPlot, is better for displaying
data during the simulation (such as nbody simulations).
Contact me for any questions.
All projects are of similar difficulty and complexity, but they may
differ in numerical effort (I give a rough indication by * to ***** from
'straightforward' to 'tough', but all have their numerical delicacies!).
All leave space for exploration of the physics or the numerics, and all
are suitable for modifications and additions, depending on your personal
interest. And all can be done within the time available.
Subjects supervised by Hubert Baty:
 Exploration numérique du problème à 3 corps restreint en interaction
gravitationnelle (**)
 Développement d'un code hydrodynamique 1D: traitement des chocs (**)
 Etude numérique de l'équation de GinzburgLandau: instabilité convective
et absolue (**)
 Prise en main et tests d'un code 1D/2D MHD (***)
a basic applet for the GinzburgLandau equation
Subjects supervised by Joachim Köppen:
 Onedimensional Hydrodynamics: compute the
time evolution of a fluid in a tube, starting with given initial
profiles of density and velocity. Showing the formation of shocks
and rarefraction waves. With the addition of the energy equation,
the temperature profile can be followed, giving rise to the contact
discontinuity. The program can be modified for spherical geometry
and with the additional force of selfgravity, the collapse of a
protostar or galaxy can be calculated.
(***  Explicit solution of a system of partial differential equations,
artifical viscosity, conservation of mass and momentum). The basic
project can be directed to several different applications and
extensions, as described in the dossier.
the applet
 Smoothed Particle Hydrodynamics: compute the
time evolution of a fluid in a tube, starting with given initial
profiles of density and velocity. This will show the formation of shocks
and rarefraction waves. With the addition of the energy equation,
the temperature profile can be followed, giving rise to the contact
discontinuity. This problem is treated by a
particlebased method: A large number of particles moves in space following
certain acceleration laws in such a way that their smoothed density
behaves as closely as possible the behaviour of the density in a real
fluid. Because there are no fixed spatial grids on which the functions
are evaluated, the program can easily be modified for two and three
space dimensions, without any constraint on the geometrical symmetry.
(***  solution of a system of ordinary differential equations).
The basic project can be directed to several different applications
and extensions, depending on the time available.
the applet
 NBody Simulations: compute the evolution
of the positions of one or several planets or stars due to their
gravitational interaction with the central sun or the other stars.
(*  Solution of a system of ordinary differential equations with constant
and adaptive time steps, treating the singularity at zero separation).
The basic project can be directed towards several different applications,
such as mentioned in the dossier. For example, one can verify the Jeans
criterion for the gravitational collapse of a cloud of stars or galaxies.
Orbits in JavaScript
Two Body applet
Three Body applet
Special cases: Three Stars of Equal Mass
 Saturn's Ring: write a code to perform full
nbody simulations which compute the evolution of the satellites in orbit
around a planet due to their gravitational interaction with the planet and
among each other. This will check a criterion for stability of such a ring
of moons, found analytically by J.C.Maxwell in 1856.
(**  Solution of a system of ordinary differential equations with constant
and adaptive time steps).
The basic project can be extended towards several different more general
cases, such as mentioned in the dossier.
links to a JAVA applet also doing the full
nbody problem
 Orbits under MOND:
The postulate of mysterious dark matter  which nobody has identified yet 
to explain the observed rotation curves of galaxies can be avoided if we
accept that Newtonian dynamics is modified under the large dimensions
in the realm of galaxies. Let us write a code that explores what type of
orbits are possible if we suppose the MOND formulation, and compare them
with the usual Kepler orbits of standard Newtonian gravity. We can treat
the onebody problem, where small test particle moves in a fixed potential,
and the twobody problem of two masses orbiting each other. More complicated
cases are unfortunately beyond the simple methods we can employ!
(**  Solution of a system of ordinary differential equations with constant
and adaptive time steps).
links to JAVA applets which do the one and twobody
problem in Newtonian gravity and MOND
 Circular and Escape Velocities from a
Galactic Disk: derive these velocities from a given mass distribution
of stars in a galaxy. (***  Integration of the contributions to the
force acting on a test star by all other mass elements. Treating the
singularity of Newton's gravity law at zero distance. Integration
of the force out to infinity to get the potential)
 Finding the Gravitational Potential
of a Galaxy: compute the gravitational potential in and around
a given mass distribution. (****  Solution of Poisson's equation by
relaxation methods (Jacobi, GaussSeidel, Simultaneous OverRelaxation);
also possible: by solution of system of linear equations by Gaussian
elimination). This project can be split up in several projects with
different methods and/or models.
 Ram Pressure Stripping of Disk Galaxies:
in a simplified model we consider the fate of a point mass held in a
potential which is subjected to an impulsive external force, to
determine under which conditions it escapes from the potential.
(*  Solution of ordinary differential equations, systematic exploration
of the parameter space)
the applet
 Monte Carlo Radiative Transfer:
We follow the flight of a photon through a medium, starting at
its emission at the source, going through all absorptions and
reemissions, until it escapes and reaches the observer. If we do
this with a sufficiently high number of photons, we can compute
the observable intensity for any arbitrary geometry of the source
and the interacting medium and for any amount of scattering 
and any law of scattering. Let us try out this powerful method to
solve radiative transfer problems!
(**  Monte Carlo simulation, transformation of uniform random
numbers into those distributed with given functions)
the applet
 Temperature Fluctuations in Ionized
Nebulae: compute the temperaturesensitive intensity ratio of
the forbidden lines [O III] 5007/4363 Angstrom from a nebula in which
the temperature fluctuates with a certain amplitude about a central
value, and derive the apparent "mean temperature". Comparison with the
true mean value illustrates that in the presence of such fluctuations
the interpretation of nebular spectra is systematically affected.
(**  Solution of the ionic rate equations by Gaussian elimination, Monte
Carlo simulation of the temperature, transformation of uniform random
numbers into those distributed with a given function, fitting a formula
to numerical results)
the applet
 Timedependent Ionization of Gaseous
Nebula: HIIregions are the ionized interstellar gas around hot
main sequence stars which live only a few million years. Planetary
nebulae are the expanding ejecta from a star which becomes very
hot before it becomes a white dwarf. The hot stars ionize the gas,
and normally we assume ionization equilibrium to interpret the
emission line spectrum. But what happens when such a star switches
on (or off) very rapidly? Under which conditions we should expect
deviations from equlibrium? Let us follow the evolution of the ionization
for any given timehistory of the ionizing stellar flux, including those
from published stellar evolutionary tracks.
(*  Solution of ordinary differential equations, adaptive time step,
interpolation, systematic exploration of parameter space)
The basic project can be extended towards other
cases, such as mentioned in the dossier.
the applet
 Deflection of Cosmic Rays by Earth magnetic field:
compute the trajectories of cosmic ray particles in the magnetic field of
the Earth, starting from various initial positions  far from the Earth 
and different initial velocities, in order to find out under which conditions
the particle can reach the Earth surface. With these simulations one can
show which parts of the Earth are more hazardous to these rays, and one can
make a world map map of the probability, i.e. the cosmic ray flux as a function
of energy. These simulations can be compared to the results of the pioneering
study of this subject, done in the 1930s, long before computers
(***  Solution of ordinary differential equations, automatic time stepping,
systematic exploration or Monte Carlo simulation for the initial positions
and velocities, transformation of uniform random numbers into those distributed
with a given function)
a basic applet
 Excitation in and Emission from a Molecular Cloud:
The CO molecules in an interstellar cloud are excited by collisions
with hydrogen molecules into excited states, from which they relax
by producing line emission in the microwave region. The emission can
be absorbed by other CO molecules, leading to their excitation.
The emission from the entire cloud is the result of the balance between
excitations and deexcitations of the 'average' molecule. Under
certains conditions of gas density and size of the cloud the upper
levels of the CO molecule can be overpopulated, and the cloud will
become a MASER, which is able to amplify any radiation that comes
in from the exterior, e.g. the cosmic microwave background.
(**  Solution of system of linear equations; iteration between
radiation field and molecular states)
the applet
 Perturbations of planetary orbits: Newton's law of gravity
predicts that the orbits of two bodies moving around each other, such as a planet
around the Sun, are perfect conic sections. However, the presence of another planet
will cause the first planet to deviate from the Keplerian orbit. In our Solar System,
the massive planets Jupiter and Saturn have such an effect on all other planets, and
also on any spaceprobe that we send away. The presence of Neptune was predicted by
LeVerrier from its perturbing effects on Uranus' orbit. With our present computers it
is easy to compute numerically any planetary orbits, and also include the forces from
the other planets. Thus we can compute these rather small perturbations by comparing
the actual trajectory with the hypothetical one, where we had neglected these additional
forces. We shall do this for a simplified system of two planets moving around the
sun, in order to focus on the essential effects.
(**  Solution of ordinary differential equations, automatic time stepping,
assuring low numerical errors)
the applet
 Precession of orbits: The extraordinary property of
Newton's law of gravity is that all orbits (in the twobody problem) are closed
orbits, perfect conic sections. However, any deviation from this perfect case
causes elliptical orbits to precess, that is their major axis performs a slow
revolution about the central body. Newton himself had worried about deviations
from the inversesquare law of the force, and shown how the precession rate is
linked to such a deviation. We know that precession can occur due to perturbations
by other planets, deviations of the central body's shape from a sphere, any
external tidal forces and due to relativistic effects ... In more recent times,
observations of galactic rotation curves suggest either that we have to accept the
presence of enigmatic dark matter if we want to stick to Newton's gravity law or
we have to modify Newtonian dynamics. The latter approach (MOND) implies that orbits
on galactric scales would deviate from Kepler orbits and show a precession. In our
exercise, we want to explore under which situations this precession would make itself
detectable. Could we measure it by observing certain binary stars with high precision?
Which kind of binaries and which precision would be needed?
(**  Solution of ordinary differential equations, automatic time stepping,
assuring low numerical errors)
a basic applet
 Bayesian Estimation of Parameters:
Measuring the Radio Flux from the Moon:
The temperature on the surface of the moon can be measured by observations
at radio waves: The flux is proportional to the temperature. This can be
done by executing a drift scan, where the moon passes through the main
lobe of the radio telescope's antenna, due to the rotation of the earth.
The signal goes up, then down again, in a bellshaped curve. Its height
would give the flux ... however, the data are noisy, as observational data
often is! A good and efficient method to extract as much information as
possible from a data set is to compare it with a model for the physics,
the observational process, and the noise. In a systematic way, we scan
through all possible combinations of relevant parameters. By assuming
a certain type of noise (Gaussian) we can compute the likelihood that
the observed data is the result adding noise to the model with a
given set of parameters, and hence we can construct maps of the
probability for the model parameters. In this way, we can find the
most probable height of the curve  which will give the radio flux 
as well as its error bar, i.e. the range of its value for a given
percentage.
(*  generation of random numbers, efficient scanning of multidimensional
parameter space, interpolation, numerical integration)
an applet (fit Gaussian to noisy data)
 Ionization structure of an HII region:
Stars hotter than about 30000 K emit plenty of photons with energies
above 13.6 eV, the ionization energy of hydrogen. As a consequence,
any gas around the star will be ionized up to some distance which is
determined by the actual number of hard photons emitted per second.
We place such a star in a large gas cloud of uniform density, and
compute the state of ionization of hydrogen as a function of
distance from the star. As the distance increases the photon flux
decreases, and so does the ionization, until the gas becomes neutral.
It is possible to compute also the ionization of helium or other heavy
elements, and to compute the strength of emissions lines.
(**  Solution of ordinary differential equations, automatic radial
stepping, integration over stellar spectrum, fix point iteration)
 The quantum states of the Hydrogen atom:
Solving Schrödinger's equation for an electron in the electrical
field of a proton yields the wave function of the electron, and gives
the bound states of the electron. With the wave functions one can
compute the transition probabilities for the spectral lines, such
as the Balmer lines.
(*  Solution of ordinary differential equations, automatic radial
stepping, scanning through the energies and iteration to find the
eigenstates)
a basic JAVA applet
 Internal structure of a star:
stars are spherical gas clouds in hydrostatic equilibrium,
held together and compressed by their own gravity, and with
densities and temperatures high enough for thermonuclear reactions.
With these known physical processes and using the known properties
of the matter (equation of state, opacity, energy production rate)
it is possible to compute the structure of stars of given mass.
With guess values for the density and temperature at the centre
one integrates the equations radially outwards, and changes the
guess values until the conditions of low density and temperature
at the surface are met.
(**  Solution of ordinary differential equations, automatic radial
stepping, manual or automatic adjustment of initial conditions in
order to meet surface conditions)
a simple JAVA applet
 Evolution of viscous disks:
When matter falls into the gravitational well of a compact object,
like a black hole, a star, or the central region of a galaxy, it
forms a disk, due to its initial angular momentum. In this disk
of dense gas, hydrodynamic interactions between the different
layers apparently give rise to viscosity, and conversion of kinetic
energy into thermal energy. In this way the gas can lose the angular
momentum by emission of radiation. We shall consider the thin disk
model as it is used to interpret the behaviour of Xray binaries.
It is possible to add the recipes that explain the periodic
outbursts of cataclysmic variables.
(***  Solution of partial differential equations)
 Interpretation of emission lines:
radial velocity and internal kinematics of a planetary nebula.
Emission lines from ionized nebulae are usually optically thin,
therefore their shape tells about the distribution of radial
velocities on the line of sight. The wavelength of the line centre
give the radial velocity of the nebula with respect to us, which is
useful information of its motion within the galactic context. The
line width gives information about the internal motions of the gas
in the nebula, i.e. thermal and turbulent motions as well as any
systematic motion of the gas with respect to the centre of the nebula,
such as expansion of the nebular shell. We have some real observational
data, consisting of nine spectra taken at different positions over
the face of a planetary nebula. To obtain the parameters of an
observed emission line, we compare it with gaussian profiles
of systematically changed parameters (width, height,
central wavelength etc) to compute the likelihood of this fit.
This Bayesian techniques allows to determine the probability for
e.g. the width of the observed line, hence its most probable
value and the confidence interval.
(*  Generation of random numbers, Generation of simulated
noisy observational data, Scanning through a large parameter space,
Integration of marginal probability distributions)
Observational data: part 1, part 2
 Radio continuum emission from the Sun:
Given a simple model for the solar atmosphere  i.e. the density and
temperature as a function of height above the photosphere  one can
compute the emission from the sun at various radio wavelengths. Since
the atmosphere is a plasma, the frequencydependent opacity due to the
electrons need to be taken into account. The emission itself is due to
freefree emission by the electrons (or "Bremsstrahlung"). The emitted
intensity is obtained by integrating over various lines of sight, thus
allowing to construct a radioimage of the sun. This project could be
developed further by taking into account the refraction of the waves in the solar
atmosphere ...
(*  Computation of opacities and emissivities, integration along
lines of sight through a spherical atmosphere)
a simple JAVA applet
 Evolution of a spherical gas cloud embedded in hot gas:
What happens to a gas cloud in a thin hot medium? Under certain
conditions the gas will get heated and the cloud will evaporate.
Under other conditions the hot gas at the surface could cool because
of thermal conduction, and condense onto the cloud. There are a number
of situations in the interstellar medium and in winds from stars where
the fate of such a blob of gas is of interest. Here we want to solve
the fluid equations for a onedimensional system in order to see
whether we can characterize the physical conditions for condensation
and evaporation, and compare our findings with criteria by other
authors.
(***  Solution of partial differential equations, automatic time
stepping)
 Formation of absorption lines in a stellar atmosphere:
The decrease of temperature with height in a stellar atmosphere gives
rise to spectral lines to appear in "absorption" in the spectrum of
the star. For a simple formulation of the stellar atmosphere  the
heightdependency of the temperature  and a simple representation
of the optical properties of the gas  continuum and line absorption 
one can compute the emergent spectrum. This allows to study how lines
are formed and how their strength depends on the parameters of the
atmosphere and the atom.
(**  Computation of absorption coefficients, optical depths,
Integration along the line of sight through the atmosphere,
efficient scanning of the line parameter)
simple JAVA applet
 Dynamical evolution of an HII region:
When a massive star is born and lights up in the interstellar gas
that surrounds it, the emergent photons will start ionizing this
gas. The gas becomes transparent and therefore the photons can
ionize gas further away. This progressive ionization of its
neighbourhood continues until an equilibrium is reached between
the star's photon output and the recombinations taking place in
the ionized region, called a Strömgren sphere.
(****  Solution of ordinary differential equations, numerical
integration over stellar flux, iteration for ionization,
automatic predictor/corrector method for radial integration,
automatic time stepping)
 Chemical evolution of the solar neighbourhood:
By thermonuclear fusion in their central parts, stars are not only
generating energy that makes them shine, but they also produce heavy
elements from hydrogen and helium. While massive stars produce most
elements like oxygen and iron, intermediate mass stars also contribute
significantly to helium, carbon, and nitrogen, which are ejected into
the interstellar medium with an appreciable time delay because of the
longer life times of these stars. We want to construct a code that
computes the evolution of the chemical composition of the interstellar
gas dues to the birth and death of stars, taking into account the details
of the production of elements in stars of various mass ranges.
Other ingredients, such as the inflow of gas of arbitray composition
or the outflow of interstellar gas due to supernova explosions etc
can be added.
(***  Solution of ordinary differential equations, numerical integration
of the contributions to the elements, interpolation in tables)
Applet: One Zone Model
Simple Infall Model
One Zone Model with Gas Flows
Applet with Exponential Infall
One Zone Model: Oxygen and Nitrogen
Applet with fully detailed treatment
(a bit complex)
 Big Bang Nucleosynthesis (underway)
 to my HomePage

last update: 12 March 2013 J.Köppen