For the problems below, you need to write programs using C language. Turn in your own source programs written independently (by email) together with reports.


It has often been customary to describe the initial mass function (IMF) of stars in a cluster using the Salpeter function of the form:

\begin{align} N(M) dM \propto M^{- \alpha} dM \quad \mathrm{with} \quad \alpha = 2.35 \end{align}

where $N (M )dM$ denotes the number of stars with mass between $M$ and $M + dM .$ Use $10^6$ random numbers to draw a Monte-Carlo sample from the Salpeter IMF with mass in the range of $1 - 100 M_\odot$. Plot the resulting probability distribution function as a function of $\log M$, and fit it to a power law, and compare the fitted power-index with $\alpha = 2.35.$


Purely Scattering Atmosphere: Suppose a photon incident on the bottom of a plane-parallel slab with optical depth $\tau_\mathrm{max}$. We assume that the slab is infinite in the $x$ and $y$ directions: $z$-axis is normal to the slab. The photon can be scattered (with no absorption) at any point within the slab. We begin with $z = 0$ (bottom of the slab) and follow the photon’s trajectories up to $z_\mathrm{max} = 1$ (top of the slab) by taking following steps:

  1. We inject a photon from below whose flux is isotropic in any direction. The probability for a certain injection angle at $z = 0$ (with respect to the z-axis), is given by $p(\mu) d \mu = 2 \mu d \mu$ with $0 \le \mu \le 1$, where $\mu \equiv \cos \theta$. Use a uniform deviate $\xi_1$ to sample $\mu$. Use another uniform deviate $\xi_2$ to obtain $\phi = 2 \pi \xi_2$. Calculate the initial direction of the photon $( \sin \theta \cos \phi , \sin \theta \sin \phi, \cos \theta).$ Note that the initial position of the photon is $(x, y, z) = (0, 0, 0).$
  2. When the optical depth of a photon is $\tau$, the probability $P$ that a photon interacts with the medium is given by $P = 1 − e^{−τ}.$ Noting that $P$ is random over $[0, 1)$, $\tau$ can be sampled from a uniform deviate $\xi_3$ as $\tau = - \ln (\xi_3)$. Pick up a value for $\tau$.
  3. The distance traveled by a photon along a ray is given by $L = \tau z_\mathrm{max} / \tau_\mathrm{max}.$ Update the photon’s new position as $x = x + L \sin \theta \cos \phi, y = y + L \sin \theta \sin \phi,$ and $z = z + L \cos \theta$.
  4. If $z < 0$, add one to the number of reflected photons. If $z > z_\mathrm{max}$ , add one to the number of transmitted photons. In either case skip to Step 6 below.
  5. Assume that the photons are scattered uniformly into $4 \pi$ steradians. Generate the new direction by sampling uniformly for $\phi$ in the range $0$ to $2 \pi$ and $\mu$ in the range $−1$ to $1$: $\phi = 2 \pi \xi_4$ and $\mu = 2 \xi_5 -1$, where $\xi_4$ and $\xi_5$ are two independent uniform deviates in $[0, 1)$.
  6. Repeat Steps 2–5 until the fate of the photon has been determined.
  7. Repeat Steps 1–6 with additional incident photons until sufficient data has been obtained.


Calculate the transmission and reflection probabilities for $\tau_\mathrm{max} = 0.01, 0.1, 1,$ and $10.$ Begin with $10^3$ incident photons and increase this number until satisfactory statistics are obtained. Give a qualitative explanation of your results.


Draw some typical paths of photons in the $x–z$ plane.


Consider a 2D Poisson equation in the $x \times y = [−L/2, L/2] \times [−L/2, L/2]$ plane

\begin{align} \left( {{\partial^2} \over {\partial x^2}} + {{\partial^2} \over {\partial y^2}} \right) u(x, y) = f(x, y), \end{align}

subject to the periodic boundary conditions: $u(x+L/2, y) = u(x−L/2, y)$ and $u(x, y + L/2) = u(x, y − L/2).$ Perhaps the most convenient way to solve this equation is to use Fourier transforms.


Use the indices $j$ and $k$ to indicate the grid points in the $x-$ and $y-$directions, respectively. (That is, $x_j = (−1/2 + j/N )L$ and $y_k = (−1/2 + k/N )L$ for $j, k = 0, \cdots , N .$) Write a finite difference representation of equation (2).



\begin{align} u_{j,k} = \sum_{m=0}^{N-1} \sum_{m=0}^{N-1} U_{m,n} e^{- 2 \pi i (m x_j + n y_k ) / L}, \end{align}
\begin{align} f_{j,k} = \sum_{m=0}^{N-1} \sum_{m=0}^{N-1} F_{m,n} e^{- 2 \pi i (m x_j + n y_k ) / L}, \end{align}

with integers $m$ and $n$. Note that i is the imaginary unit (i.e., $i^2 = −1$). Show that the result of Part (a) is reduced to

\begin{align} U_{m,n} = F_{m,n} / \lambda_{m, n}, \end{align}

where $\lambda_{m, n}$ is the 2D gravitational kernel. Find an expression for $\lambda_{m, n}$.


Calculate $u(x, y)$ for $f = e^{ −(x +y ) }$ on $x \times y = [−1, 1] \times [−1, 1].$ Make a contour plot of $u$ in the $x–y$ plane as well as a cut profile $u(x, 0)$ along the $x$-axis $(y = 0)$.