diff --git a/demos/HST/S5_wfc3_template.ecf b/demos/HST/S5_wfc3_template.ecf index 28bea82b..e5d9e875 100644 --- a/demos/HST/S5_wfc3_template.ecf +++ b/demos/HST/S5_wfc3_template.ecf @@ -10,7 +10,7 @@ allapers False # Run S5 on all of the apertures considered in S4? Otherw fit_par ./S5_fit_par_wfc3_template.epf # What fitting epf do you want to use? fit_method [dynesty] #options are: lsq, emcee, dynesty (can list multiple types separated by commas) -run_myfuncs [batman_tr,polynomial,hstramp] # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) +run_myfuncs [batman_tr,polynomial,hstramp] # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) compute_ltt False # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect) # Manual clipping in time diff --git a/demos/JWST/S5_template.ecf b/demos/JWST/S5_template.ecf index 5669dce1..991598ea 100644 --- a/demos/JWST/S5_template.ecf +++ b/demos/JWST/S5_template.ecf @@ -13,7 +13,7 @@ manual_clip None # A list of lists specifying the start and end integrati fit_par ./S5_fit_par_template.epf # What fitting epf do you want to use? fit_method [dynesty] #options are: lsq, emcee, dynesty (can list multiple types separated by commas) -run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) +run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) compute_ltt False # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect) force_positivity False # Optional boolean for sinusoid_pc and poet_pc models. Set True to force positive phase variations. diff --git a/docs/media/S5_template.ecf b/docs/media/S5_template.ecf index 23e991a4..728f4c2a 100644 --- a/docs/media/S5_template.ecf +++ b/docs/media/S5_template.ecf @@ -13,7 +13,7 @@ manual_clip None # A list of lists specifying the start and end integrati fit_par ./S5_fit_par_template.epf # What fitting epf do you want to use? fit_method [dynesty] #options are: lsq, emcee, dynesty (can list multiple types separated by commas) -run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) +run_myfuncs [batman_tr,polynomial] # options are: batman_tr, batman_ecl, catwoman_tr, sinusoid_pc, quasilambert_pc, poet_tr, poet_ecl, poet_pc, fleck_tr, expramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, and GP (can list multiple types separated by commas) compute_ltt False # options are: True (correct model for the light travel time effect), or False (ignore the light travel time effect) force_positivity False # Optional boolean for sinusoid_pc and poet_pc models. Set True to force positive phase variations. diff --git a/docs/source/ecf.rst b/docs/source/ecf.rst index 472fc71b..b87dc1e0 100644 --- a/docs/source/ecf.rst +++ b/docs/source/ecf.rst @@ -887,7 +887,7 @@ run_myfuncs ''''''''''' Determines the astrophysical and systematics models used in the Stage 5 fitting. For standard numpy functions, this can be one or more (separated by commas) of the following: -[batman_tr, batman_ecl, fleck_tr, poet_tr, poet_ecl, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, GP]. +[batman_tr, batman_ecl, catwoman_tr, fleck_tr, poet_tr, poet_ecl, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth, lorentzian, damped_osc, GP]. For theano-based differentiable functions, this can be one or more of the following: [starry, sinusoid_pc, quasilambert_pc, expramp, hstramp, polynomial, step, xpos, ypos, xwidth, ywidth], where starry replaces both the batman_tr and batman_ecl models and offers a more complicated phase variation model than sinusoid_pc that accounts for eclipse mapping signals. @@ -915,6 +915,17 @@ manual_clip ''''''''''' Optional. A list of lists specifying the start and end integration numbers for manual removal. E.g., to remove the first 20 data points specify [[0,20]], and to also remove the last 20 data points specify [[0,20],[-20,None]]. If you want to clip the 10th integration, this would be index 9 since python uses zero-indexing. And the manual_clip start and end values are used to slice a numpy array, so they follow the same convention of *inclusive* start index and *exclusive* end index. In other words, to trim the 10th integrations, you would set manual_clip to [[9,10]]. +Catwoman Convergence Parameters +''''''''''''''''''''''''''''''' +The following two parameters can help in the case of convergence issues when using catwoman_tr. + +catwoman_max_err +^^^^^^^^^^^^^^^^ +The ``max_err`` parameter used by catwoman; defaults to ``1.0``. For more information, see the relevant location of catwoman's `readthedocs page `_ or the relevant location of `catwoman's API `_. + +catwoman_fac +^^^^^^^^^^^^ +The ``fac`` parameter used by catwoman; defaults to ``None``. For more information, see the relevant location of catwoman's `readthedocs page `_ or the relevant location of `catwoman's API `_. Limb Darkening Parameters ''''''''''''''''''''''''' @@ -1091,6 +1102,8 @@ Available fitting parameters are: - Transit and Eclipse Depth Parameters - ``rp`` or ``rprs`` - planet-to-star radius ratio, for the transit models. - ``fp`` or ``fpfs`` - planet-to-star flux ratio, for the eclipse models. + - ``rp2`` or ``rprs2`` - an additional planet-to-star radius ratio for use with the catwoman transit model to model transit limb-asymmetries. + - ``phi`` - the angle (in degrees) of the line separating the semi-circles defined by ``rp`` and ``rp2`` in the catwoman transit model. If ``phi`` is set to 90 degrees (the parameter's default value), the ``rp`` is the trailing hemisphere and ``rp2`` is the leading hemisphere. If ``phi`` is set to 0, then ``rp`` is the northern hemisphere and ``rp2`` is the southern hemisphere. When fitting for multiple planets, add ``_pl#`` after the parameter (e.g., ``rprs``, ``rprs_pl1``, ``rprs_pl2``, etc.). This also applies to the planetaty orbital parameters below. Also be sure to set the ``num_planets`` parameter in your ECF (not EPF) to specify the number of planets being modelled simultaneously. - Orbital Parameters - ``per`` - orbital period (in days) diff --git a/environment.yml b/environment.yml index 83db21a0..e424f8ea 100644 --- a/environment.yml +++ b/environment.yml @@ -36,6 +36,7 @@ dependencies: - tqdm - pip: - astraeus@git+https://github.com/kevin218/Astraeus@main + - catwoman@git+https://github.com/taylorbell57/catwoman@master - celerite2 # Needed for GP - crds<12 # Upper limit needed to avoid bugs with crds.get_context_name - exotic-ld==3.2.0 # Lower limit needed for updated JWST sensitivity files, upper limit needed for breaking changes diff --git a/environment_pymc3.yml b/environment_pymc3.yml index 64a7ecb4..e5f5a8ac 100644 --- a/environment_pymc3.yml +++ b/environment_pymc3.yml @@ -39,6 +39,7 @@ dependencies: - pip: - arviz==0.12.1 - astraeus@git+https://github.com/kevin218/Astraeus@main + - catwoman@git+https://github.com/taylorbell57/catwoman@master - celerite2 # Needed for GP - crds<12 # Upper limit needed to avoid bugs with crds.get_context_name - exoplanet diff --git a/setup.cfg b/setup.cfg index 08fc0fc2..45250d87 100644 --- a/setup.cfg +++ b/setup.cfg @@ -29,6 +29,7 @@ install_requires = astroquery batman-package bokeh<3.0 + catwoman@git+https://github.com/taylorbell57/catwoman@master ccdproc celerite2 # Needed for GP corner diff --git a/src/eureka/S5_lightcurve_fitting/models/AstroModel.py b/src/eureka/S5_lightcurve_fitting/models/AstroModel.py index a06e5967..71170d59 100644 --- a/src/eureka/S5_lightcurve_fitting/models/AstroModel.py +++ b/src/eureka/S5_lightcurve_fitting/models/AstroModel.py @@ -68,6 +68,9 @@ def __init__(self, model, pid=0, channel=0, eval=True): self.t0 = None self.rprs = None self.rp = None + self.rprs2 = None + self.rp2 = None + self.phi = 90. self.inc = None self.ars = None self.a = None @@ -158,6 +161,23 @@ def __init__(self, model, pid=0, channel=0, eval=True): if eval: value = value.value setattr(self, 'rp', value) + # Allow for rp2 or rprs2 + if (self.rprs2 is None) and ('rp2' in model.parameters.dict.keys()): + item0 = 'rp2' + self.pid_id + if model.parameters.dict[item0][1] == 'free': + item0 += self.channel_id + value = getattr(parameterObject, item0) + if eval: + value = value.value + setattr(self, 'rprs2', value) + if (self.rp2 is None) and ('rprs2' in model.parameters.dict.keys()): + item0 = 'rprs2' + self.pid_id + if model.parameters.dict[item0][1] == 'free': + item0 += self.channel_id + value = getattr(parameterObject, item0) + if eval: + value = value.value + setattr(self, 'rp2', value) # Allow for a or ars if (self.ars is None) and ('a' in model.parameters.dict.keys()): item0 = 'a' + self.pid_id diff --git a/src/eureka/S5_lightcurve_fitting/models/CatwomanModel.py b/src/eureka/S5_lightcurve_fitting/models/CatwomanModel.py new file mode 100644 index 00000000..9d9f7650 --- /dev/null +++ b/src/eureka/S5_lightcurve_fitting/models/CatwomanModel.py @@ -0,0 +1,34 @@ +from functools import partial +try: + import catwoman +except ImportError: + print("Could not import catwoman. Functionality may be limited.") + +from .BatmanModels import BatmanTransitModel + + +class CatwomanTransitModel(BatmanTransitModel): + """Transit Model""" + def __init__(self, **kwargs): + """Initialize the transit model + + Parameters + ---------- + **kwargs : dict + Additional parameters to pass to + eureka.S5_lightcurve_fitting.models.Model.__init__(). + Can pass in the parameters, longparamlist, nchan, and + paramtitles arguments here. + """ + # Inherit from BatmanTransitModel class + super().__init__(**kwargs) + self.name = 'catwoman transit' + # Define transit model to be used + self.transit_model = partial(catwoman.TransitModel, + max_err=kwargs['max_err'], + fac=kwargs['fac']) + + if ('rp2' not in self.longparamlist[0] + and 'rprs2' not in self.longparamlist[0]): + raise AssertionError('You must include an rp2 parameter in your ' + 'EPF when using catwoman.') diff --git a/src/eureka/S5_lightcurve_fitting/models/__init__.py b/src/eureka/S5_lightcurve_fitting/models/__init__.py index 22b1519d..40701b26 100644 --- a/src/eureka/S5_lightcurve_fitting/models/__init__.py +++ b/src/eureka/S5_lightcurve_fitting/models/__init__.py @@ -1,6 +1,7 @@ from .Model import Model, CompositeModel from .AstroModel import AstroModel, PlanetParams from .BatmanModels import BatmanTransitModel, BatmanEclipseModel +from .CatwomanModel import CatwomanTransitModel from .FleckModel import FleckTransitModel from .CentroidModel import CentroidModel from .DampedOscillator import DampedOscillatorModel diff --git a/src/eureka/S5_lightcurve_fitting/s5_fit.py b/src/eureka/S5_lightcurve_fitting/s5_fit.py index f217531e..c5416099 100644 --- a/src/eureka/S5_lightcurve_fitting/s5_fit.py +++ b/src/eureka/S5_lightcurve_fitting/s5_fit.py @@ -569,6 +569,7 @@ def fit_channel(meta, time, flux, chan, flux_err, eventlabel, params, else: BatmanTransitModel = m.BatmanTransitModel BatmanEclipseModel = m.BatmanEclipseModel + CatwomanTransitModel = m.CatwomanTransitModel PoetTransitModel = m.PoetTransitModel PoetEclipseModel = m.PoetEclipseModel PoetPCModel = m.PoetPCModel @@ -670,6 +671,27 @@ def fit_channel(meta, time, flux, chan, flux_err, eventlabel, params, nints=lc_model.nints, num_planets=meta.num_planets) modellist.append(t_eclipse) + if 'catwoman_tr' in meta.run_myfuncs: + t_transit = CatwomanTransitModel(parameters=params, + fmt='r--', log=log, time=time, + time_units=time_units, + freenames=freenames, + longparamlist=lc_model.longparamlist, + nchannel=chanrng, + nchannel_fitted=nchannel_fitted, + fitted_channels=fitted_channels, + paramtitles=paramtitles, + ld_from_S4=meta.use_generate_ld, + ld_from_file=meta.ld_file, + ld_coeffs=ldcoeffs, + recenter_ld_prior=meta.recenter_ld_prior, # noqa: E501 + compute_ltt=meta.compute_ltt, + multwhite=lc_model.multwhite, + nints=lc_model.nints, + num_planets=meta.num_planets, + fac=meta.catwoman_fac, + max_err=meta.catwoman_max_err) + modellist.append(t_transit) if 'fleck_tr' in meta.run_myfuncs: t_transit = m.FleckTransitModel(parameters=params, fmt='r--', log=log, time=time, diff --git a/src/eureka/S5_lightcurve_fitting/s5_meta.py b/src/eureka/S5_lightcurve_fitting/s5_meta.py index 8f153eac..63b83c74 100644 --- a/src/eureka/S5_lightcurve_fitting/s5_meta.py +++ b/src/eureka/S5_lightcurve_fitting/s5_meta.py @@ -110,6 +110,10 @@ def set_defaults(self): # Set this to False if not relevant self.recenter_ld_prior = getattr(self, 'recenter_ld_prior', False) + # Catwoman convergence-aiding parameters + self.catwoman_fac = getattr(self, 'catwoman_fac', None) + self.catwoman_max_err = getattr(self, 'catwoman_max_err', 1.0) + # General fitter, fitparams CSV file to resume from self.old_fitparams = getattr(self, 'old_fitparams', None) diff --git a/src/eureka/S6_planet_spectra/s6_spectra.py b/src/eureka/S6_planet_spectra/s6_spectra.py index f3933aa7..b503f4bc 100644 --- a/src/eureka/S6_planet_spectra/s6_spectra.py +++ b/src/eureka/S6_planet_spectra/s6_spectra.py @@ -224,7 +224,8 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None): continue # Manipulate fitted values if needed - if meta.y_param_basic in ['rp^2', 'rprs^2']: + if (meta.y_param_basic[:2] == 'rp' and + meta.y_param_basic[-2:] == '^2'): meta = compute_transit_depth(meta) elif meta.y_param_basic in ['1/r1', '1/r3']: meta = compute_timescale(meta) @@ -234,11 +235,12 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None): if meta.y_label is None: # Provide some default formatting - if meta.y_param_basic in ['rp^2', 'rprs^2']: + if (meta.y_param_basic[:2] == 'rp' and + meta.y_param_basic[-2:] == '^2'): # Transit depth suffix = planetSuffix+channelSuffix meta.y_label = '$(R_{\\rm p'+suffix+'}/R_{\\rm *})^2$' - elif meta.y_param_basic in ['rp', 'rprs']: + elif meta.y_param_basic[:2] == 'rp': # Radius ratio suffix = planetSuffix+channelSuffix meta.y_label = '$R_{\\rm p'+suffix+'}/R_{\\rm *}$' @@ -361,8 +363,7 @@ def plot_spectra(eventlabel, ecf_path=None, s5_meta=None, input_meta=None): # Should we also make the scale_height version of the figure? make_fig6301 = (meta.isplots_S6 >= 3 and meta.has_fig6301reqs - and meta.y_param in ['rp', 'rp^2', - 'rprs', 'rprs^2']) + and meta.y_param[:2] == 'rp') if make_fig6301: # Make spectrum plot with scale height on the 2nd y-axis scale_height = compute_scale_height(meta, log) @@ -414,7 +415,8 @@ def parse_s5_saves(meta, log, fit_methods, channel_key='shared'): errs The uncertainties from a sampling algorithm like dynesty or emcee. """ - if meta.y_param_basic in ['rp^2', 'rprs^2']: + if (meta.y_param_basic[:2] == 'rp' and + meta.y_param_basic[-2:] == '^2'): y_param = meta.y_param[:-2] elif meta.y_param_basic in ['1/r1', '1/r3']: y_param = meta.y_param[2:] @@ -1567,7 +1569,8 @@ def compute_scale_height(meta, log): """ if meta.planet_Rad is None: meta.planet_Rad = meta.spectrum_median - if meta.y_param_basic in ['rp^2', 'rprs^2']: + if (meta.y_param_basic[:2] == 'rp' and + meta.y_param_basic[-2:] == '^2'): meta.planet_Rad = np.sqrt(meta.planet_Rad) meta.planet_Rad = np.nanmean(meta.planet_Rad) meta.planet_Rad *= (meta.star_Rad*constants.R_sun / @@ -1829,7 +1832,8 @@ def transit_latex_table(meta, log): rows = int(np.ceil(nvals/meta.ncols)) # Figure out the labels for the columns - if meta.y_param_basic in ['rp^2', 'rprs^2']: + if (meta.y_param_basic[:2] == 'rp' and + meta.y_param_basic[-2:] == '^2'): suffix = getPlanetSuffix(meta)+getChannelSuffix(meta) colhead = '\\colhead{Transit Depth'+suffix+'}' elif meta.y_param_basic in ['rp', 'rprs']: @@ -1851,7 +1855,7 @@ def transit_latex_table(meta, log): out += "CC|" out = out[:-1]+"}\n" # Give the table a caption based on the tabulated data - if meta.y_param_basic in ['rp', 'rp^2', 'rprs', 'rprs^2']: + if meta.y_param_basic[:2] == 'rp': out += "\\tablecaption{\\texttt{Eureka!}'s Transit Spectroscopy " out += "Results \\label{tab:eureka_transit_spectra}}\n" elif meta.y_param_basic in ['fp', 'fpfs']: diff --git a/src/eureka/lib/citations.py b/src/eureka/lib/citations.py index 1986a421..a16fdd10 100644 --- a/src/eureka/lib/citations.py +++ b/src/eureka/lib/citations.py @@ -735,6 +735,34 @@ adsnote = {Provided by the SAO/NASA Astrophysics Data System} }""" ], + "catwoman": [r"""@ARTICLE{catwoman:joss, + author = {{Jones}, Kathryn and {Espinoza}, N{\'e}stor}, + title = "{catwoman: A transit modelling Python package for asymmetric light curves}", + journal = {The Journal of Open Source Software}, + keywords = {C, exoplanets, Python, transit, C++, astronomy, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - Earth and Planetary Astrophysics}, + year = 2020, + month = nov, + volume = {5}, + number = {55}, + eid = {2382}, + pages = {2382}, + doi = {10.21105/joss.02382}, + archivePrefix = {arXiv}, + eprint = {2106.15643}, + primaryClass = {astro-ph.IM}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2020JOSS....5.2382J}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} + }""", + r"""@software{catwoman:ascl, + author = {{Jones}, Kathryn and {Espinoza}, N{\'e}stor}, + title = "{catwoman: Transit modeling Python package for asymmetric light curves}", + howpublished = {Astrophysics Source Code Library, record ascl:2108.007}, + year = 2021, + month = aug, + eid = {ascl:2108.007}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2021ascl.soft08007J}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} + }"""], "fleck": [r"""@software{2022ascl.soft03009M, author = {{Morris}, Brett M.}, title = "{fleck: Fast starspot rotational modulation light curves}", diff --git a/src/eureka/lib/util.py b/src/eureka/lib/util.py index 9b48f404..f9cc577e 100644 --- a/src/eureka/lib/util.py +++ b/src/eureka/lib/util.py @@ -895,6 +895,8 @@ def make_citations(meta, stage=None): # check if batman or GP is being used for transit/eclipse modeling if "batman_tr" in meta.run_myfuncs or "batman_ecl" in meta.run_myfuncs: other_cites.append("batman") + if "catwoman_tr" in meta.run_myfuncs: + other_cites.append("catwoman") if "fleck_tr" in meta.run_myfuncs: other_cites.append("fleck") if "starry" in meta.run_myfuncs: