User Guide

This comprehensive guide provides detailed tutorials, advanced usage examples, and best practices for using PUFFIN to model externally irradiated protoplanetary disks with photoevaporative winds.

Tutorial 1: Basic 1D Model

Overview

The 1D model computes the radial density structure along the midplane. This is the simplest and fastest option, ideal for quick parameter studies.

Basic Usage

import puffin_disk
import numpy as np
import matplotlib.pyplot as plt

# Define physical parameters
m_star    = 1.0      # Stellar mass (M_sun)
r_d       = 100.0    # Disk outer radius (AU)
sigma_1au = 100.0    # Surface density at 1 AU (g cm^-2)
F_FUV     = 1000.0   # External FUV field (G_0)

# Create and compute model
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV)
r_array, rho_array = model.compute()

# Plot the result
plt.figure(figsize=(8, 6))
plt.loglog(r_array, rho_array, linewidth=2, color='#2E86AB')
plt.xlabel('Radius (AU)', fontsize=14)
plt.ylabel(r'Midplane Density (g cm$^{-3}$)', fontsize=14)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('density_1d.png', dpi=300)
plt.show()

Understanding the Output

The compute() method returns two arrays:

r_array: Radial grid in AU (default: 200 logarithmically-spaced points)
rho_array: Gas density at the midplane in g cm⁻³

The model automatically combines three components:

Disk interior: Power-law surface density with exponential taper
Photoevaporative wind: Spherically diverging outflow
Transition region: Smooth connection using ‘plateau’ prescription (see Keyte & Haworth 2026)

You can also access individual components:

rho_disk    = model.rho_disk      # Disk component only
rho_wind    = model.rho_wind      # Wind component only
rho_plateau = model.rho_plateau   # Transition region

Key Features

Automatic mass-loss rate interpolation: If you don’t specify m_dot, PUFFIN automatically interpolates from the FRIED grid based on your input parameters:

# Mass-loss rate interpolated automatically
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV)
r_array, rho_array = model.compute()

# Check what was used
print(f"Interpolated mass-loss rate: {model.m_dot:.2e} M_sun/yr")

Console output: By default, PUFFIN prints a formatted table showing all model parameters:

▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓ INITIALIZATION ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓

 > Model initialized
 ┌────────────────────┬───────────────┬──────────┐
 │ Parameter          │ Value         │ Unit     │
 ├────────────────────┼───────────────┼──────────┤
 │ Stellar mass       │ 1.00          │ M_sun    │
 │ Disk radius        │ 100.0         │ AU       │
 │ Sigma (1AU)        │ 100.0         │ g cm⁻²   │
 │ FUV field          │ 1.00 × 10³    │ G0       │
 │ N grid cells       │ 200           │          │
 │ Grid size          │ 800.0         │ AU       │
 │ Mass-loss rate     │ 1.23 × 10⁻⁸   │ M_sun/yr │
 │ Gamma              │ 2.45          │          │
 │ p                  │ 0.20          │          │
 │ q                  │ 0.40          │          │
 └────────────────────┴───────────────┴──────────┘

To disable console output:

model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, verbose=False)

Tutorial 2: Basic 2D Model

Overview

The 2D model expands the framework by incorporating the disk’s vertical structure under hydrostatic equilibrium, the photoevaporative wind, and a transition region that smoothly links the disk to the wind.

Basic Usage

import puffin_disk
import numpy as np
import matplotlib.pyplot as plt

# Define physical parameters (same as 1D)
m_star    = 0.5      # Stellar mass (M_sun)
r_d       = 100.0    # Disk outer radius (AU)
sigma_1au = 1000.0   # Surface density at 1 AU (g cm^-2)
F_FUV     = 3000.0   # External FUV field (G_0)

# Create and compute 2D model
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV)
r_array, z_array, rho_array = model.compute()

# Create contour plot
fig, ax = plt.subplots(figsize=(10, 6))
levels = np.arange(-20, -11, 0.2)
contour = ax.contourf(r_array, z_array, np.log10(rho_array), 
                      levels=levels, cmap='Spectral_r', extend='both')
ax.set_xlabel('Radius (AU)', fontsize=14)
ax.set_ylabel('Height (AU)', fontsize=14)
cbar = plt.colorbar(contour, ax=ax, label=r'log$_{10}$ Density (g cm$^{-3}$)')
ax.set_aspect('equal')
plt.tight_layout()
plt.savefig('density_2d.png', dpi=300)
plt.show()

Understanding the Output

The compute() method returns three arrays:

r_array: Radial grid in AU (default: 1000 points)
z_array: Vertical grid in AU (default: 1000 points)
rho_array: 2D gas density array in g cm⁻³ with shape (nz, nr)

Important: The density array uses row ordering where the first dimension is vertical (z) and the second is radial (r). When plotting, use contourf(r_array, z_array, rho_array).

Access model components:

rho_disk = model.rho_disk        # Disk component with hydrostatic structure
rho_wind = model.rho_wind        # Combined wind (spherical + transition)

The typical runtime for a 2D model with n_points=1000 is ~2-3 minutes on a standard laptop.

Tutorial 3: Customizing Model Parameters

Grid Resolution

1D models: Default is 200 points. Higher resolution is recommended for capturing sharper density gradients:

# Standard resolution
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, n_points=200)

# High resolution for detailed structure
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, n_points=1000)

2D models: Default is 1000×1000. Resolution ≥1000 is required for accurate hydrostatic equilibrium:

# Standard resolution
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, n_points=1000)

# High resolution for detailed structure
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, n_points=2000)

If you need a lower-resolution model (e.g., for use as input to thermochemical calculations), first compute the PUFFIN model with at least 1000 grid points, and then regrid it to the desired resolution.

Grid Size

The computational domain extends from 0 to gridsize (in AU). Default is min(8 × r_d, 800) AU. Generally, 8 × r_d is required to capture the full extent of dense photoevaporative winds.

# Modify the gidsize (in AU)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, gridsize=1500)

Mass-Loss Rate

Override automatic interpolation with your own value:

# Specify mass-loss rate directly (M_sun/yr)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, m_dot=1e-8)

# For 2D models, note that PUFFIN applies 2× scaling internally
# to match systematic differences in 2D hydrodynamical simulations
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, m_dot=1e-8)

Taper Parameter (γ)

The exponential tapering parameter gamma controls how steeply the disk truncates. PUFFIN uses an empirical scaling by default, but you can override it:

# Use default empirical prescription
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV)

# Sharper cutoff (higher gamma)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, gamma=4.0)

# Gentler cutoff (lower gamma)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, gamma=2.0)

Transition Region Parameters

Fine-tune the disk-wind transition:

# Standard transition (default)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, p=0.2, q=0.4)

# Slower transition (more extended plateau)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, p=0.1, q=0.2)

# Faster transition (sharper disk-wind interface)
model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, p=0.3, q=0.6)

Temperature Gradient Steepness (2D only)

Control the vertical temperature transition sharpness with parameter k:

# Standard gradient (default)
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, k=1.75)

# Sharper temperature gradient
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, k=3.0)

# Gentler gradient
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, k=1.0)

Hydrostatic Equilibrium Iterations (2D only)

Increase iterations if the model fails to converge:

# Standard iteration count (default: 20)
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, N_ITER=20)

# More iterations for better convergence
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, N_ITER=50)

Most models converge within the default 20 iterations.

Tutorial 4: Parameter Space Exploration

Efficiently explore how disk properties vary across parameter space:

import puffin_disk
import numpy as np
import matplotlib.pyplot as plt

# Define parameter grid
masses = [0.3, 0.7, 2.0]       # M_sun
FUV_fields = [1e2, 1e3, 1e4]   # G_0

# Fixed parameters
r_d = 100.0
sigma_1au = 100.0

# Storage for results
results = {}

# Loop over parameter combinations
for m_star in masses:
    for F_FUV in FUV_fields:
        key = f"M{m_star}_F{F_FUV:.0e}"
        
        # Compute model (use verbose=False to reduce output)
        model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, F_FUV, 
                                        verbose=False)
        r, rho = model.compute()
        
        results[key] = {
            'radius': r,
            'density': rho,
            'm_star': m_star,
            'F_FUV': F_FUV,
            'm_dot': model.m_dot
        }

# Plot family of density profiles
fig, ax = plt.subplots(figsize=(10, 6))

for key, data in results.items():
    label = f"M={data['m_star']} M☉, F={data['F_FUV']:.0e} G₀"
    ax.loglog(data['radius'], data['density'], label=label, linewidth=2)

ax.set_xlabel('Radius (AU)', fontsize=14)
ax.set_ylabel(r'Density (g cm$^{-3}$)', fontsize=14)
ax.legend(fontsize=10)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('parameter_scan.png', dpi=300)
plt.show()

Tutorial 5: Integration with Chemical Models

Using `PUFFIN` density structure as input to `DALI` chemical model

DALI (Bruderer et al. 2012, 2013) requires a specific setup_grid.dat file format that includes not only the gas density structure, but also grid cell boundaries, gas-to-dust ratios, and dust settling prescriptions.

Because PUFFIN computes only the gas density structure, converting PUFFIN outputs to DALI format requires you to specify additional disk properties:

Gas-to-dust ratio in the shielded disk interior
Gas-to-dust ratio in the dust-depleted wind
Dust settling parameters (scale height and grain size distribution)
Physical criterion for the disk-wind boundary (e.g., optical depth threshold, density threshold)

This tutorial demonstrates an example workflow based on the implementation in Keyte & Haworth (2025), where the transition from disk to wind is defined using an arbitrary density threshold n=10⁸ cm⁻³. Note that these properties are science choices that depend on your specific research questions—there are no universal “correct” values.

Step 1: Save PUFFIN Output

First, run PUFFIN and save the density structure:

import puffin_disk
import numpy as np
from scipy.interpolate import RegularGridInterpolator

# Generate 2D model at high resolution
m_star = 1.0
r_d = 100.0
sigma_1au = 1000.0
F_FUV = 1000.0

# Compute at high resolution for accuracy
model = puffin_disk.DiskModel2D(m_star, r_d, sigma_1au, F_FUV, n_points=1000)
r_array_hr, z_array_hr, rho_array_hr = model.compute()

# Convert density to number density
mu = 1.4  # Mean molecular weight
m_p = 1.67262192e-24  # Proton mass in g
n_gas_hr = rho_array_hr / (mu * m_p)

# ========================================================================
# REGRID TO LOWER RESOLUTION FOR DALI
# ========================================================================

# Define lower resolution grid for DALI (200x200)
n_r_dali = 200
n_z_dali = 200

# Create grid
r_array = np.logspace(np.log10(r_array_hr[0]), np.log10(r_array_hr[-1]), n_r_dali)
z_min = 1e-2  # Small offset to avoid log(0)
z_array_log = np.logspace(np.log10(z_min), np.log10(z_array_hr[-1]), n_z_dali - 1)
z_array = np.concatenate([[0.0], z_array_log])  # z needs to start at 0 in DALI

# Create interpolator for high-resolution density
interpolator = RegularGridInterpolator(
    (r_array_hr, z_array_hr), 
    n_gas_hr.T,  # Transpose to (nr, nz) for interpolator
    method='linear',
    bounds_error=False,
    fill_value=1e-30
)

# Interpolate to DALI grid
n_gas = np.zeros((n_z_dali, n_r_dali))
for i in range(n_r_dali):
    for j in range(n_z_dali):
        n_gas[j, i] = interpolator([r_array[i], z_array[j]])

print(f"Regridded from {len(r_array_hr)}x{len(z_array_hr)} to {n_r_dali}x{n_z_dali}")

# ========================================================================

# Define disk-wind boundary using density threshold
# Cells above this threshold are "disk", below are "wind"
n_threshold = 1e8  # cm^-3 (typical disk-wind boundary)

# Flatten arrays for saving
n_r, n_z = len(r_array), len(z_array)
r_flat = []
z_flat = []
n_flat = []

for i in range(n_r):
    for j in range(n_z):
        r_flat.append(r_array[i])
        z_flat.append(z_array[j])
        n_flat.append(n_gas[j, i])

# Save to file
output_file = 'puffin_outputs.dat'
np.savetxt(output_file, 
           np.column_stack([r_flat, z_flat, n_flat]),
           header=f'r(AU) z(AU) n_gas(cm^-3) | Grid: {n_r}x{n_z} | n_threshold={n_threshold:.1e} cm^-3',
           fmt='%.6e')

print(f"PUFFIN output saved to {output_file}")
print(f"Grid dimensions: {n_r} x {n_z}")
print(f"Disk-wind boundary set at n = {n_threshold:.1e} cm^-3")

Step 2: Create DALI Setup File

Now create the DALI setup_grid.dat file. This includes additional physics choices about dust properties:

import numpy as np

# Load PUFFIN output
data = np.loadtxt('puffin_outputs.dat', unpack=True)

# Grid dimensions (200x200 after regridding)
n_r = 200  # Number of radial grid points
n_z = 200  # Number of vertical grid points

# Reshape to 2D grids
rr = np.reshape(data[0], (n_r, n_z))
zz = np.reshape(data[1], (n_r, n_z))
ngas = np.reshape(data[2], (n_r, n_z))

# ========================================================================
# DUST PARAMETERS - These are user choices based on your science goals
# ========================================================================

gd_disk = 100        # Gas-to-dust ratio in disk (standard ISM value)
gd_wind = 1/3e-4     # Gas-to-dust in wind (dust-depleted, only small grains)
chi = 0.2            # Scale height parameter for settling
f = 0.9              # Fraction of large grains at midplane
n_threshold = 1e8    # Density threshold for disk-wind boundary (cm^-3)

# ========================================================================

# Create DALI setup_grid.dat file
outfile = 'setup_grid.dat'

with open(outfile, "w") as f_new:
    # Write header
    f_new.write("# number of cells in r-direction \n")
    f_new.write("%i \n" % (n_r - 1))
    f_new.write("# number of cells in z-direction \n")
    f_new.write("%i \n" % (n_z - 1))
    f_new.write("#  i_r    i_z               ra               rb               za               zb            n_gas      gas-to-dust             f_ls region \n")
    
    # Loop over grid cells
    for i in range(n_r - 1):

        # Since PUFFIN directly provides the gas density structure, we first compute
        # h(r) from the gas density at each radius, then apply the 'standard' DALI
        # dust prescriptions (Bruderer 2013, equations 2 and 3)
        
        # Calculate scale height at this radius
        # (defined where density drops to 1/sqrt(e) of midplane value)
        max_dens = np.max(ngas[i, :])
        try:
            scaleheight_idx = np.where(ngas[i, :] < max_dens * np.exp(-0.5))[0][0]
            scaleheight_au = zz[i, scaleheight_idx]
        except:
            scaleheight_idx = 0
            scaleheight_au = zz[i, 0]
        
        for j in range(n_z - 1):
            
            # Cell boundaries
            ra = rr[i][j]
            rb = rr[i+1][j]
            
            if j == 0:
                za = 0.0
                zb = zz[i][j+1]
            else:
                za = zz[i][j]
                zb = zz[i][j+1]
            
            # Cell midpoints
            r_mid = 0.5 * (ra + rb)
            z_mid = max(0.5 * (za + zb), 1e-10)  # Avoid division by zero
            
            # Current cell density
            n_gas = ngas[i][j]
            
            # ===============================================
            # DISK vs WIND: Use density threshold
            # ===============================================
            
            is_disk = (n_gas > n_threshold)
            
            # ===============================================
            # LARGE GRAIN FRACTION (f_ls) - Dust Settling
            # ===============================================

            # This part is a formulation of Eq. 2 and 3 from Bruderer 2013, using the
            # scale height determined above,
            
            # Scale height angle and polar angle
            h = np.arctan(scaleheight_au / r_mid)
            h = np.maximum(h, 1e-12)
            theta = np.arctan(r_mid / z_mid)  # Polar angle from z-axis
            
            if is_disk:
                # Inside dense disk: calculate settling
                prefactor = f / ((1 - f) * chi)
                term1 = -0.5 * (((np.pi/2) - theta) / (chi * h))**2
                term2 = 0.5 * (((np.pi/2) - theta) / h)**2
                expo = np.exp(term1 + term2)
                
                # f_ls = fraction of large grains
                X_ratio = prefactor * expo
                f_ls = X_ratio / (1.0 + X_ratio)
                f_ls = np.maximum(1e-10, f_ls)
            else:
                # Wind region: no large grains
                f_ls = 1e-10
            
            # ===============================================
            # GAS-TO-DUST RATIO
            # ===============================================
            
            if is_disk:
                gtd = gd_disk  # Dense disk
            else:
                gtd = gd_wind  # Tenuous wind (dust-depleted)
            
            # Region identifier (typically 1 for entire disk)
            region = 1
            
            # Write row
            f_new.write("%6i %6i %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %15.10e %6i \n" 
                       % (i+1, j+1, ra, rb, za, zb, n_gas, gtd, f_ls, region))

print(f"DALI setup file created: {outfile}")
print(f"Grid dimensions: {n_r-1} x {n_z-1} cells")
print(f"Disk-wind boundary: n > {n_threshold:.1e} cm^-3")
print("Ready to run DALI with this grid")

Understanding the DALI Format

The setup_grid.dat file has the following structure:

Header:

Number of radial cells
Number of vertical cells
Column descriptions

Data columns (for each cell):

i_r: Radial cell index
i_z: Vertical cell index
ra: Inner radial boundary (AU)
rb: Outer radial boundary (AU)
za: Lower vertical boundary (AU)
zb: Upper vertical boundary (AU)
n_gas: Gas number density (cm⁻³)
gas-to-dust: Gas-to-dust mass ratio
f_ls: Large grain fraction (0-1)
region: Region identifier (typically 1)

Key Physics Choices

When creating DALI input, you must specify dust parameters based on your science goals:

Gas-to-Dust Ratio:

gd_disk (e.g., 100): Standard ISM value for shielded disk interior
gd_wind (e.g., 3000 or higher): Dust-depleted wind where grains have settled/been lost (eg. Facchini et al. 2016)
Disk/wind boundary based on eg. density threshold (see Keyte & Haworth 2025) or optical depth τ = 1 threshold (see Keyte & Haworth 2026)

Dust Settling Parameters:

f (e.g., 0.9): Fraction of large grains at the midplane
chi (e.g., 0.2): Ratio of dust scale height to gas scale height
These control vertical stratification of grain sizes

The f_ls (large grain fraction) calculation implements a settling prescription where:

Large grains settle toward the midplane
Small grains remain well-mixed
The vertical distribution depends on the local scale height angle

Important Notes:

Temperature fields: As always, do NOT use PUFFIN temperatures. DALI will compute self-consistent temperatures through radiative transfer.
Dust choices are yours: The gas-to-dust ratios and settling parameters shown here are examples. Adjust based on your specific disk model and science questions.

Verification

After creating the setup file, verify it looks correct:

# Quick visual check
import matplotlib.pyplot as plt

# Load the created setup file
setup_data = np.loadtxt('setup_grid.dat', skiprows=4, unpack=True)

# Reshape to 2D (note: n_r-1 and n_z-1 cells)
n_cells_r = n_r - 1
n_cells_z = n_z - 1

setup_r = np.reshape(setup_data[2], (n_cells_r, n_cells_z))  # rb values
setup_z = np.reshape(setup_data[4], (n_cells_z, n_cells_z))  # zb values  
setup_ngas = np.reshape(setup_data[6], (n_cells_r, n_cells_z))
setup_fls = np.reshape(setup_data[8], (n_cells_r, n_cells_z))

fig, axes = plt.subplots(1, 2, figsize=(15, 5))

# Plot density
levels = np.arange(2, 14, 0.2)
im1 = axes[0].contourf(setup_r, setup_z, np.log10(setup_ngas), 
                       levels=levels, cmap='Spectral_r', extend='both')
axes[0].set_xlabel('Radius (AU)')
axes[0].set_ylabel('Height (AU)')
axes[0].set_title('Gas Number Density')
plt.colorbar(im1, ax=axes[0], label='log₁₀ n (cm⁻³)')

