Tutorial for custom use of exorad_opac
[1]:
import exorad_opac
import astropy.units as u
import astropy.constants as const
import numpy as np
The standard preprocessing routine that is used in exorad_opac_create is called wrap_standard_preprocessing. If you want to do any special things in your preprocessing, you can create a python file or notebook and start from this routine. Here is what the routine does:
[2]:
help(exorad_opac.wrap_standard_preprocessing)
Help on function wrap_standard_preprocessing in module exorad_opac.wrapper:
wrap_standard_preprocessing(path, R, abunds=None, MMW=None, cp=None, stellar_intensity=None, flux_scaling=None, update_datafile=True)
The default workflow.
If you want to do something specefic, you may need to change either of the steps
Steps:
0. Read and parse config file (load_data and parse_config) - done in __init__
1. setup the pressure (p_array and p_center)
2. setup the temperature grid for the premixed quantities (tgrid)
2. generate the Radtrans object - will return and set the Radtrans object
3. Calculate Temperature profile
4. Calculate Chemistry
5. Calculate the stellar intensity
6. Adapt Temperature profile for deep adiabat (if wished for)
7. update the datafile - optional
8. Write output
Note on units:
- Inside exorad_opac we use cgs+bar exclusively
- opac.yaml can be mixed units that are converted during readin (see docs)
Parameters
----------
path: str
Path to exorad simulation folder
R: int or [S0,S1,S2]
Wavelength Resolution
update_datafile: bool
Weather you want to update the datafile with pressure, temperature, cp and Rs
abunds (optional): list of dictionaries (one for each temperature in the tgrid)
Abundances to be used for the opacity grid calculation.
Defaults to the calculated value from equilibrium chemistry
MMW (optional): float
A mean molecular weight to be used.
Defaults to the calculated value from equilibrium chemistry
cp (optional): float
A heat capacity to be used
Defaults to the calculated value from equilibrium chemistry
stellar_intensity (optional): array of same length as Radtrans.freqi
A custom stellar intensity (Flux/pi) scaled to the position of the planet
See Radtrans.get_star_spectrum for a comparable spectrum
Note that the default changes to the phoenix spectrum calculated by
Radtrans.get_star_spectrum
flux_scaling (optional): "BB", None or value
Scale the stellar flux. If value, will use the value for scaling
Have a look into exorad_opac.wrapper to see what it does.
Note: everything is in cgs (+bar for pressure) inside exorad_opac. Only the opac.yaml file can deviate from it.
There are three basic parts:
A parser, that holds some configuration (e.g., from the
opac.yaml)A preprocessing wrapper, that defines the steps for preprocessing (generating opacities and so on)
A output writer class, that writes the output so that
exPERT/MITgcmcan read it
We will go through the steps here:
The parser
The parser objects are meant to hold the configuration of the preprocessing (e.g., Tstar, Rstar, tmin, …). You can either define a custom parser or use the opacyaml parser which loads a configuration from the opacyaml file
There are two parsers: - exorad_opac.OpacYamlParser parses the opacyaml file - exorad_opac.CustomParser sets the most important things
Note: The CustomParser is initialized with all the default values and will only overwrite a few things when its initialized (see help(exorad_opac.CustomParser())). Its inteaded to be used together with the PreprocessingWrapper classes, where non standard parameters should be set in the preprocessing.
This example will use the CustomParser to showcase the preprocessing.
[3]:
Rstar = (1*u.R_sun).cgs.value # cm
Tstar = 5000 # K
gravity = 900 # cm/s2
semimajoraxis = (0.01 * u.au).cgs.value # cm
R = "S1" # wavelength resolution
parser = exorad_opac.CustomParser(R=R, Tstar=Tstar, Rstar=Rstar, semimajoraxis=semimajoraxis, gravity=gravity)
The preprocessing wrapper
There is currently only one preprocessing wrapper which is the PRTPreprocessingWrapper that uses petitRADTRANS to generate k-tables.
[4]:
pwrap = exorad_opac.PRTPreprocessingWrapper(parser)
We now follow through the steps of preprocessing.
First step is to setup the temperature grid to be used for the k-tables. You can also use the default values for tgrid, or if your parser is the OpacYamlParser you can use the tmin, tmax, tresolution ion the opac.yaml to do the same.
[5]:
tgrid = np.logspace(np.log10(100), np.log10(10000), 100)
pwrap.setup_tgrid(tgrid=tgrid)
# pwrap.setup_tgrid() # Use default values, or if provided values from opac.yaml
Next is the pressure grid in the GCM and the ktables.
[6]:
# pwrap.setup_pressure(p0=700, np_log=41, zero_eps=2e-5, dp_lin=100) # setup using specific values
# pwrap.setup_pressure() # Use default values, or if provided values from opac.yaml
pwrap.setup_pressure(p_array=np.logspace(-5, 3, 45)) # setup using specific pressure grid
We can continue and set up the custom petitRADTRANS object
[7]:
pRT = pwrap.generate_pRT_obj(line_species=['H2O_exomol'], rayleigh_species=[], continuum_opacities=[]) # supply any parameter for petitRADTRANS in here!
we are using the showman grid
/Users/schneider/anaconda3/envs/exorad/lib/python3.10/site-packages/petitRADTRANS/radtrans.py:100: FutureWarning: pRT_input_data_path was set by an environment variable. In a future update, the path to the petitRADTRANS input_data will be set within a .ini file that will be automatically generated into the user home directory (OS agnostic), inside a .petitradtrans directory
warnings.warn(f"pRT_input_data_path was set by an environment variable. In a future update, the path to "
Read line opacities of H2O_exomol_R_S1...
Done.
For the initial temperature profile we have quite some options. Head over to exorad_opac/ic.py to see which profiles are implemented!
[8]:
# Right now we just want to use the standard temperature profile that depends on the orbital parameters
pwrap.calc_init_temperature()
# check help(pwrap.calc_init_temperature()) to see your options
Using T_int from Thorngren2019. Update data.exoprt with T_int = 566.8137912496366, if you need consistancy.
[8]:
array([2601.36051658, 2601.36021872, 2601.35975608, 2601.35903472,
2601.3579046 , 2601.356124 , 2601.35329966, 2601.34878448,
2601.34150086, 2601.32963127, 2601.31006954, 2601.27743664,
2601.22229705, 2601.12789741, 2600.96416056, 2600.67657134,
2600.16557362, 2599.24849997, 2597.59005549, 2594.57894981,
2589.1225532 , 2579.35439074, 2562.39286139, 2534.8001666 ,
2495.66650133, 2455.96121931, 2454.60335859, 2556.47943245,
2781.20587561, 3009.62287888, 3122.48101103, 3216.92286116,
3271.2443721 , 3285.03103877, 3297.82717373, 3322.58539926,
3371.68509809, 3472.43157394, 3696.87050679, 4014.30881616,
4339.7374049 , 4670.28978811, 5002.88890446, 5334.35772722])
Now we calculate the chemistry
[9]:
# Use the default equilibrium chemistry package:
pwrap.calc_chemistry()
# ... or specific values
abunds = []
for ti,t in enumerate(pwrap.tgrid):
abunds.append({'H2O':np.ones_like(pwrap.p_array)})
pwrap.calc_chemistry(MMW=2.3, cp=5000, abunds=abunds)
Correct the temperature profile for the deep adiabat. This needs to be done once we have set a MMW and cp in calc_chemistry. Because the adiabat depends on these values.
[10]:
pwrap.correct_init_temperature_for_deep_adiabat(theta_deep=1400, p_max_deep=10, p_min_deep=10)
# pwrap.correct_init_temperature_for_deep_adiabat() # Use default values
[10]:
array([2601.36051658, 2601.36021872, 2601.35975608, 2601.35903472,
2601.3579046 , 2601.356124 , 2601.35329966, 2601.34878448,
2601.34150086, 2601.32963127, 2601.31006954, 2601.27743664,
2601.22229705, 2601.12789741, 2600.96416056, 2600.67657134,
2600.16557362, 2599.24849997, 2597.59005549, 2594.57894981,
2589.1225532 , 2579.35439074, 2562.39286139, 2534.8001666 ,
2495.66650133, 2455.96121931, 2454.60335859, 2556.47943245,
2781.20587561, 3009.62287888, 3122.48101103, 3216.92286116,
3271.2443721 , 3285.03103877, 3297.82717373, 3322.58539926,
3371.68509809, 3472.43157394, 3696.87050679, 4014.30881616,
4339.7374049 , 4670.28978811, 5002.88890446, 5334.35772722])
We can now update the pressure, temperature, atm_Rd and atm_cs in the datafile to be consistent
[11]:
# pwrap.update_datafile(datafile='path/to/data')
We finally setup the stellar intensity. We have a few options here. Either provide the stellar intensity here, or use the inbuild phoenix models. You can decide if you want the stellar intensity to be scaled to a black body or a specific total intensity (Flux/pi). Check help(pwrap.setup_stellar_intensity) for more details on your options.
[12]:
stellar_intensity = np.ones_like(pRT.freq) # use a random stellar intensity
flux_scaling = 3457345
pwrap.setup_stellar_intensity(stellar_intensity=stellar_intensity, flux_scaling=flux_scaling)
# or use the default phoenix stellar model for the given orbital and stellar parameters
# pwrap.setup_stellar_intensity()
The stellar flux is scaled by 3.7359526388671136e-09 to match the input.
Finally calculate the opacity grid
[13]:
pwrap.calc_opa_grid()
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
Region of zero opacity detected, setting the photon destruction probability in this spectral range to 1.
The output
And write out the stuff to your experiment directory
[14]:
output_writer = exorad_opac.pRTOuputWriter(pRT)
# output_writer.write_output("path/to/experiment")