Relaxation Module

For examples of how to use these functions, see Quickstart.

relax_wrapper.py:: Contains the class for wrapping the relaxation code. Handles loading, saving, running, and converging boundary conditions for the relaxation simulations. The goal is find the solution of interest to the end-user and automate the choices one must make to be intelligently picked and consistent, i.e., the boundary conditions and ramping thru parameter space. Handling of the solution is done by the other major class wind_solution found in wind_ae/wrapper/wrapper_utils/windsoln.py.

Loading and interacting with wind solutions

wind_simulation.load_planet(csv_file, calc_postfacto=True, expedite_postfacto=False, name='Loaded Planet', print_atmo=True, print_warnings=True)

Loads a planet solution file into inputs/guess.inp as a new guess and updates input parameter folders. Ramping to a new planet with ramp_to() will rewrite parameter folders for the desired new planet, but the guess will remain the same.

Parameters:

csv_file (str) – Path to the csv file containing the planet solution.
calc_postfacto (bool, optional) – If True, calculates post-facto variables (e.g., photoionization heating, etc.). Defaults to True.
name (str, optional) – Name of the planet, used for identification. Defaults to ‘Loaded Planet’.
print_atmo (bool, optional) – If True, prints the atmosphere composition and mass fractions. Defaults to True.
print_warnings (bool, optional) – If True, prints warnings from the wind_solution analysis. Defaults to True.

Returns:

None

wind_simulation.run_wind(expedite=False, calc_postfacto=True, verbose=False, retry=False)

Runs the C relaxation code. If successful, returns 0, saves solution in saves/windsoln.csv, and copies this solution to be the next guess in inputs/guess.inp.

Parameters:

expedite (bool, optional) – If True, will also run outward integration. Will return 4 if relaxation code has solved, but outward integration fails (does not affect profile or mass loss accuracy). Defaults to False.
calc_postfacto (bool, optional) – If True, performs post-facto calculations (e.g., heating/cooling, number density, etc.). Defaults to True.
verbose (bool, optional) – If True, prints verbose output. Defaults to False.
retry (bool, optional) – If True, retries on analysis errors. Defaults to False.

Returns:

Status code.: 0 - Both relaxation and outward integration have succeeded (if expedite=False), or relaxation has succeeded (if expedite=True). 4 - Relaxation has succeeded, but outward integration has failed. 1 - Relaxation has failed. 2 - Error in copying saves/windsoln.csv to inputs/guess.inp (rare). 3 - Analysis error when executing wind_solution() to load current solution into windsoln object.

Return type:

int

wind_simulation.load_uservars(csv_file=None, expedite=False)

Loads all user variables for plotting and analysis. Does not rewrite input parameter files, so can be run concurrently with simulations.

Parameters:

csv_file (str, optional) – Path to the csv file containing the planet solution. Defaults to None.
expedite (bool, optional) – If True, will not calculate expensive post-facto variables (e.g., photoionization heating, etc.) to save time. Defaults to False.

Returns:

None

wind_simulation.generate_rate_coeffs()

Generates rate coefficients for secondary ionization and populates array in src/rate_coeffs.h. Rate coefficients are generated from interpolation over the Dere (2007) table.

Returns:: None

wind_simulation.load_spectrum(generate=False, wl_norm=1e-07, print_warnings=False)

Loads a high resolution spectrum and, if generate=True, smooths it, keeping the flux at ionization edges for the species in the current windsoln accurate, and conserving flux locally in the vicinity of spectral peaks. The spectrum specified in the windsoln is the one loaded and smoothed. One can change the spectrum using ramp_to_user_spectrum() to ramp to a user-input spectrum or change the wavelength range and integrated flux via ramp_spectrum(). To change the total flux but maintain spectrum shape and range, use ramp_flux().

Parameters:

generate (bool, optional) – If True, generates a new smoothed spectrum stored in inputs/spectrum.inp. Defaults to False.
wl_norm (float, optional) – Converts nm to cm. Should not need to be changed. Defaults to 1e-7.
print_warnings (bool, optional) – If True, prints warnings from McAstro/planets/insolation/glq_spectrum.csv. Defaults to False.

Returns:

None

wind_simulation.save_planet(filename=None, filepath='', overwrite=False, polish=False)

Copies wind_ae/saves/windsoln.csv to the desired folder. Note: It is not necessary to define a folder. All files will be saved in the saves/ repo by default.

Parameters:

filename (str, optional) – Name and/or path+name. .csv not required. Defaults to None.
filepath (str, optional) – Path to the folder where the file will be saved. Defaults to ‘’.
overwrite (bool, optional) – If True, will overwrite existing file with the same name. Defaults to False.
polish (bool, optional) – If True, will polish the boundary conditions before saving. Defaults to False.

Returns:

None