# Plot settling
im2 = axes[1].contourf(setup_r, setup_z, setup_fls,
                       levels=np.arange(0, 1, 0.1), cmap='Spectral_r', extend='both')
axes[1].set_xlabel('Radius (AU)')
axes[1].set_ylabel('Height (AU)')
axes[1].set_title('Large Grain Fraction')
plt.colorbar(im2, ax=axes[1], label='f_ls')

plt.tight_layout()
plt.savefig('dali_setup_verification.png', dpi=200)
plt.show()

Best Practices

Temperature and Radition Fields: Critical Warning

⚠️ NEVER use PUFFIN temperature/radiation fields directly for chemistry or radiative transfer.

The temperatures computed by PUFFIN are used only to establish the density structure through hydrostatic equilibrium. They are not self-consistent with the radiation field and should not be trusted for any downstream analysis.

Always:

Use PUFFIN to generate density structures only
Input these densities into proper radiative transfer codes (DALI, RADMC-3D, ProDiMo)
Use the self-consistent temperatures and radiation fields from radiative transfer for chemistry

Troubleshooting

Common Issues

Issue: “Disk mass < 1e-5 M_sun” abort message

Cause: Surface density too low or disk radius too small

Solution: Increase sigma_1au or r_d

Issue: 2D model produces unusual structure

Cause: Grid resolution too low or parameters produce extreme conditions

Solutions:

Increase n_points to >>1000
Increase N_ITER for better convergence
Check parameters are within validated range
Vary k parameter for gentler/steeper temperature gradients
Manually specify gamma, p, and q to adjust disk-wind transition

Issue: Wind densities seem too high/low

Cause: Mass-loss rate interpolation or manual specification

Solutions:

Check interpolated m_dot value: print(model.m_dot)
Manually specify m_dot if needed
For 2D, remember the 2× internal scaling factor
Verify you’re within FRIED grid coverage

Issue: Unexpected density jumps or discontinuities

Cause: Insufficient grid resolution or numerical artifacts

Solutions:

Increase n_points
Adjust gamma parameter for smoother taper
Modify p and q for gentler transition

Issue: Results don’t match hydrodynamical simulations exactly

Expected: PUFFIN is a parametric approximation with typical factor-of-a-few accuracy

When to worry: Deviations >5× suggest:

Parameters outside validated range
Unusual physics not captured by parametric model
Need for full hydrodynamical simulation

Further Resources

Code Repository

GitHub: https://github.com/lukekeyte/PUFFIN
Report bugs, request features, or contribute improvements

Documentation

Full documentation: https://puffin.readthedocs.io
Physical model details: See “Physical Model” section
API reference: See source code docstrings

Key References

PUFFIN framework:

Keyte & Haworth (2026), arXiv:2602.02011

FRIED hydrodynamical grid:

Haworth et al. (2018), MNRAS 481, 452
Haworth et al. (2023), MNRAS 526, 4315

Getting Help

Check the documentation: Most questions answered in Physical Model and FAQ sections
Email support: l.keyte@qmul.ac.uk for technical questions
GitHub Issues: For bug reports and feature requests

User Guide

Tutorial 1: Basic 1D Model

Overview

Basic Usage

Understanding the Output

Key Features

Tutorial 2: Basic 2D Model

Overview

Basic Usage

Understanding the Output

Tutorial 3: Customizing Model Parameters

Grid Resolution

Grid Size

Mass-Loss Rate

Taper Parameter (γ)

Transition Region Parameters

Temperature Gradient Steepness (2D only)

Hydrostatic Equilibrium Iterations (2D only)

Tutorial 4: Parameter Space Exploration

Tutorial 5: Integration with Chemical Models

Using PUFFIN density structure as input to DALI chemical model

Step 1: Save PUFFIN Output

Step 2: Create DALI Setup File

Understanding the DALI Format

Key Physics Choices

Verification

Best Practices

Temperature and Radition Fields: Critical Warning

Troubleshooting

Common Issues

Further Resources

Code Repository

Documentation

Key References

Getting Help

Using `PUFFIN` density structure as input to `DALI` chemical model