Project 3: Starter Code Structure

COMP 536 | Short Projects

Author

Dr. Anna Rosen

Published

April 22, 2026

This page provides the starter scaffold for Project 3. Copy these files into your repo and fill in the implementations.

Quick Start

From your repo root:

python run.py test          # Run pytest
python run.py validate      # Run numerical checks
python run.py make-figures  # Generate plots

If those run cleanly, you’re most of the way to a good grade.

Repo Structure

Your submission must follow this structure:

project3_mcrt/
├── run.py                  # CLI entry point (required)
├── src/
│   ├── __init__.py
│   ├── constants.py        # Physical constants (CGS)
│   ├── star.py             # ZAMS star with band luminosities/opacities
│   ├── zams.py             # From Project 1
│   ├── utils.py            # Planck function, integration helpers
│   ├── dust.py             # Draine opacity file reader + Planck mean
│   ├── grid.py             # 3D Cartesian grid
│   ├── photon.py           # Photon packet + ray marching
│   ├── transport.py        # Main MCRT loop
│   ├── detectors.py        # Escape tracking, analysis
│   └── mcrt_viz.py         # Visualization functions
├── data/
│   └── kext_albedo_WD_MW_5.5A_30_D03_all.txt
├── outputs/
│   ├── figures/            # Generated plots (not committed)
│   └── results/            # CSV results (not committed)
├── tests/
│   └── test_validation.py
├── research_memo.pdf
├── growth_memo.pdf
├── README.md
└── requirements.txt

Docstring Convention

Use numpy-style docstrings for all functions and classes. See examples in the Python Fundamentals module.

run.py

"""COMP 536 Project 3 — Command-line interface.

This file dispatches to your implementation modules.
Keep it simple: real logic belongs in your modules, not here.
"""
import argparse
import subprocess
import sys
from pathlib import Path


def main():
    # ------------------------------------
    # 1. Set up argument parsing
    # ------------------------------------
    parser = argparse.ArgumentParser(
        description="COMP 536 Project 3: Monte Carlo Radiative Transfer"
    )
    parser.add_argument(
        "command",
        choices=["test", "validate", "make-figures"],
        help="test | validate | make-figures"
    )
    args = parser.parse_args()

    # ------------------------------------
    # 2. Ensure output directories exist
    # ------------------------------------
    Path("outputs/figures").mkdir(parents=True, exist_ok=True)
    Path("outputs/results").mkdir(parents=True, exist_ok=True)

    # ------------------------------------
    # 3. Dispatch to the appropriate action
    # ------------------------------------
    if args.command == "test":
        # Run pytest and exit with its return code
        sys.exit(subprocess.call(["pytest", "-v"]))

    elif args.command == "validate":
        # TODO: Import your modules and run validation checks
        # Example:
        #   from src.transport import run_mcrt
        #   from src.grid import Grid
        #   ... run single-star V-band, check f_esc against e^{-tau} ...
        #   ... check energy conservation ...
        raise NotImplementedError("Implement your validation checks")

    elif args.command == "make-figures":
        # TODO: Import your plotting module and generate figures
        # Example:
        #   from src.mcrt_viz import plot_convergence_analysis, plot_sed
        #   plot_convergence_analysis("outputs/figures/convergence_analysis.png")
        raise NotImplementedError("Implement your figure generation")


if __name__ == "__main__":
    main()

src/constants.py

"""Physical constants in CGS units.

Copy from Project 1. Add any new constants you need for Project 3.
Document each constant with units and source in your README.
"""

# Solar reference values
MSUN = None  # Solar mass [g]
RSUN = None  # Solar radius [cm]
LSUN = None  # Solar luminosity [erg/s]
TSUN = None  # Solar effective temperature [K]

# Fundamental constants
G = None        # Gravitational constant [cm^3 g^-1 s^-2]
SIGMA_SB = None # Stefan-Boltzmann constant [erg cm^-2 s^-1 K^-4]
C = None        # Speed of light [cm/s]
H = None        # Planck's constant [erg s]
K_B = None      # Boltzmann constant [erg/K]

# Unit conversions
PC = None  # Parsec [cm]

Module Responsibilities

You design the classes, functions, and interfaces. The Planning Worksheet will help you think through the contracts before you code. Below is what each module should handle.

star.py

Enhanced from Project 1. Each star needs:

  • ZAMS properties from zams.py (mass, \(L_\mathrm{bol}\), \(R\), \(T_\mathrm{eff}\))
  • Position \((x, y, z)\) in cm
  • Band luminosities \(L_\mathrm{band}\) for B, V, K (fraction of \(L_\mathrm{bol}\) in each band)
  • Planck mean opacities \(\langle\kappa\rangle_\mathrm{band}\) for each band (depends on \(T_\mathrm{eff}\))

utils.py

Physics helper functions:

  • Planck function \(B_\lambda(T)\) in CGS
  • Band integration (use scipy.integrate.simpson or similar)

dust.py

Read and use the Draine opacity data:

  • Parse the data file (watch the header format)
  • Compute Planck-weighted mean opacity for a given stellar temperature and wavelength band
  • Define band wavelength limits (B, V, K)

grid.py

3D Cartesian grid:

  • Uniform cells spanning \([-L/2, L/2]\) in each dimension
  • Cell indexing: position \(\to\) \((i_x, i_y, i_z)\)
  • Boundary checking
  • Gas/dust density storage
  • Absorption tracking arrays

photon.py

Photon packet creation and propagation:

  • Packet attributes: position, direction, luminosity, opacity, band
  • Isotropic direction sampling
  • Distance to next cell boundary (the hardest function — plan it carefully)
  • Ray marching through the grid

transport.py

Main MCRT simulation loop:

  • Luminosity-weighted star selection for packet emission
  • Run simulation separately for each band (B, V, K)
  • Energy conservation checking
  • \(L_\mathrm{packet} = L_\mathrm{band,total} / N_\mathrm{packets}\)

detectors.py

Track escaping packets and compute results:

  • Escape fractions per band
  • Escaped/absorbed luminosity totals

mcrt_viz.py

All required plots:

  • Opacity validation (Draine curve with band-averaged values)
  • Convergence analysis (\(f_\mathrm{esc}\) vs. \(N_\mathrm{packets}\))
  • Absorption maps (2D projections for each band)
  • SED comparison (intrinsic vs. observed)

tests/test_validation.py

Validation tests (run via python run.py test):

  • Empty box: all packets escape (\(f_\mathrm{esc} = 1\))
  • Uniform medium, single star: \(f_\mathrm{esc} \approx e^{-\tau}\)
  • Energy conservation: \(|L_\mathrm{in} - (L_\mathrm{abs} + L_\mathrm{esc})|/L_\mathrm{in} < 0.001\)
  • Convergence scaling: error \(\propto 1/\sqrt{N}\)

Conceptual Approach: Band-Separated MCRT

The Key Insight

Run the Monte Carlo simulation independently for each wavelength band (B, V, K) rather than mixing all wavelengths in one complex simulation. Each band gets its own complete MCRT run with its own packets.

Why This Approach?

  1. Simpler Implementation: No complex multi-dimensional weighting needed
  2. Easier Debugging: Get one band working perfectly before adding others
  3. Clearer Physics: Each packet represents “monochromatic” radiation
  4. Natural Development: Implement V-band in Phase 1, then extend to all bands in Phase 2

Important Considerations

  • Each star emits different amounts of energy in each band (based on its temperature)
  • Each star has different opacity in each band (Planck mean weighted by stellar spectrum)
  • Within a single band’s simulation, all packets carry equal energy
  • Star selection is weighted by that star’s luminosity in the current band only
  • The three band simulations are completely independent — they can even run in parallel

Development Strategy

Phase 1: Implement complete MCRT for V-band only Phase 2: Loop over all three bands, storing results separately Analysis: Compare escape fractions across bands to see differential extinction

Required Output Format

Data Table (save as outputs/results/escape_fractions.csv)

# band, opacity_cm2_g, L_input_Lsun, L_escaped_Lsun, f_escape, mean_tau
B,VALUE,VALUE,VALUE,VALUE,VALUE
V,VALUE,VALUE,VALUE,VALUE,VALUE
K,VALUE,VALUE,VALUE,VALUE,VALUE

Figure Naming Convention

All figures must be saved in outputs/figures/ with these names:

  • opacity_validation.png
  • convergence_analysis.png
  • absorption_map_B.png
  • absorption_map_V.png
  • absorption_map_K.png
  • sed_comparison.png
  • escape_directions.png (optional)

Important Notes

Structure is suggested, not rigid

This document provides a suggested code structure to help you organize your project. You may deviate if you have good reasons, but document significant changes in your README.

Implementation reminders

  • Run MCRT separately for each band (B, V, K) — do NOT mix bands in a single simulation
  • Each band simulation is independent: packets in the B-band run know nothing about V or K
  • All packets within a band carry the SAME luminosity: \(L_\mathrm{band,total} / N_\mathrm{packets}\)
  • Use dust density (\(\rho_\mathrm{gas} \times f_\mathrm{dust/gas}\)) for optical depth calculations
  • Sample \(\tau = -\ln(\xi)\) with negative sign
  • Grid coordinates: \(x, y, z \in [-L/2, L/2]\)
  • Calculate separate Planck mean opacities for each star/band combination
  • Phase 1: Get V-band working perfectly first, then add B and K bands
  • Track statistics, not individual packet objects (memory efficiency for \(10^6\) packets). Storing \(10^6\) Photon objects will consume ~GB of RAM. Use scalar accumulators (L_escaped, n_absorbed) and numpy arrays for spatial maps instead.

What to Do Next

  1. Complete the Planning Worksheet — units, contracts, hand trace
  2. Fill in constants.py with CGS values
  3. Build your grid and test cell indexing
  4. Implement single-star, constant-\(\kappa\) V-band (Phase 1A)
  5. Validate: \(f_\mathrm{esc} = e^{-\tau}\), energy conservation, convergence
  6. Add multiple stars (Phase 1B), then multiple bands (Phase 2)
  7. Run python run.py validate and python run.py make-figures to check your work

See the full assignment for requirements and the Expectations & Grading for how your work will be evaluated.