wind_simulation.easy_output_file(outputs=['v', 'T', 'rho'], output_file='/home/docs/checkouts/readthedocs.org/user_builds/wind-ae/envs/latest/lib/python3.13/site-packages/wind_ae/saves/simplified_outputs/output.dat', comments='', overwrite=False)

Writes desired solution variables to a csv file that can be easily shared.

Parameters:

outputs (list of str, optional) – Desired solution columns. ‘r’ not required. Defaults to [‘v’, ‘T’, ‘rho’].
output_file (str, optional) – Output file name. Can include path. Defaults to ‘output.dat’.
comments (str, optional) – Additional comments to include in the header. Defaults to ‘’.
overwrite (bool, optional) – If True, will overwrite existing file. Defaults to False.

Returns:

None

Ramping functions

wind_simulation.ramp_to(system=None, intermediate_converge_bcs=False, final_polish=False, integrate_out=True, static_bcs=False, make_plot=False)

Ramps the current wind solution to a new system. Wrapper for ramp_Ftot(), ramp_grav(), ramp_star().

Parameters:

system (system, optional) – The target system to ramp towards. Syntax: system = system(Mp, Rp, Mstar, a, Ftot, Lstar) in cgs. Defaults to None.
intermediate_converge_bcs (bool, optional) – If True, converges boundary conditions during and after each ramp function. Increases runtime, but may improve likelihood of convergence for sensitive solutions or solutions distant from the initial guess. Defaults to False.
final_polish (bool, optional) – If True, converges boundary conditions (i.e. polish_bcs()) after all ramping is done. Defaults to False.
integrate_out (bool, optional) – If True, integrates final solution out past sonic point to Coriolis radius. Defaults to True.
static_bcs (bool, optional) – If True, keeps current boundary conditions static during ramping. Defaults to False.
make_plot (bool, optional) – If True, generates velocity, temperature, density, and ionization fraction plots of steps in the ramping process. Defaults to False.

Returns:

Status code. 0 for success, other values for failure modes.

Return type:

int

wind_simulation.ramp_Ftot(F_goal, converge_bcs=False, expedite=False, integrate_out=True, static_bcs=False, make_plot=True)

Ramps the current wind solution to a new Ftot in ergs/s/cm2.

Parameters:

F_goal (float) – The target Ftot value to ramp towards.
converge_bcs (bool, optional) – If True, converges boundary conditions after ramping Ftot. Defaults to False.
expedite (bool, optional) – If True, expedites the ramping process, but may reduce likelihood of convergence for sensitive solutions. Defaults to False.
integrate_out (bool, optional) – If True, integrates out past the sonic point. Defaults to True.
static_bcs (bool, optional) – If True, uses static boundary conditions. Defaults to False.
make_plot (bool, optional) – If True, produces a four panel plot of velocity, temperature, density, and ionization fraction for steps in the ramping process. Defaults to True.

Returns:

Status code from ramp_var.

Return type:

int

wind_simulation.ramp_var(var, var_end, var_class=None, delta=0.02, delta_additive=False, converge_bcs=False, make_plot=True, expedite=False, integrate_out=True, static_bcs=False)

Ramps a variable in the system or physics class to a target value with adaptive stepsizes.

Parameters:

var (str) – Name of the variable to ramp. Options: var_class=’system’: var = ‘Ftot’, ‘Mp’, ‘Rp’, ‘Mstar’, ‘semimajor’ var_class=’physics’: var = ‘HX’, ‘species_list’, ‘atomic_masses’, ‘Z’, ‘ne’, ‘molec_adjust’
var_end (float) – Target value for the variable.
var_class (str, optional) – Class of the variable (‘system’ or ‘physics’). Defaults to None.
delta (float, optional) – Default stepsize. Defaults to 0.02.
delta_additive (bool, optional) – If True, uses additive steps. Defaults to False.
converge_bcs (bool, optional) – If True, converges boundary conditions during ramping. Defaults to False.
make_plot (bool, optional) – If True, plots ramping progress. Defaults to True.
expedite (bool, optional) – If True, expedites ramping. Defaults to False.
integrate_out (bool, optional) – If True, integrates out after ramping. Defaults to True.
static_bcs (bool, optional) – If True, uses static boundary conditions. Defaults to False.

Returns:

Status code. 0 for success, other values for failure modes.

Return type:

int

wind_simulation.ramp_grav(system, delta=0.02, converge_bcs=False, make_plot=True, expedite=False, integrate_out=True, static_bcs=False)

Ramps planet mass and radius along lines of constant surface gravity with adaptive stepsizes.

Parameters:

system (system) – System object with Mp, Rp, Mstar, semimajor, Ftot, Lstar in cgs units.
delta (float, optional) – Default stepsize. Defaults to 0.02.
converge_bcs (bool, optional) – If True, runs polish_bcs(), which can be costly but may improve likelihood of convergence. Defaults to False.
make_plot (bool, optional) – If True, plots temperature, velocity, density, and ionization fraction profiles at steps throughout the ramp. Defaults to True.
expedite (bool, optional) – If True, skips integrating out to improve speed. Defaults to False.
integrate_out (bool, optional) – If True, requires integrate out to R_cori at the end of the ramp. Defaults to True.
static_bcs (bool, optional) – If True, skips running ramp_base_bcs() when ramper is having difficulty converging. Useful when not wanting to use default BCs. Defaults to False.

Returns:

Status code.: 0 - Successfully ramped (and integrated out and/or converged BCs) 101 - Failed to ramp to desired Mp, Rp 1 - (Unusual) Successfully ramped to Mp, Rp then encountered trouble solving. Last working solution should be sufficient. 4 - Successfully ramped to Mp, Rp, but unable to integrate out. Try running converge_Rmax() or integrate_out(). 5 - Successfully ramped to Mp, Rp, but unable to polish all BCs. See printout for more details.

Return type:

int

wind_simulation.ramp_star(system, delta=0.02, converge_bcs=True, make_plot=True, expedite=False, integrate_out=True, static_bcs=False)

Ramps stellar mass (Mstar), semimajor axis (semimajor), and stellar bolometric luminosity (Lstar). If Mstar and semimajor axis change, ramps linearly along Hill radius rate of change. If semimajor axis and Lstar change, ramps linearly along optical flux (essentially along skin temperature since T_skin propto Fopt) rate of change.

Parameters:

system (system) – System object that takes Mp, Rp, Mstar, semimajor, Ftot, Lstar in cgs units.
delta (float, optional) – Default stepsize. Defaults to 0.02. (Unlikely to need to change)
converge_bcs (bool, optional) – If True, runs polish_bcs(), which can be costly but may improve likelihood of convergence. Defaults to True.
make_plot (bool, optional) – If True, plots temperature, velocity, density, and ionization fraction profiles at steps throughout the ramp. Defaults to True.
expedite (bool, optional) – If True, skips integrating out to improve speed. Defaults to False.
integrate_out (bool, optional) – If True, requires integrate out to R_cori at the end of the ramp. Defaults to True. Note: R_cori will still need to be computed by running converge_Rmax().
static_bcs (bool, optional) – If True, skips running ramp_base_bcs() when ramper is having difficulty converging. Useful when not wanting to use default BCs. Defaults to False.

Returns:

Status code.: 0 - Successfully ramped (and integrated out and/or converged BCs) 101 - Failed to ramp to desired Mstar, Lstar, semimajor 1 - Successfully ramped to Mstar, Lstar, semimajor then encountered trouble solving. Last working solution should be sufficient (self.path+’saves/windsoln.csv’). 4 - Successfully ramped to Mstar, Lstar, semimajor but unable to integrate out. Try running converge_Rmax() or integrate_out(). 5 - Successfully ramped to Mstar, Lstar, semimajor but unable to polish all BCs. See printout for more details.

Return type:

int

Metals functions

wind_simulation.calc_metallicity(metals_list=None, Z=1)

Returns an array of the Lodders (2009) Z x solar metallicity MASS fractions for the metals in metals_list. Mass fractions must sum to 1, so the difference from 1 is added to the hydrogen mass fraction (always at index 0 in the array).

Parameters:

metals_list (list of str) – List of metals, e.g., [‘H I’, ‘He I’, ‘Mg II’].
Z (float) – Multiple of solar metallicity.

Returns:

Array of mass fractions for each species in metals_list.

Return type:

np.ndarray

wind_simulation.add_metals(desired_species_list, Z=1, custom_mfs=[], Ncol_sp_tot=0, integrate_out=True, converge_Ncol_sp=True)

Adds NEW metals to existing solutions. To ramp mass fractions or metallicity of species already in the solution, use ramp_metallicity(). Skips adding if species is already present. This function will converge to a self-consistent Ncol_sp for each species at the final step if converge_Ncol_sp=True.

Parameters:

desired_species_list (list of str) – List of all species desired in the solution. Give element + ionization state (e.g., ‘HI’, ‘h1’, ‘h I’ for neutral hydrogen).
Z (float, optional) – Metallicity in units of solar metallicity. Defaults to 1.
custom_mfs (list or array, optional) – User-specified mass fractions. Defaults to [].
Ncol_sp_tot (int, optional) – Total sonic point column density across all species. Defaults to 0.
integrate_out (bool, optional) – If True, integrates past sonic point to Coriolis radius. Defaults to True.
converge_Ncol_sp (bool, optional) – If True, converges Ncol_sp for each species. Defaults to True.

Returns:

Status code. 0 for success, other values for failure modes.

Return type:

int

wind_simulation.remove_metals(remove_species_list, run_wind=True)

Removes a metal from the list. Note: Helium will take a significant amount of time to remove, so it is advisable to ramp directly to your desired solution from an existing pure-H solution.

Parameters:

remove_species_list (list of str) – List of species to remove.
run_wind (bool, optional) – If True, runs wind relaxation after removal. Defaults to True.

Returns:

None

wind_simulation.ramp_metallicity(Z=1, custom_mfs=[], static_bcs=False, integrate_out=True, converge_Ncol_sp=True, polish_bcs=False)

Ramps up the metallicity of the species present in the simulation. Can ramp in multiples of solar Z or set custom mass fractions.

Parameters:

Z (float, optional) – Multiples of Lodders (2008) solar metallicity. Defaults to 1.
custom_mfs (list of float, optional) – Custom mass fractions. If empty, defaults to Z. If not empty, uses custom_mfs. Defaults to [].
static_bcs (bool, optional) – If True, does not ramp lower boundary conditions (may affect convergence). Defaults to False.
integrate_out (bool, optional) – If True, integrates outwards past the sonic point to the Coriolis radius. Defaults to True.
converge_Ncol_sp (bool, optional) – If True, attempts to converge the Ncol_sp parameter to self-consistency after solving (if integrate_out=False). Defaults to True.
polish_bcs (bool, optional) – If True, polishes all boundary conditions after solving. Defaults to False.

Returns:

Status code. 0 for success, other values for failure modes.

Return type:

int

Boundary Conditions

Polishing (Converging to self-consistent BCs)

wind_simulation.polish_bcs(converge_Rmax=True, static_bcs=None, user_override_press=False, base_press=1, bolo_on=True)

Polishes upper and lower boundary conditions to self-consistency.

Parameters:

converge_Rmax (bool, optional) – If True, converges Rmax (more costly). Defaults to True.
static_bcs (bool, optional) – If True, skips running ramp_base_bcs() and maintains the current boundary conditions. Defaults to None (i.e., the value of self.static_bcs).
user_override_press (bool, optional) – If True, allows user to override pressure boundary conditions. Defaults to False.
base_press (float, optional) – Base pressure in microbars (barye) to use for boundary conditions. Defaults to 1.

Returns:

0 if successfully polished all BCs, 5 if failed at polishing one or more BCs. See printout for more details. If “Molec. to atomic transition” has failed, plot energy_plot() to assess whether the current solution is satisfactory.

Return type:

int

wind_simulation.ramp_base_bcs(user_override_press=False, base_press=1, static_bcs=None, Kappa_opt=0.004, Kappa_IR=0.01, molec_adjust=None, adiabat=False, polish=False, intermediate=False, tolerance=0.1)

Ramps Rmin, mass density at Rmin (rho_rmin), and temperature at Rmin (T_rmin) to the physical values computed in base_bcs(). Default is the values at base_press = 1 microbar. If the loaded solution base BC is at higher pressure, base_press will stay that value in subsequent iterations.

Parameters:

base_press (int or float, optional) – Desired pressure at base of simulation in microbars. Defaults to 1. Default will be overridden if loaded solution has different base press, unless user_override_press=True.
user_override_press (bool, optional) – To lower or raise pressure of base, set True. Defaults to False.
static_bcs (bool, optional) – If True, skips running ramp_base_bcs() and maintains the current boundary conditions. Defaults to None (i.e., the value of self.static_bcs).
Kappa_opt (float, optional) – Optical opacity. Sets bolometric heating/cooling in molecular region below wind. Defaults to 4e-3.
Kappa_IR (float, optional) – IR opacity. Defaults to 1e-2.
molec_adjust (float, optional) – Dimensionless mean molecular weight. Mu of H2 is 2.3*mH. Defaults to self.windsoln.molec_adjust if =None. If molec_adjust <= 0, turns off molecular layer.
adiabat (bool, optional) – If True, computes base BCs assuming atmosphere is adiabatic below wind. Defaults to False.
polish (bool, optional) – If True, skips ramping T when self.windsoln.bolo_heat_cool=0 to avoid costly iteration. Defaults to False.
intermediate (bool, optional) – If True, returns early if average difference is below tolerance. Defaults to False.
tolerance (float, optional) – Tolerance for intermediate ramping. Defaults to 0.1.

Returns:

0 if successfully ramped base boundary conditions, 1 if failed.

Return type:

int

wind_simulation.converge_mol_atomic_transition(width_factor=0.0, _called_in_ramp_bcs=False, polish=False, bolo_on=True)

_Formerly run_isotherm()_

Ramps the complementary error function that governs the transition from the molecular to atomic regions in the atmosphere. Both the mean molecular weight and bolometric heating/cooling are transitioned using the same error function.

If the bolometric heating/cooling flag is off (sim.windsoln.bolo_heat_cool=0.0) there will be no bolometric heating/cooling, but the mean molecular weight will still transition from molecular to atomic unless sim.windsoln.molec_adjust = 0.0. To turn off the bolometric heating/cooling run sim.turn_off_bolo().

Criterion for transition: when photoionization heating begins to dominate over the PdV cooling and a wind will launch. Below the wind, for an average planet, molecules have not photodissociated so mu should be mean molecular weight instead of mean atomic weight and the erf enforces this. Molecular opacities mean that bolometric heating and cooling dominate the energy budget and create an isotherm in the molecular region below the wind. In the optically thin atomic wind, bolometric heating and cooling are negligible, so the erf also enforces the drop off of bolometric heating and cooling.

If a user needs to call this function for some reason, width_factor is the only argument they should need to change. For numerical or physical reasons it may be necessary for the transition to occur over more than 1 scaleheight. In this case, increase width_factor.

Parameters:

width_factor (float, optional) – Factor to adjust the width in scaleheights of the transition region. Defaults to 0.
_called_in_ramp_bcs (bool, optional) – Set True if called within ramp_base_bcs to avoid recursion issues. Defaults to False.
polish (bool, optional) – If True, computes full postfacto heat_ion & cool_PdV for more accurate transition point. Defaults to False.
bolo_on (bool, optional) – If polish=True, bolo_on=True turns on bolometric heating/cooling if appropriate. Defaults to True.

Returns:

Index of transition if return_idx is True, otherwise status code.

Return type:

int or tuple

wind_simulation.converge_Ncol_sp(expedite=False, integrate_out=False, quiet=False)

Description:: Self-consistently converges Ncol at the sonic point, such that the column density boundary condition for each species, Ncol_sp, matches the number density at the sonic point.

Parameters:

integrating (expedite - bool; default = False. True speeds solutions by not) – past Sonic Point.
Ncol_sp (integrate_out - bool; default = False. Forces outward integration by raising) – incrementally. Final Ncol_sp may be slightly higher than self-consistent.
if (quiet - bool; default = False. Prints warnings about Ncol_sp being only estimates) – integrate_out = False.

wind_simulation.converge_Rmax(final_converge_Ncol=True)

Self-consistently sets Rmax to be the Coriolis length. Does not worry about converging other boundaries. Cannot be expedited as the Coriolis length must be calculated.

Parameters:: final_converge_Ncol (bool, optional) – If True, converges Ncol_sp self-consistently again after Rmax has been set to Rcori. Defaults to True.
Returns:: Status code. 0 for success, other values for failure modes.
Return type:: int

wind_simulation.integrate_out(quiet=False)

Forces solution to integrate out past sonic point to the Rmax. May not preserve self-consistency of per-species column density at sonic point (Ncol_sp) as it may be necessary to raise Ncol_sp in order to integrate outwards. By default Rmax is set to currently-computed R_cori, but can be manually set Rmax to desired value via: .. rubric:: Example

sim.inputs.write_planet_params(sim.windsoln.Rmin,Rmax,*sim.windsoln.bcs_tuple[2:]) sim.run_wind()

Parameters:: quiet (bool, optional) – If True, suppress output messages.
Returns:: 0 if successful, 1 if relaxation error, 4 if Ncol_sp had to be raised
Return type:: int

Computing self consistent boundary conditions

wind_simulation.base_bcs(molec_adjust=None, Kappa_opt=0.004, Kappa_IR=0.01, adiabat=False, base_press=1, user_override_press=False)

Sets base of simulation below XUV tau=1 surface and assumes bolometric heating & cooling of molecules dominates there.

Computes the density and temperature at either R_ubar (microbar pressure radius) or R_IR (vertical tau=1 surface to IR radiation). Between this radius and the nanobar radius where wind is launched, the temperature profile is an isotherm at the skin temperature. If R_IR is below Rp, it takes R_mbar to be the base of the simulation.

Assumes that Rp is the slant path optical tau = 1 surface. Computes skin temperature by balancing bolometric heating and cooling. Computes effective temperature from stellar bolometric luminosity. For highly irradiated planets (most planets for which atmospheric escape is significant), it is an isotherm at T_eff=L_star/(16*pi*sigma_sb*a^2) between R_IR and Rp or an isotherm at T_skin between R_mbar and Rp. If adiabat=True, computes the R_IR assuming an adiabat between T_eff and T_skin.

For high metallicity atmospheres, users should increase the molecular adjustment factor (molec_adjust). Likewise, they should change Kappa_opt and Kappa_IR.

Parameters:

molec_adjust (float, optional) – Dimensionless mean molecular weight. Mu of H2 is 2.3*mH. Defaults to self.windsoln.molec_adjust if None. If molec_adjust == 0, turns off molecular layer.
Kappa_opt (float, optional) – Optical opacity. Defaults to 4e-3.
Kappa_IR (float, optional) – IR opacity. Defaults to 1e-2.
adiabat (bool, optional) – If True, computes IR BCs assuming an adiabat. Defaults to False.
base_press (int or float, optional) – Pressure at base of simulation in microbars. Defaults to 1.
user_override_press (bool, optional) – If True, overrides base pressure. Defaults to False.

Returns:

(R, Rmax, rho, T): R (float): Radius at base of simulation in units of Rp Rmax (float): Set by original boundary conditions, units of Rp rho (float): Density in units of RHO0 from defs.h (1e-15 ~nBAR*MH/K/T0) T (float): Effective temperature in units of T0 (1e4 K)

Return type:

tuple

If R_IR > Rp:: Returns: R_IR, Rmax, rho_IR, T_eff
Else:: Returns: Rp, Rmax, rho_Rp, T_skin

wind_simulation.self_consistent_Ncol(method=1, warning=True)

Computes the self-consistent column density sonic point boundary condition from the neutral number density of a given species.

Parameters:

method (int, optional) – Method for calculation. Defaults to 1.
warning (bool, optional) – If True, prints warnings. Defaults to True.

Returns:

(currents, goals): currents (np.ndarray): Array of current Ncol_sp values. goals (np.ndarray): Array of self-consistent Ncol_sp values.

Return type:

tuple

wind_simulation.turn_off_bolo()

Turns off bolometric heating and cooling AND the mean molecular weight adjustment in the molecular region. These values are currently coupled by the error function defined in converge_mol_atomic_transition().

Returns:: 0 if successful, 1 if failed to turn off bolometric heating/cooling.
Return type:: int

wind_simulation.erf_velocity(return_idx=False, _called_in_ramp_bcs=False, polish=False, width_factor=0, called_in_mol_layer=False)

Defines the drop-off radius of the complementary error function that governs the drop-off of bolometric heating/cooling and the mean molecular weight in the isothermal part of the wind as photoionization heating begins to dominate and the atmosphere becomes atomic and non-isothermal.

Parameters:

return_idx (bool, optional) – If True, returns radial index where drop-off occurs. Defaults to False.
_called_in_ramp_bcs (bool, optional) – Avoids recursion of ramping bases if called in ramp_base_bcs. Defaults to False.
polish (bool, optional) – If True, computes full bolometric heating/cooling and photoionization heating curves for improved accuracy. Defaults to False.
width_factor (float, optional) – Sets width in radius space over which the erfc drops off as a factor of scaleheight at the transition index (Hsc). If 0, uses current width of transition in Hsc. Defaults to 0.

Returns:

Velocity at location of drop-off in units of cm/s. float: Rate of drop-off. tuple (optional): If return_idx=True, also returns (width_factor, drop_idx).

Return type:

float

wind_simulation.raise_Ncol_sp(to_total=0.8, by_factor=0, expedite=False)

Increase per-species column density at the sonic point. This is distinct from converging Ncol_sp. Too low of Ncol_sp may cause numerical instability, so if ramping is stuck, try raising Ncol_sp. Can also diagnose by plotting six_panel_plot—if Ncol_sp drops off sharply after the sonic point for several species, this may be the source of numerical instability.

Parameters:

to_total (float, optional) – If by_factor=0, raises sum(Ncol_sp) until it equals to_total. Defaults to 0.8.
by_factor (float, optional) – If not 0, multiplies Ncol_sp by this factor. Defaults to 0.
expedite (bool, optional) – If True, does not integrate out. Defaults to False.

Returns:

Status code from run_wind.

Return type:

int

Explicitly Ramping Lower Boundary Conditions

wind_simulation.ramp_T_rmin(goal_T, integrate_out=False, _called_in_polish=False)

Ramps normalized temperature at Rmin.

Parameters:

goal_T (float) – Goal temperature at Rmin in units of 1e4 K.
integrate_out (bool, optional) – If True, integrates out after ramping. Defaults to False.

Returns:

0 if successful, 1 if failed to converge.

Return type:

int

wind_simulation.ramp_Rmin(goal_Rmin, integrate_out=False, _called_in_polish=False)

Ramps normalized Rmin, where Rmin is the base of the simulation and should ideally be a radius “below the wind”, i.e., in the molecular region below ~1 microbar.

Parameters:

goal_Rmin (float) – The desired Rmin value in normalized units of Rp.
integrate_out (bool, optional) – If True, integrates out to Coriolis radius after ramping. Defaults to False.

Returns:

0 if successful, 1 if failed to converge.

Return type:

int

wind_simulation.ramp_rho_rmin(goal_rho, integrate_out=False, _called_in_polish=False)

Ramps to normalized mass density at Rmin.

Parameters:

goal_rho (float) – The desired rho value in normalized units of RHO0 (found in src/defs.h, default 1e-15).
integrate_out (bool, optional) – If True, integrates out to Coriolis radius after ramping. Defaults to False.

Returns:

0 if successful, 1 if failed to converge.

Return type:

int

wind_simulation.turn_off_tidal_grav()

Turns off the tidal gravity.

Returns:: None

Ramping spectrum

wind_simulation.flux_norm(goal_flux_in_range, eV_range=[13.6, 100], ramp=False, plot=True, integrate_out=True, converge_bcs=True)

Computes the total flux (across the loaded spectrum range) necessary to achieve the goal_flux_in_range in the eV_range identified. If ramp=True, also ramps the solution to that Ftot.

Parameters:

goal_flux_in_range (float) – Desired flux in ergs/s/cm2 in range given by eV_range.
eV_range (list or array, optional) – Range to normalize to in electron volts (13.6-100 is EUV). Defaults to [13.6, 100].
ramp (bool, optional) – If True, ramps the solution to the computed Ftot. Defaults to False.
plot (bool, optional) – If True, plots ramping progress. Defaults to True.
integrate_out (bool, optional) – If True, integrates out to Coriolis radius after ramping. Defaults to True.
converge_bcs (bool, optional) – If True, polishes boundary conditions after ramping. Defaults to True.

Returns:

0 if ramp successful or already done, error codes (1, 4, etc.) if ramp unsuccessful. float: goal_total_flux if ramp=False (does not ramp).

Return type:

int

wind_simulation.ramp_spectrum(Fnorm=0.0, norm_spec_range=[], goal_spec_range=[], units='eV', normalize=True, kind='full', plot=False)

Ramps stellar spectrum wavelength/energy range to new wavelength/energy range.

Parameters:

Fnorm (float, optional) – Flux in ergs/s/cm2 AT SEMIMAJOR AXIS OF PLANET. If 0.0, flux is normalized to the current value in norm_spec_range. Otherwise, ramps to given Fnorm. Defaults to 0.0.
norm_spec_range (list or array, optional) – Desired range over which to normalize the flux in units of ‘units’. Defaults to [].
goal_spec_range (list or array, optional) – Custom upper and lower limits of spectrum in units of ‘units’. Defaults to [].
units (str, optional) – Units of range values. Options are ‘cm’, ‘nm’, ‘eV’. Defaults to ‘eV’.
normalize (bool, optional) – If True, flux in new range will be normalized to Fnorm in norm_spec_range. Defaults to True.
kind (str, optional) – ‘full’ or ‘mono’ - spectrum frequency type. Defaults to ‘full’.
plot (bool, optional) – If True, plots ramping progress. Defaults to False.

Returns:

0 if ramping successful, 1 if ramping with smaller stepsize unsuccessful, 2 if need to provide Ftot or set norm_flux=True.

Return type:

int

Example

>>> ramp_spectrum(1000,norm_spec_range=[13.6,100],goal_spec_range=[13.6,2000],units='eV',normalize=True)
Will generate a solution with XUV spectrum over 13.6-2000 eV and total integrated flux over 13.6-2000 eV such that the flux in the norm_spec_range (13.6-100 eV) is 1000 ergs/s/cm^2.

wind_simulation.ramp_to_user_spectrum(spectrum_filename, species_list=[], updated_F=0.0, norm_spec_range=[], goal_spec_range=[], units='eV', normalize=True, plot=True, ramp_range=True, hires_savgol_window=None)

Ramping stellar spectrum wavelength/energy range to new wavelength/energy range.

Parameters:

spectrum_filename (str) – Name of the formatted spectrum file saved in McAstro/stars/additional_spectra/.
species_list (list of str, optional) – Species list for spectrum binning. Defaults to [].
updated_F (float, optional) – Summed flux over norm_spec_range in ergs/s/cm² at planet’s semimajor axis. Defaults to 0.0.
norm_spec_range (list or array, optional) – Desired range over which to normalize, in units of ‘units’. Defaults to [].
goal_spec_range (list or array, optional) – Custom upper and lower limits of spectrum, in units of ‘units’. Defaults to [].
units (str, optional) – Units for range values (‘eV’, ‘cm’, ‘nm’). Defaults to ‘eV’.
normalize (bool, optional) – If True, flux in new range will be normalized to Fnorm in norm_spec_range. Defaults to True.
kind (str, optional) – Spectrum frequency type (‘full’ or ‘mono’).
plot (bool, optional) – If True, plot the spectrum. Defaults to False.
ramp_range (bool, optional) – If False, will not ramp spec range or flux. Defaults to True.
hires_savgol_window (int, optional) – Window length for Savitzky–Golay smoothing if user spectrum has >10,000 bins. Must be odd. Defaults to 501 (Higher = more smoothing). Once you have identified an appropriate value, it will be automatically saved to the McAstro/stars/additional_spectra/spectrum_filename for future loads of that spectrum.

Returns:

None

wind_simulation.format_user_spectrum(wl, flux_1au, wl_units='cm', spectrum_name='', comment='', overwrite=False)

Users may import a custom spectrum. This function formats the spectrum to be readable by the code.

Parameters:

wl (list or array) – Wavelengths.
flux_1au (list or array) – Flux per wavelength bin at 1 AU from star in ergs/s/cm2. To scale from Earth, multiply by (distance to star in au / 1 au)^2. Not flux density (to convert from flux density, multiply by delta(wl)).
wl_units (str, optional) – Units for wavelength. Options are ‘cm’, ‘nm’, ‘m’, ‘A’ (angstroms). Defaults to ‘cm’.
spectrum_name (str, optional) – Name of file. Will be saved in solution spectrum. Defaults to ‘’.
comment (str, optional) – Any comments (name of spectrum and other info). No need to add newline character. Defaults to ‘’.
overwrite (bool, optional) – If True, overwrites existing file of same name. Defaults to False.

Returns:

None

Plotting Aliases

wind_simulation.energy_plot(ax=0, alpha=0.8, all_terms=False, CII_line_cool=False, CIII_line_cool=False, OII_line_cool=False, OIII_line_cool=False, legend=True, sub_sonic=True)

Plots energy balance terms used in the energy equation (Broome et al. 2025)

Parameters:

windsoln – The wind solution object containing the simulation data. (“sim.windsoln”)
ax – The axis to plot on (default is 0, which creates a new figure)
alpha – Transparency level for the plot lines (default is 0.8). Useful when overplotting multiple on same axes
all_terms – If True, plot terms not included in Wind-AE, e.g., free-free cooling, conduction (default is False)
CII_line_cool – If True, include CII line cooling terms (default is False)
CIII_line_cool – If True, include CIII line cooling terms (default is False)
OII_line_cool – If True, include OII line cooling terms (default is False)
OIII_line_cool – If True, include OIII line cooling terms (default is False)
legend – If True, display the legend (default is True)
sub_sonic – If True, sets x-axis upper limit at sonic point radius

Returns:

None

wind_simulation.six_panel_plot(Mdot_legend=True, c='k', ls='-', label='', label_dim=[0, 1.3, 2], ion_label=True, first_plotted=True, ax=0)

Plots density (g/cm3), temperature (K), velocity (10 km/s), ionization fraction, column density (g/cm2), and number density (1/cm2), as a function of r (Rp).

Parameters:

object (soln - windsoln)
True (Mdot_legend - Bool; if)
Else (put Mdot in legend of plot.)
prints. (just)
color (c - str; line)
style (ls - str; line)
label (label - str; line)
default=[0 (label_dim - list;)
1.3
[x (2]. Location of label and ncols)
y
ncols].
SixPlot (first_plotted - Bool; True if this the first of many OR the ONLY) – to be plotted on the same axes.
first_plotted=False (ax - matplotlib axis obj; if) – will be be plotted on desired figure with other simulations for comparison
this (provide axis object so) – will be be plotted on desired figure with other simulations for comparison

Returns:

ax - axes object (if first_plotted=True)

Example

ax1 = SixPlot(sim1.windsoln, first_plotted=True) SixPlot(sim2.windsoln, ax=ax1) SixPlot(sim3.windsoln, ax=ax1)

wind_simulation.quick_plot(Mdot_legend=True, c='k', ls='-', label='', label_dim=[0, 1.3, 2], ion_label=True, first_plotted=True, ax=0)

Produces a velocity, density, temperature, and neutral fraction plot from intermediate solutions while ramping. For a more aesthetic four-panel plot use quickplot().

Parameters:

windsoln – windsoln object to be plotted.
ax (np.ndarray) – The ax array on which plots are placed (must have shape = (2,2)).
ax_Ys (matplotlib.axes.Axes, optional) – Axis for plotting ionization fraction. Defaults to None.
label (str, optional) – Label to be added to ax[0][0]. Defaults to None.
alpha (float, optional) – Alpha of lines plotted. Defaults to 0.8.
pvar (list of lists of str, optional) – Variables to plot. Defaults to [[‘v’,’rho’],[‘T’,’Ys_HI’]].
norm (list of lists of floats, optional) – Values to norm (divide) pvars with. Defaults to [[1e5,1e0],[1e0,1e0]].
radius_prefix (str, optional) – SI prefix for x-axis (None defaults to units of Rp). Defaults to None.
sub_sonic (bool, optional) – Only plot sub-sonic region. Defaults to False.
past_rmin (bool, optional) – Only plot past Rmin region. Defaults to False.
sonic_vert (bool, optional) – Add vertical line at sonic point. Defaults to True.

Returns:

None