Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Modeling the Universe
ASTR 596 - Fall 2025

Project 3: Monte Carlo Radiative Transfer in a Dusty Star Cluster

San Diego State University

ASTR 596: Modeling the Universe Instructor: Dr. Anna Rosen Due: Friday October 17, 2025 by 11:59 PM
Working Mode: Pair Programming

Pair Programming Assignments:

Core algorithms must be independently implemented. You may share debugging strategies, test cases, and plotting code.

Learning Objectives

Upon successful completion of this project, you will be able to:

Core Physics Understanding:

Computational Skills:

Validation & Analysis:

Scientific Interpretation:

Project Overview & Scientific Context

Mission Brief

You’re modeling a young star cluster embedded in a dusty molecular cloud. Your goal: build a Monte Carlo radiative transfer code to predict how dust extinction affects the observed starlight. This project connects microscopic physics (dust absorption) to macroscopic observables (colors and extinction).

Required Background from Module 4

This project builds directly on Module 4’s comprehensive radiative transfer framework. Please ensure you have thoroughly reviewed the following materials from the course website’s Statistical Thinking - Module 4: From Photons to Information:

Key Concepts You Should Understand

Before starting this project, you should be comfortable with the following concepts:

The Physical Picture: Monte Carlo Radiative Transfer

Monte Carlo Radiative Transfer (MCRT) solves the radiative transfer equation through a stochastic process, meaning it uses randomness and probability to represent physical reality. Each “photon packet” represents a bundle of photons carrying a fraction of the source luminosity. The stochastic nature manifests in two key ways: (1) the random sampling of interaction distances from an exponential distribution, and (2) the random emission directions. Through the law of large numbers, these random processes converge to the deterministic solution of the radiative transfer equation.

Discrete Absorption: In this project, you’ll implement discrete absorption where each packet either:

  1. Deposits ALL its energy at a single interaction point (when the sampled optical depth τ\tau is reached), or

  2. Escapes with ALL its energy (if it reaches the boundary first)

There is NO gradual energy loss - it’s all or nothing! As N→∞N \to \infty, this stochastic Monte Carlo approach converges to the exact solution of the radiative transfer equation.

The Physical System

Grid Setup

Star Cluster Configuration

5 ZAMS stars with the following properties (Z=0.02Z=0.02):

StarMass (M⊙M_\odot)TeffT_{\text{eff}} (K)LL (L⊙L_\odot)RR (R⊙R_\odot)Spectral Type
130~38,500~119,000~7.7O7V
220~34,000~43,000~6.0O9V
310~25,000~5,500~3.9B0V
45~17,000~530~2.6B5V
52~9,100~16~1.6A5V

Note: Use your ZAMS models from Project 1 to calculate exact values.

Dusty Medium

Density Calculations:

Dust Model: Draine (2003a,b) Milky Way dust with RV=5.5R_V = 5.5. All data here: https://www.astro.princeton.edu/~draine/dust/dustmix.html

Note on units: Use CGS units as your default throughout this project.

Reading the Draine opacity file:

Band-Averaged Opacities

You’ll calculate Planck mean opacities for each band by integrating the Draine data weighted by each star’s Planck function:

⟨κ⟩band,star=∫λ1λ2κ(λ)Bλ(Teff,star)dλ∫λ1λ2Bλ(Teff,star)dλ\langle\kappa\rangle_{\text{band,star}} = \frac{\int_{\lambda_1}^{\lambda_2} \kappa(\lambda) B_\lambda(T_{\text{eff,star}}) d\lambda}{\int_{\lambda_1}^{\lambda_2} B_\lambda(T_{\text{eff,star}}) d\lambda}
(1)

Why Planck mean opacity? Each star emits radiation following its Planck function Bν(Teff)B_\nu(T_\text{eff}). The Planck mean weights the dust opacity by the actual spectrum of light emitted by the star. This gives the effective opacity “seen” by photons from that specific star. Since hot stars emit more blue light and cool stars emit more red light, each star experiences a different effective dust opacity in each band.

Important: Calculate separate opacities for each star-band combination. A 38,500 K star will have different band-averaged opacities than a 9,100 K star, even using the same dust model, because their emission spectra weight the wavelength-dependent κλ\kappa_{\lambda} differently.

Tip: When creating each star, compute and store ⟨κ⟩band\langle\kappa\rangle_\text{band} for each band in a dictionary or array as a star attribute for quick lookup during packet emission.

Primary Bands:

BandWavelength RangeCentral λ\lambdaNotes
B390-500 nm445 nmBlue optical
V500-600 nm551 nmVisual
K1.95-2.40 μ\mum2.19 μ\mumNear-infrared

Important: Convert all wavelengths to cm for CGS consistency (1 nm = 10-7 cm, 1 μ\mum = 10-4 cm)

Conversion from Draine file: The file provides the wavelength-dependent κdust\kappa_{\text{dust}} (per unit dust mass), so τ=κdustρdust\tau = \kappa_\text{dust} \rho_\text{dust}.

Stellar Sources and Band Luminosities

Update your star.py class from Project 1 to calculate band-specific luminosities. Each star emits differently in each wavelength band based on its temperature:

⟨L⟩band=Lbol×∫λ1λ2Bλ(Teff)dλ∫0∞Bλ(Teff)dλ\langle L \rangle_{\text{band}} = L_{\text{bol}} \times \frac{\int_{\lambda_1}^{\lambda_2} B_\lambda(T_{\text{eff}}) d\lambda}{\int_{0}^{\infty} B_\lambda(T_{\text{eff}}) d\lambda}
(2)

where the denominator equals σTeff4/π\sigma T_{\text{eff}}^4 / \pi (from integrating the Planck function over all wavelengths).

The Planck function is:

Bλ(T)=2hc2λ51ehc/(λkT)−1B_\lambda(T) = \frac{2hc^2}{\lambda^5} \frac{1}{e^{hc/(\lambda kT)} - 1}
(3)

Important: Each star’s band luminosity depends on its temperature. Hot stars emit predominantly in B-band while cool stars emit more in K-band. Calculate ⟨L⟩band\langle L \rangle_{\text{band}} for each star individually.

Key Physics You’ll Implement

Optical Depth and Absorption

The heart of Monte Carlo radiative transfer is determining where photons interact with dust:

Lpacket=Lband,totalNpackets whereLband,total=∑starsLstar,bandL_{\text{packet}} = \frac{L_{\text{band,total}}}{N_{\text{packets}}} \text{ where} L_{\text{band,total}} = \sum_{\text{stars}} L_{\text{star,band}}
(5)

(sum over all stars for THIS band only). This ensures uniform Monte Carlo statistics within each band’s simulation - every packet in that band contributes equally to the error regardless of which star emits it.

Critical Point: There is no partial absorption! Each packet either deposits 100% of its energy or 0%. This binary outcome, averaged over millions of packets, reproduces the continuous radiation field.

Packet Emission Geometry

While stellar radii are tiny compared to the box size (Rstar/Lbox∼10−7R_{\text{star}}/L_{\text{box}} \sim 10^{-7}), packets must start from the stellar surface:

Initial position: Uniformly sampled on stellar sphere

Initial direction: Your independently sampled isotropic direction ϕ\phi and θ\theta will be the packet’s propagation direction (i.e., unit vector n^\hat{n}):

n^=(sin⁡θdircos⁡ϕdir,sin⁡θdirsin⁡ϕdir,cos⁡θdir)\hat{n} = (\sin\theta_{\text{dir}}\cos\phi_{\text{dir}}, \sin\theta_{\text{dir}}\sin\phi_{\text{dir}}, \cos\theta_{\text{dir}})
(6)

Important note: Store n^\hat{n} as a photon attribute!

Packet Propagation

Since we’re modeling pure absorption (no scattering), photon packets travel in straight lines:

r⃗new=r⃗old+Δs⋅n^\vec{r}_\text{new} = \vec{r}_\text{old} + \Delta s \cdot \hat{n}
(7)

where Δs\Delta s is the distance to the next cell boundary or interaction point.

Tip: Handle the special case where a direction component is exactly zero (e.g., dx = 0) by setting the distance to that boundary to infinity rather than dividing by zero.

Expected Opacity Values

Band⟨κ⟩\langle\kappa\rangle (cm2^2/g)Why This Value?
B∼×104\sim \times 10^4Peak dust extinction near 220 nm
V∼7.3×103\sim 7.3 \times 10^3Still strong extinction
K∼1.5×103\sim 1.5 \times 10^3Wavelength ≫\gg grain size

These apply to dust mass, not gas mass!

Important Note: Your exact escape fraction values will differ from others due to random stellar positions within the cluster. The relative ordering (K>V>B)(K > V > B) must always hold, but absolute values depend on where stars are placed relative to the box boundaries.

Computational Strategy

Grid Resolution Recommendations

Why this strategy: Start with low resolution to debug your algorithm quickly. Once working, increase resolution to see physical trends emerge. For final science analysis, depending on performance, use high packet counts to minimize Monte Carlo noise and high resolution (643 or 1283) resolution balances memory usage with spatial resolution.

Performance Measurement

Before running 106 packets:

  1. Time your code with 1000 packets

  2. Calculate packets/second

  3. Estimate time for full runs

  4. If estimated time > 24 hours, profile your code to find bottlenecks (see Appendix A on using cProfile)

Implementation Phases

Phase 0: Core Infrastructure - Get Your Framework Working

Build your modular MCRT framework with proper project structure:

project3_mcrt/
├── src/
│   ├── __init__.py
│   ├── constants.py      # Copy from Project 1, ensure CGS units
│   ├── star.py           # Update from Project 1 with band luminosities
│   ├── zams.py           # Copy from Project 1 (Tout et al. formulae)
│   ├── utils.py          # Planck function, integration utilities
│   ├── dust.py           # Draine opacity processing
│   ├── grid.py           # 3D grid structure (CGS units)
│   ├── photon.py         # Packet propagation
│   ├── transport.py      # Main MCRT engine
│   ├── detectors.py      # Escape tracking, observables
│   └── mcrt_viz.py       # Visualization and plotting functions
├── data/
│   └── kext_albedo_WD_MW_5.5B_30
├── outputs/
│   ├── figures/
│   └── results/
├── docs/
│   ├── research_memo.md
│   └── growth_memo.md
├── tests/
│   └── test_validation.py
├── README.md
├── requirements.txt
└── project3_analysis.py  # Main analysis script

Important:

Create utils.py with:

Create mcrt_viz.py with:

Validation: Verify grid indexing and coordinate transformations are correct

Phase 1: Single Band Implementation

Start with one band (suggest V-band), then expand:

Part A - Single Star:

Part B - All Stars:

Phase 2: Multi-Band Analysis

Expand to all three bands and increase packet counts:

Phase 3: Science Analysis and Extensions

Complete required outputs and implement extensions (see Grading Rubric).

Quick Sanity Checks

Before diving into analysis, verify your code with these tests:

  1. Empty box test: Set ρdust=0\rho_\text{dust} = 0 everywhere →\to 100% of packets should escape

  2. Extreme opacity test: Set κ×ρdust\kappa \times \rho_\text{dust} very high →\to nearly 0% escape (Note: high opacity reduces runtime during debugging since packets will be absorbed quickly/close to sources)

  3. Energy conservation: ∣Lin−(Labs+Lesc)∣/Lin<0.001|L_{in} - (L_{abs} + L_{esc})| / L_{in} < 0.001 for all runs (note this may be higher for very low packet counts due to Monte Carlo noise and grid resolution).

  4. Statistical convergence: Error should scale as 1/Npackets1/\sqrt{N_{packets}}

  5. Single star centering: Star at origin in uniform medium should give fesc=e−τf_{esc} = e^{-\tau}

Common Debug Checks

If your results seem wrong:

  1. Verify τ=−ln⁡(ξ)\tau = -\ln(\xi) not +ln⁡(ξ)+\ln(\xi)

  2. Check you’re using ρdust\rho_\text{dust} not ρgas\rho_{gas}

  3. Confirm packets that escape are counted in LescapedL_\text{escaped}

  4. Test with one packet and print every step

  5. Verify your random number generator is working (should give uniform distribution)

Required Outputs & Analysis for Your Research Memo

  1. Opacity Validation Plot:

Validate your opacity calculations from the Draine data:

  1. Convergence Analysis:

Plot fescf_{\text{esc}} vs NpacketsN_{\text{packets}} for each band:

  1. Spectral Energy Distribution

Show how dust reddens the star cluster:

  1. Spatial Aborption Maps (2D Projections)

Create absorption maps for each band (B, V, K) by integrating through the full box:

$∑zLabsorbed, band(x,y,z)\sum_z L_\text{absorbed, band}(x,y,z)

Use log color scale to show dynamic range, choose your limits wisely.

  1. Escape Direction Map (optional)

2D histogram showing the angular distribution of escaping light:

6. Data Table

QuantityB-bandV-bandK-band
Band-averaged opacity κ\kappa (cm2^2/g)
Input luminosity (L⊙L_\odot)
Escaped luminosity (L⊙L_\odot)
Escape fraction
Mean optical depth

Analysis Discussion for Research Memo

Your memo should address both the physics and computational aspects. These topics are suggestions to help guide your analysis - pursue the aspects you find most interesting in your results.

Physics Results:

Numerical Analysis:

Interpretation:

Critical Implementation Notes

Key Points for Success

CGS Unit Conversions

Common Pitfalls to Avoid

  1. Using gas density instead of dust density (100% — error!)

  2. Sign error in tau sampling (packets never absorb)

  3. Lost packets at boundaries (energy not conserved)

  4. Variable packet luminosities (breaks statistics)

  5. Not integrating opacities properly (missing temperature weighting)

  6. Wrong units in wavelength conversion (nm/μ\mum to cm)

Validation Requirements

Before analyzing results, your code MUST pass:

Test 1: Uniform Sphere

Single star at center, uniform medium:

Test 2: Energy Conservation

∣Lin−(Labs+Lesc)Lin∣<0.001\left|\frac{L_{in} - (L_{abs} + L_{esc})}{L_{in}}\right| < 0.001
(9)

Test 3: Convergence Scaling

Standard error ∝N−1/2\propto N^{-1/2} for N∈[103,106]N \in [10^3, 10^6]

Grading Rubric

ComponentPointsFocus
Core Implementation40Correct MCRT algorithms, multi-band treatment, energy conservation
Code Design & Quality20Modular structure, documentation, efficiency, readability
Research Memo & Analysis25Physics interpretation, required visualizations, quality of discussion
Extension10Implementation quality and comparison with baseline
Validation Tests5Energy conservation, uniform sphere test, convergence verification

Extension Requirement (Choose 1)

Implement ONE extension and compare results with the baseline. The list below provides ideas to get you started, but you’re encouraged to propose your own extension that interests you. If you have an idea or want to discuss feasibility, come chat with me during office hours.

Physics Extensions:

Computational Extensions:

Observational Extensions:

Analysis Extensions:

Your own idea: Propose something that connects to your research interests or explores an aspect of radiative transfer you find intriguing.

Your extension must include quantitative comparison with baseline and discussion of physical implications.

Implementation Tips

Before you begin coding, these practical tips will save you hours of debugging:

Performance & Efficiency:

Debugging Strategy:

Numerical Robustness:

Getting Started

  1. Review Module 4 lecture notes

  2. Copy your constants.py, star.py, and zams.py from Project 1

  3. Start simple: Single packet, single star, uniform medium

  4. Build incrementally: Add complexity only after validation

  5. Test constantly: Every new feature needs a test

What Success Looks Like

Appendix A: Using Python’s cProfile for Performance Analysis

What is cProfile?

cProfile is Python’s built-in profiler that measures where your code spends time. It tracks every function call and measures execution time, helping identify bottlenecks.

Basic Usage

Method 1 - Command line:

python -m cProfile -s cumtime project3_analysis.py > profile_output.txt

The -s cumtime sorts by cumulative time spent in each function.

Method 2: In your script:

import cProfile
import pstats

# Profile a specific function
profiler = cProfile.Profile()
profiler.enable()

# Your MCRT code here
results = run_mcrt(n_packets=1000)

profiler.disable()

# Print statistics
stats = pstats.Stats(profiler)
stats.sort_stats('cumtime')
stats.print_stats(20)  # Show top 20 functions

Interpreting Output

bash - Unknown Directive
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1000    0.234    0.000    45.678   0.046   transport.py:45(propagate_packet)
128000  12.345   0.000    12.345   0.000   grid.py:23(get_cell_index)

What to Look For

  1. Functions with high cumtime - these are your bottlenecks

  2. Functions called millions of times with small percall time - consider vectorization

  3. Unexpected functions taking significant time - may indicate bugs

Example Optimization Workflow

# Before optimization: profile your code
# Identify that get_cell_index takes 30% of runtime
# Optimize that specific function
# Re-profile to verify improvement

Remember: Profile before optimizing! Don’t guess where the bottlenecks are.

Appendix B: Parallel Processing for Bands

Optional Speedup: Band Parallelization

Since B, V, and K bands are independent, you can run them simultaneously for ~3× speedup using Python’s multiprocessing module (click here for docs).

Key concept: Each band runs as a separate process with its own random seed and absorption map. After all bands complete, combine the results.

Considerations:

Starting point: Look into multiprocessing.Pool and its map() function to distribute bands across cores.

This optimization is entirely optional - focus on getting correct physics first.

Appendix C: Troubleshooting Checklist

If your results seem incorrect, check these common issues:

“Monte Carlo: Solving intractable integrals by rolling dice since 1946”

Project 2
Project 2: Science Background
Project 3
Project 3 MCRT - Starter Code Structure