Frequently Asked Questions

Installation & Setup

How do I install PUFFIN?

PUFFIN can be installed directly from PyPI using pip:

pip install puffin_disk

For development or local installation from source:

git clone https://github.com/lukekeyte/PUFFIN.git
cd PUFFIN
pip install -e .

What are the system requirements?

PUFFIN requires:

Python 3.8 or higher
numpy >= 1.20.0
scipy >= 1.7.0
matplotlib >= 3.3.0

The package has been fully tested on Mac/Linux, with limited testing on Windows.

I’m getting import errors. What’s wrong?

First, ensure all dependencies are installed:

pip install numpy scipy matplotlib

If you’re still experiencing issues, try updating your packages:

pip install --upgrade numpy scipy matplotlib

For development installations, make sure you’re in the correct directory and using the -e flag.

Model Selection & Usage

How do I choose the grid resolution (`n_points`)?

Minimum requirement: 200 grid points for 1D models and 1000x1000 grid points for 2D models (enforced)

Recommended for most applications: >1000 grid points

Increasing resolution improves accuracy but increases runtime, especially for 2D models.

NOTE: If you need a lower-resolution 2D 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 (see below).

What grid size should I use?

The default grid size is min(8 × r_d, 800) AU, which captures the disk and immediate wind region.

Increase gridsize when:

Studying extended photoevaporative flows beyond 8×r_d
The wind component dominates at large radii
You need to match observations covering large angular extents

Decrease gridsize when:

Focusing on the inner disk structure
Wind structure is minimally extended

Parameter Ranges & Validation

What happens if I use parameters outside the validated range?

PUFFIN will raise a ValueError if you attempt to use parameters outside these ranges:

Stellar mass: 0.3 to 3.0 M☉
Disk radius: 20 to 150 AU
Surface density: 10 to 10,000 g cm⁻²
FUV field: 100 to 100,000 G₀

These limits reflect the calibration against the FRIED grid. Using parameters outside these ranges may produce physically unrealistic results.

Can I extrapolate beyond the validated ranges?

While you could technically modify the code to extrapolate, this is strongly discouraged unless you:

Carefully validate output structures against reasonable expectations
Have independent constraints (e.g., hydrodynamical simulations) to compare against
Clearly document that you’re operating outside the validated regime

The parametric prescriptions were calibrated specifically for the validated parameter space and may not generalize well.

Mass-Loss Rates

How is the mass-loss rate determined?

If you don’t provide m_dot explicitly, PUFFIN interpolates from the FRIED grid (Haworth et al. 2018, 2023) using 4D linear interpolation across stellar mass, disk radius, surface density, and FUV field strength.

For 2D models, the interpolated value is scaled by a factor of 2 to account for systematic differences between 1D and 2D hydrodynamical simulations.

Can I specify my own mass-loss rate?

Yes! Simply pass the m_dot parameter (in M☉ yr⁻¹):

model = puffin_disk.DiskModel1D(m_star, r_d, sigma_1au, FFUV_G0, m_dot=1e-7)

NOTE: If using a user-defined mass-loss rate, take care to ensure it is physically justified given the model input parameters.

Output Interpretation & Usage

What do the output arrays represent?

1D Model returns:

radius: Radial grid in AU (shape: [n_points])
density: Total midplane density in g cm⁻³ (shape: [n_points])

2D Model returns:

r_array: Radial grid in AU (shape: [nr])
z_array: Vertical grid in AU (shape: [nz])
density: Total gas density in g cm⁻³ (shape: [nz, nr])

Both models also store component densities (disk, wind, transition) as instance attributes for diagnostic purposes.

Can I use the temperature arrays for chemical modeling?

No. This is explicitly stated in the documentation but bears repeating:

The temperature profiles computed during model construction are used only to establish the density structure through hydrostatic equilibrium. They do not represent self-consistent temperature fields for a given density distribution.

For chemical modeling or synthetic observations, you must:

Export PUFFIN density structures
Pass them to a dedicated radiative transfer code (e.g., DALI, RADMC-3D, ProDiMo)
Use the self-consistent temperature and radiation fields from that code

How do I convert PUFFIN output for use with other codes?

Most radiative transfer codes accept density grids in similar formats. An detailed example using PUFFIN with the DALI code is provided in the User Guide.

Consult your radiative transfer code’s documentation for specific input requirements (units, grid orientation, format).

What units should I use?

PUFFIN uses the following units internally:

Length: AU (radii) and cm (internal calculations)
Density: g cm⁻³
Temperature: K
Mass-loss rate: M☉ yr⁻¹
FUV field: G₀ (Habing units, 1.6 × 10⁻³ erg cm⁻² s⁻¹)

All outputs use these same units, so you may need to convert for specific applications.

Troubleshooting & Performance

My 2D model is taking a very long time to run. Is this normal?

Typical runtimes depend on resolution:

1000 grid points: ~2-3 minutes
1500 grid points: ~5 minutes
2000 grid points: ~10 minutes

If your model is taking significantly longer:

Check your grid resolution: Is n_points unnecessarily high?
Verify convergence: Look for console messages about iteration progress

The model returned an error message. What does it mean?

Common error messages:

“Stellar mass must be between 0.3 to 3.0 M_sun”: Your m_star parameter is outside the validated range.

“Need at least 200 grid points”: Increase n_points to ensure numerical accuracy.

“Grid size must be larger than disk radius”: Your gridsize is too small; increase it or use the default.

“ABORTED: Disk mass < 1e-5 M_sun”: The disk is too low-mass for the model to handle reliably. Increase sigma_1au or r_d.

“Hydrostatic equilibrium did not converge”: The iterative solver failed. Try:

Increasing N_ITER (default: 20)
Adjusting the temperature gradient parameter k (default: 1.75)
Checking for extreme parameter combinations

How do I access individual density components?

After running the model, you can access:

1D model:

model.rho_disk: Disk component
model.rho_wind: Wind component
model.rho_plateau: Plateau transition component

2D model:

model.rho_disk: Disk component with hydrostatic equilibrium
model.rho_wind_sph: Spherical wind component
model.rho_bowl: Disk-wind transition (“halo”) component
model.rho_wind: Combined wind (max of spherical and bowl)

These are useful for understanding which component dominates in different regions.

Physical Model & Limitations

What physics is included in PUFFIN?

PUFFIN models:

Power-law disk surface density with exponential tapering
Vertical hydrostatic equilibrium (2D model)
FUV heating in photodissociation region (PDR)
Spherically diverging photoevaporative wind
Smooth disk-wind transition regions

What physics is NOT included?

PUFFIN does not model:

Internal photoevaporation from the central star
External EUV radiation
Magnetic fields or MHD effects
Time evolution or dynamical changes
Dust dynamics (assumes well-mixed gas and dust for temperature structure only)
Detailed chemical networks (use radiative transfer codes for this)
Lots of other things!

If these effects are important for your system, you’ll need hydrodynamical simulations or specialized codes.

How accurate are the density structures?

PUFFIN reproduces hydrodynamical simulations with typical factor-of-a-few accuracy across ≥98% of the validated parameter space. Larger deviations can occur:

Near sharp density gradients in the disk-wind interface
In regions with strong FUV irradiation gradients
At the edges of the validated parameter space

For most chemical modeling applications, this accuracy is sufficient given other uncertainties in disk properties.

What is the `gamma` parameter and when should I change it?

The gamma parameter controls the steepness of the exponential disk taper:

Σ(r) ∝ exp[−(r/r_d)^γ]

By default, PUFFIN calculates γ using:

γ = 0.77 × (M_*/M☉)^0.10 × (r_d/AU)^0.78 × (F_FUV/100 G₀)^0.07

where the coefficients were obtained using MCMC fitting to the FRIED grid of hydrodynamical simulations.

PUFFIN also allows for user-defined values of gamma.

Best Practices & Recommendations

What’s the recommended workflow for a chemical modeling project?

Define your disk parameters based on observations or theoretical expectations
Generate the density structure using PUFFIN 2D model
Validate the structure by plotting and checking for physical reasonableness
Export to radiative transfer (e.g., DALI, ProDiMo) for self-consistent temperature/radiation
Run chemistry on the full 2D structure with proper temperature/radiation fields
Generate observables and compare with data

How do I cite PUFFIN in my research?

If you use PUFFIN, please cite:

Primary reference:

Keyte & Haworth (2026), arXiv:2602.02011
"A parametric model for externally irradiated protoplanetary disks with photoevaporative winds"

For mass-loss rates from FRIED:

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

Where can I get help or report bugs?

Documentation: https://puffin.readthedocs.io
GitHub Issues: https://github.com/lukekeyte/PUFFIN
Email: l.keyte@qmul.ac.uk

When reporting issues, please include:

Your operating system and Python version
The complete error message
A minimal example that reproduces the problem
Your parameter values

Advanced Usage

Can I modify the source code for my specific application?

Yes! PUFFIN is open-source under the MIT License. You’re welcome to:

Modify prescriptions for your specific needs
Add new physical components
Integrate with your own codes

If you make improvements that might benefit the community, please consider contributing them back via a pull request.

How do I run PUFFIN in batch mode for parameter surveys?

Suppress verbose output and wrap in a loop:

import puffin_disk
import numpy as np

# Parameter ranges
masses = np.linspace(0.3, 3.0, 10)
fuv_fields = np.logspace(2, 5, 10)

results = []
for m_star in masses:
    for FFUV_G0 in fuv_fields:
        model = puffin_disk.DiskModel1D(
            m_star, r_d=100, sigma_1au=1000, FFUV_G0=FFUV_G0,
            verbose=False
        )
        r, rho = model.compute()
        results.append((m_star, FFUV_G0, r, rho))

How do I interpolate PUFFIN output onto a different grid?

Use scipy interpolation tools:

from scipy.interpolate import RegularGridInterpolator

# Define lower resolution grid 200x200)
n_r_new = 200
n_z_new = 200

# Create grid
r_array = np.logspace(np.log10(r_array_highres[0]), np.log10(r_array_highres[-1]), n_r_new)
z_min = 1e-2  # Small offset to avoid log(0)
z_array_log = np.logspace(np.log10(z_min), np.log10(z_array_highres[-1]), n_z_new - 1)
z_array = np.concatenate([[0.0], z_array_log])  # z needs to start at 0

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

# Interpolate to DALI grid
rho_new = np.zeros((n_z_new, n_r_new))
for i in range(n_r_new):
    for j in range(n_z_new):
        rho_new[j, i] = interpolator([r_array[i], z_array[j]])

Can I access the FRIED grid data directly?

Yes! The FRIED data is available at https://github.com/thaworth-qmul/FRIEDgrid