From 71f3513e2fe125c6fc17ff2f29bdd52fdb8e4374 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:34:31 +0100 Subject: [PATCH 01/28] Add kilonova-heating-rate to optional requirements --- optional_requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/optional_requirements.txt b/optional_requirements.txt index b7e93647..26ab1ebb 100644 --- a/optional_requirements.txt +++ b/optional_requirements.txt @@ -2,4 +2,5 @@ gwemlightcurves sncosmo lalsimulation nestle -sherpa \ No newline at end of file +sherpa +kilonova-heating-rate \ No newline at end of file From 5f99210b94b019c81060586f709db8c7a6b4064a Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:35:14 +0100 Subject: [PATCH 02/28] Updated transient docs --- docs/transients.txt | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/docs/transients.txt b/docs/transients.txt index 32ebd3c4..de4e5085 100644 --- a/docs/transients.txt +++ b/docs/transients.txt @@ -2,7 +2,7 @@ Transients ============ -A transient in :code:`redback` is implemented as a :code:`redback.transient` object. This class implements the required functionality for fitting all other analysis. +A transient in :code:`redback` is implemented as a :code:`redback.transient` object. This class implements the required functionality for fitting and all other analysis. It also provides a homogenous way to plot the data, or do any processing such as converting integrated flux to luminosity. General transient object @@ -20,7 +20,7 @@ These parent classes are inherited by some transient specific classes, which pro - Afterglow: afterglow of a gamma-ray burst - - SGRB/LGRB: short and long gamma-ray bursts. + - SGRB/LGRB: short and long gamma-ray bursts - Kilonova: For kilonovae - Prompt: For prompt gamma-ray bursts - Supernova: For supernovae of different varieties @@ -29,7 +29,7 @@ These parent classes are inherited by some transient specific classes, which pro These classes come with additional functionality and lookup tables which provide metadata useful for further analysis, such as redshift, T90, start time, etc. They also allow other processing such as converting flux to luminosity. -We provide two method for the former, a simple analytical method, and a more involved method using `sherpa `_ +We provide two methods for converting integrated flux to luminosity, a simple analytical method, and a more involved method using `sherpa `_. For each of the transients we have different :code:`data_modes` which determines what data to fit, plot labels, type of likelihood to use etc. We note that the latter two can be changed by users if desired. @@ -53,7 +53,8 @@ Loading and creating a transient object is simple. .. code:: python kne = 'at2017gfo' - kilonova = redback.Kilonova.from_open_access_catalogue(name=kne, data_mode="flux_density", active_bands=np.array(["g", "i"])) + kilonova = redback.Kilonova.from_open_access_catalogue(name=kne, + data_mode="flux_density", active_bands=np.array(["g", "i"])) This loads the data from at2017gfo that was previously downloaded and creates a kilonova object, with the flux_density data mode. Here we have also specified :code:`active_bands=np.array(['g', 'i')`. This sets the rest of the data to be inactive, i.e., not used in the fitting. @@ -66,10 +67,13 @@ We can use this transient object to create plots of the data. kwargs = {} kilonova.plot_data() fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(12, 8)) - kilonova.plot_multiband(figure=fig, axes=axes, filters=["g", "r", "i", "z", "y", "J"], **kwargs) + kilonova.plot_multiband(figure=fig, axes=axes, + filters=["g", "r", "i", "z", "y", "J"], **kwargs) Here the first `plot_data` will plot the data with all bands on one plot. While the second will plot all filters in the list in separate panels. -We have also passed kwargs here (in this case empty) but can be populated with other keyword arguments to pass to matplotlib. We note that if no figure/axes/filters +We have also passed kwargs here (in this case empty) but can be populated with other keyword arguments to pass to matplotlib. +We have also passed in a `fig` and `axes` to set up the plot in the specific way we wanted. +We note that if no figure/axes/filters are passed then redback will use the defaults. A user can create a different plot_data and plot_multiband functions and pass them in when constructing the class. @@ -87,11 +91,13 @@ Other transient objects can be constructed in a similar manner to the kilonova o afterglow = redback.SGRB.from_swift_grb(name=GRB, data_mode='flux') tidal_disruption_event = redback.TDE.from_open_access_catalogue(tde, data_mode='magnitude') prompt = redback.PromptTimeSeries.from_batse_grb_name(name=name) - supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, data_mode='flux_density', use_phase_model=True) + supernova = redback.supernova.Supernova.from_open_access_catalogue(name=sne, + data_mode='flux_density', use_phase_model=True) Which loads the SGRB, TDE, prompt, and Supernova transient objects with the data for the specific transient respectively. Note that in the supernova object, we set :code:`phase_model=True`. -This sets the time to the actual modified julian date time of the observations, to use when one wants to sample the start time of the transient. +This sets the time attribute of the class to the modified julian date time of the observations. +This is specifically for situations when users want to also sample the start time of the transient. Loading private/simulated data ------------------------- @@ -120,7 +126,8 @@ For example, flux_density_err = data['flux_density_err'].values name = '220501' - afterglow = redback.afterglow(name=name, time=time_days, flux=flux_density, flux_density_err=flux_density_err, frequency=frequency) + afterglow = redback.afterglow(name=name, time=time_days, flux=flux_density, + flux_density_err=flux_density_err, frequency=frequency) We can again plot the data and multiband data From 5e47d9b601344911a121850ac32badcb1c053036 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:35:32 +0100 Subject: [PATCH 03/28] Updated acknowledgements --- README.md | 7 +++++-- docs/acknowledgements.txt | 7 +++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 78f04716..a475d461 100644 --- a/README.md +++ b/README.md @@ -24,8 +24,11 @@ Redback is a software package for end-to-end interpretation and parameter estima ### Contributing Redback is currently in alpha with a paper in preparation. -If you are interested in contributing please join the redback [slack](https://join.slack.com/t/slack-23u4492/shared_invite/zt-14y9q1qmo-VRmc8ZxHzB3u~dB3Wi6pjw). All contributors at the alpha stage will be invited to be co-authors of the first paper. +If you are interested in contributing please join the redback +[slack](https://join.slack.com/t/slack-23u4492/shared_invite/zt-14y9q1qmo-VRmc8ZxHzB3u~dB3Wi6pjw) +and email [Nikhil Sarin](mailto:nikhil.sarin@su.se?subject=Contributing%20to%20redback). +All contributors at the alpha stage will be invited to be co-authors of the first paper. -To make changes to redback, we require users to use a pull request system. +To make changes to redback, we require users to use a merge request system. diff --git a/docs/acknowledgements.txt b/docs/acknowledgements.txt index 6e39d32f..e961dc1a 100644 --- a/docs/acknowledgements.txt +++ b/docs/acknowledgements.txt @@ -9,7 +9,9 @@ please acknowledge the :code:`redback` software, and cite the github page. Using :code:`redback` to fit transients ------------------------- -If you use :code:`redback` to fit an electromagnetic transient, please cite the `Bilby paper `_ +If you use :code:`redback` to fit an electromagnetic transient, +please cite the `Bilby paper `_, +and the :code:`redback` paper (when it is out). Alongside this we request that you also cite the paper associated with the sampler you use. Furthermore, several models implemented in :code:`redback` are introduced in previous scientific publications. @@ -24,11 +26,12 @@ For example, citation = kilonova_models.one_component_kilonova_model.citation -Here citation will be a url to the ads page for the model. +Here citation will be a url to the NASA ads page for the paper describing the model. If you use a specific model or build upon it, we request that you also cite the paper for the model. In several cases, the citation will be `redback`, in which case, only a citation to :code:`redback`, :code:`bilby`, and the sampler is necessary. +Although we recommend periodically checking the citation in the latest :code:`redback` release. Using :code:`redback` to simulate transients ------------------------- From 0a22b9990381784af13d82321145b380e3ff5682 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:35:56 +0100 Subject: [PATCH 04/28] Updated code motivation --- docs/code_motivation.txt | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/docs/code_motivation.txt b/docs/code_motivation.txt index 50d8a731..44c4fdcd 100644 --- a/docs/code_motivation.txt +++ b/docs/code_motivation.txt @@ -2,22 +2,34 @@ Code motivation ============ -How redback can be useful to you -------------------------- - The launch of new telescopes/surveys is leading to an explosion of transient observations. -Redback is a software package for end-to-end interpretation and parameter estimation of these transients. +Redback is a software package that enables end-to-end interpretation and parameter estimation of these transients. +Redback is built with an object oriented modular approach. +This ensures that users can use different parts of :code:`redback` for their own use without needing to interact with other parts. + +How :code:`redback` can be useful to you +------------------------- - Download data for supernovae, tidal disruption events, gamma-ray burst afterglows, kilonovae, prompt emission from - different catalogs/telescopes; Swift, BATSE, Open access catalogs. Users can also provide their own data or use simulated data -- Redback processes the data into a homogeneous transient object, plotting lightcurves and doing other processing. -- The user can then fit one of the models implemented in redback. Or fit their own model. Models for several different types of electromagnetic transients are implemented and range from simple analytical models to numerical surrogates. -- All models are implemented as functions and can be used to simulate populations, without needing to provide data. This way redback can be used simply as a tool to simulate realistic populations, no need to actually fit anything. -- `Bilby `_ under the hood. Can easily switch samplers/likelihoods etc. By default the choice is made depending on the data. -- Fitting returns a homogenous result object, with functionality to plot lightcurves and the posterior/evidence etc. + different catalogs/telescopes; Swift, BATSE, Open access catalogs, ZTF, etc. Users can also provide their own data or use simulated data. +- Process transient data into a homogeneous transient object, providing an interface for plotting lightcurves and doing other processing. +- Fit one of the models implemented in :code:`redback`, or fit your own model. + Models for several different types of electromagnetic transients are implemented and range from simple analytical models to numerical surrogates. +- All models are implemented as functions and can be used to simulate populations, without needing to provide data. +This way :code:`redback` can be used simply as a tool to simulate realistic populations, no need to actually fit anything. +- Fitting returns a homogenous result object, with functionality to plot lightcurves/walkers/corner and the posterior/evidence/credible interval etc. +This way :code:`redback` results can feed into hierarchical analysis of populations of transients or be used in reweighting. -Advantages of the interface to Bilby +Advantages of the interface to :code:`bilby` ------------------------- -As we are using :code:`bilby` under the hood for sampling, priors, likelihood, we have a natural interface to gravitational-wave parameter estimation. This enables easy multi-messenger analysis with redback doing the transient part, and bilby doing the GW part. -We demonstrate this in the `examples `_ \ No newline at end of file +We use `bilby `_ under the hood for inference, which has many advantages. + +- Easily change samplers, likelihoods, place constraints on priors, etc. +- Natural interface with gravitational-wave parameter estimation. + Enabling multi-messenger analysis with :code:`redback` used in fitting to the electromagnetic data, + and :code:`bilby` for gravitational-wave parameter estimation. +- A large developer/user base for core functionality. +`:code:`bilby` is adopted by the LIGO-Virgo-Kagra Collaboration +and used in all parameter estimation results by the LVK collaboration and in over 300 publications, +a testament to its ease of use and robustness. \ No newline at end of file From 22b2df46022da56855d8b4d10107b2d8248b3645 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:36:11 +0100 Subject: [PATCH 05/28] Updated contributing guide --- docs/contributing.txt | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/contributing.txt b/docs/contributing.txt index 8576a700..9097a52e 100644 --- a/docs/contributing.txt +++ b/docs/contributing.txt @@ -6,13 +6,14 @@ Redback is currently in alpha with a paper in preparation. If you are interested in contributing please join the :code:`redback` `slack `_. All contributors at the alpha stage will be invited to be co-authors of the first paper. -To make contributions to :code:`redback` at this stage, -we request that you first make an issue on :code:`redback` github `page `_. +To make contributions to :code:`redback`, +we request that you first make an issue on :code:`redback` github `page `_ describing what contribution you want to make. This minimises the risk of double effort from different people contributing the same functionality. -We also request that you use a merge request and implement some form of a unit test/or example that can be used to test your feature. -A lot of the documentation is automatically generated and this requires -that any feature you implement is well documented following the numpy or restructured text convention. +For all contributions to :code:`redback` we require a merge request system. +We also request that you implement some form of a unit test/or example that can be used to test your feature. +The API documentation is automatically generated and this requires +that every function/class you implement is well documented following the numpy or restructured text convention. Before your request is merged, one of the core developers will look through your implementation and suggest any changes if necessary. This is to ensure that the software remains consistent and easily maintainable. From b4d2e0d124ff6b8e63eca18ae1c8d432309b1505 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:36:33 +0100 Subject: [PATCH 06/28] Updated dependency injection docs --- docs/dependency_injections.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/dependency_injections.txt b/docs/dependency_injections.txt index 7b61a0c1..091c67e6 100644 --- a/docs/dependency_injections.txt +++ b/docs/dependency_injections.txt @@ -3,7 +3,7 @@ Making changes to models and plotting using dependency injections ====== Several models and classes in :code:`redback` use dependency injections. -This allows users to easily swap some functionality with their own preferred method while keeping the rest of the infrastructure intact. +This enables users to easily swap some functionality with their own preferred method while keeping the rest of the infrastructure intact. Modifying a model physics/Creating a new model ------------------------- From fc668a7fe85726f098c96975a1cccade5b4f9f14 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:36:51 +0100 Subject: [PATCH 07/28] fitting docs update --- docs/fitting.txt | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/fitting.txt b/docs/fitting.txt index 5a432de4..a4f6c94d 100644 --- a/docs/fitting.txt +++ b/docs/fitting.txt @@ -7,19 +7,19 @@ and creating a prior we now come to the exciting part; Fitting the model to data To fit our model to data we have to specify a sampler and sampler settings. The likelihood is set by default depending on the transient/data but one can use a different one -or write their own as explained in the likelihood `documentation `_ +or write their own as explained in the likelihood `documentation `_. Installing :code:`redback` with minimal requirements will install the default sampler `dynesty `_. Installing optional requirements will also install `nestle -`_. We generally find `dynesty` to be more reliable/robust but nestle is much faster. +We note that `dynesty` has checkpointing, as do many other samplers. Samplers --------------- As we use :code:`bilby` under the hood, we have access to several different samplers. Cross-checking results with different samplers is often a great way to ensure results are robust -so we encourage users to install multiple samplers and fit models with different samplers. +so we encourage users to install multiple samplers and fit with different samplers. Nested samplers @@ -62,3 +62,7 @@ Here - prior: the prior object - data_mode: type of data to fit. - kwargs: Additional keyword arguments to pass to fit_grb, such as the likelihood, or things required by the sampler, label of the result file, directory where results are saved to etc. + +We note that some samplers have multiprocessing, +which you can see how to use `here https://lscsoft.docs.ligo.org/bilby/api/bilby.core.sampler.run_sampler.html>`_. +We will soon implement some `GPU` models and :code:`parallel bilby` functionality for more rapid inference workflows. \ No newline at end of file From 7a2c15e0bd52fbcca3d62e6e039ee8a163282fc1 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:37:01 +0100 Subject: [PATCH 08/28] grammar --- docs/getting_data.txt | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/docs/getting_data.txt b/docs/getting_data.txt index 1c39d933..96bf4c0f 100644 --- a/docs/getting_data.txt +++ b/docs/getting_data.txt @@ -2,7 +2,7 @@ Get data ============ -Redback provides a simple interface to getting data from the open access, Swift, and BATSE catalogs. +Redback provides a simple interface to getting data from the open access, Swift, and BATSE catalogs, in different formats. We will soon add an interface to get data from the Zwicky Transient Facility, and LASAIR. - Swift: Prompt, X-ray afterglow [counts, flux, flux density] @@ -25,8 +25,8 @@ Or downloading the flux_density/magnitude data of the kilonova at2017gfo data = redback.get_data.get_kilonova_data_from_open_transient_catalog_data(transient=kne) -Both these commands return the data, in a pandas dataframe. They also save the raw data and the processed data in a sensible way. -In particular, the kilonova data will be saved to :code:`kilonova/at2017gfo/` and the afterglow will be saved to :code:`afterglow/GRB070809/flux/` +Both these commands return the data, in a pandas dataframe. They also save the raw and processed data in a sensible way. +In particular, the kilonova data will be saved to :code:`kilonova/at2017gfo/` and the afterglow will be saved to :code:`afterglow/GRB070809/flux/`. Please look at the API or the examples to see how we can get other data. @@ -35,18 +35,20 @@ Basic processing and metadata We do some basic processing to the raw data to make the files homogenous and save them in a homogenous format. We also get metadata about the transient, such as redshift, start time, photon index etc from publically available sources. -Users can also provide this data themselves. +Users can also provide this metadata themselves. Private data or simulated data ------------------------- We do not have to use data from the above catalogs for fitting. Redback can be used on private data or on simulated data. -We have written an example demonstrating this `here `_ +You can see this in the documentation `here `_. +We have also written an example demonstrating this `here `_ as well. + Simulating your own data ------------------------- We will soon add an automated interace to simulating data in specific bands/frequencies/flux/luminosity using the in built :code:`redback` models. For now, users can use any :code:`redback` model as functions, pass in parameters or a dictionary of parameters, and create their own simulated data. -See `API `_. \ No newline at end of file +See `Models documentation `_. \ No newline at end of file From 693001099d8ed602ded66ef14c528a1255a17c52 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:37:16 +0100 Subject: [PATCH 09/28] installation guide --- docs/installation.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/installation.txt b/docs/installation.txt index 284dd80b..742093be 100644 --- a/docs/installation.txt +++ b/docs/installation.txt @@ -23,7 +23,7 @@ Installation This will install all requirements for running :code:`redback` for general transient fitting, including bilby and our default sampler `dynesty `_. Other samplers will need to be -installed via pip or the appropriate means. Please look through :code:`redback` `documentation `_ for what samplers are available. +installed via pip or the appropriate means. Please look through :code:`redback` `fitting documentation `_ for what samplers are available. Install redback from source ------------------------- From 32d5909326d35daffa920fa4a7ffbee3a77f0ad3 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:37:33 +0100 Subject: [PATCH 10/28] likelihood more info --- docs/likelihood.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/likelihood.txt b/docs/likelihood.txt index 38e6d110..b060aa68 100644 --- a/docs/likelihood.txt +++ b/docs/likelihood.txt @@ -3,7 +3,11 @@ Likelihood ============ By default the likelihood is determined by the type of transient/data being used. -However, users can choose a different likelihood (assuming they have an edge case which requires a different likelihood). +However, users can choose a different likelihood. We note that there is typically only one `correct` choice of likelihood but +there may be edge cases such as errors in time, or non-detections, or uncertain y errors which requires users to use a different likelihood. +Many different simple to more complicated likelihoods are included in :code:`redback`, +these should cover most of the cases seen in transient data but if not, users can write their own likelihoods. +We encourage users to add such likelihoods to :code:`redback`. Regular likelihoods ------------------------- From 8726840e25b7566fabe05d6171e7396fd708cd05 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:37:51 +0100 Subject: [PATCH 11/28] Add SNcosmo and more info about datamodes --- docs/models.txt | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/models.txt b/docs/models.txt index 5fb97215..45581d7b 100644 --- a/docs/models.txt +++ b/docs/models.txt @@ -38,6 +38,7 @@ Specifically, the models already included are - magnetar + nickel - SLSN - exponential powerlaw + - SNcosmo - Magnetar boosted ejecta models: @@ -65,7 +66,8 @@ Specifically, the models already included are - 1-6 component piecewise power law - exponential_powerlaw -We note that these models can output in a range of format set by the keyword argument `output_format` or using the appropriate bolometric/flux function. +We note that these models can output in `flux_density` or `magnitude` set by the keyword argument +`output_format` or using the appropriate bolometric/flux function. Alongside these models we also include some general models which can many of the above models as a `base_model` @@ -101,8 +103,9 @@ For example: bolometric_luminosity = function(time, f_nickel=0.6, mej=30, vej=1000, kappa=2, kappa_gamma=1e2) +Here we use `all_models_dict` to provide a simple way to access the relevant function. A user could of course just import the function themselves. + Users can also use the prior objects to get a simulation of the light curves predicted by the function for randomly drawn samples from the prior. -Here we have also demonstrated the all_models_dict for easy access to all different models. .. code:: python @@ -121,12 +124,14 @@ Here we have also demonstrated the all_models_dict for easy access to all differ bolometric_luminosity = function(time, **samples.iloc[0]) Remember that the priors are simply a dictionary so users could also just pass a dictionary/dataframe they created themselves as well. + Users could also sample a lot of different draws from the prior at once (in the above we randomly drew a 100 samples) and then loop through them to simulate a population. Remember that we can only place arbitrary constraints on the priors to make a really specific population/simulation. Modifying :code:`redback` models ------------------------- -A lot of the physics in different redback models is set by default. + +A lot of the physics in different :code:`redback` models is set by default. However, several different pieces of physics in various models can be changed by either passing your own function/class (see next section), by switching the default argument with something else already implemented in redback, or changing a keyword argument. @@ -142,7 +147,7 @@ The specific physics that can be changed: - Different photosphere: See photospheres already implemented `here `_. - Different SED: See SED's already implemented `here `_. -We encourage users to add more of these, which is another easy way to contribute to :code:`redback`. +We encourage users to add more of these physics switches, which is another easy way to contribute to :code:`redback`. From 896d9016de4b39483513608ea9a882da6c1fd518 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:38:14 +0100 Subject: [PATCH 12/28] Syntax and grammar --- docs/pe_basics.txt | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/docs/pe_basics.txt b/docs/pe_basics.txt index 1deb4b65..bf9fbd54 100644 --- a/docs/pe_basics.txt +++ b/docs/pe_basics.txt @@ -2,26 +2,30 @@ Basics of Bayesian inference and parameter estimation ============ -To use redback we assume some level of familiarity with Bayesian inference and model fitting. +To use :code:`redback` we assume some level of familiarity with Bayesian inference and model fitting. However, if this is not the case, :code:`bilby` provides a basic demonstration of Bayesian inference and -how it is implemented in :code:`bilby` for a basic problem of fitting a line is available in the :code:`bilby` `documentation `_ +how it is implemented in :code:`bilby`. +An example for a basic problem of fitting a line is available in the :code:`bilby` `documentation `_. -Redback workflow +:code:`redback` workflow ------------------------- -In redback, we make this process homogenous specifically for fitting electromagnetic transients. The redback workflow for fitting is: +In :code:`redback`, we make this process homogenous specifically for fitting electromagnetic transients. The :code:`redback` workflow for fitting is: - Download the data from a public catalog, or provide your own data. Or simulate it. - Load the data into a homogenous transient object, which does the necessary processing and provides simple way to plot data. - The user then specifies a model (either already implemented in redback or their own function). - Write a prior or use the default priors. -- Specify sampler settings as in :code:`bilby` +- Specify a sampler and sampler settings as in :code:`bilby` - Fit model! - The fit returns a homogenous result object, which can be used for further diagnostics, and provides a simple way to plot the fit. -More advanced functionality +More advanced fitting functionality ------------------------- -- The likelihood is by default set by the type of transient/data used, the more advanced users can provide their own or use more complicated likelihoods implemented in redback. -- Modify the physics of a transient model by passing in different class constructors -- Place constraints on priors if necessary. \ No newline at end of file +- The likelihood is by default set by the type of transient/data used, the more advanced users can provide their own or use more complicated likelihoods implemented in :code:`redback`. +- Modify the physics of a transient model by passing in different class constructors. +- Place constraints on priors if necessary. +- Joint analysis. +- Reweighting. +- Hierarchical inference. \ No newline at end of file From 79d35f0ee30218de4a82083152900313623b4dff Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:38:31 +0100 Subject: [PATCH 13/28] grammar --- docs/priors.txt | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/priors.txt b/docs/priors.txt index f204881a..20db9852 100644 --- a/docs/priors.txt +++ b/docs/priors.txt @@ -1,7 +1,7 @@ ============ Priors ============ -Redback uses bilby priors. See `here `_ for general documentation. +Redback uses :code:`bilby` under the hood for priors. See `here `_ for general documentation of priors in :code:`bilby`. Analytical priors ------------------------- @@ -35,7 +35,7 @@ Interpolated or from file ------------------------- Users can also create a prior from a grid of values i.e., an interpolated_prior. -See documentation `here `_ +See documentation `here `_. Constrained priors ------------------------- @@ -49,6 +49,7 @@ For example, .. code:: python + import redback from bilby.core.prior import PriorDict, Uniform, Constraint priors = PriorDict(conversion_function=redback.constraints.slsn_constraint) @@ -57,4 +58,4 @@ For example, Then define our priors on all other parameters in the normal way. -You can implement your own constraints by following the constraint templates and :code:`bilby` `documentation `_ \ No newline at end of file +You can implement your own constraints by following the constraint templates and :code:`bilby` `documentation `_. \ No newline at end of file From 464f3884a0e463a1290cf8220cedf4085cc220f7 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 15:38:51 +0100 Subject: [PATCH 14/28] add format information --- docs/results.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/results.txt b/docs/results.txt index c7ebf121..f742a3af 100644 --- a/docs/results.txt +++ b/docs/results.txt @@ -10,6 +10,7 @@ such as making pp_plots, walker_plots, getting credible intervals etc. is automa See the :code:`bilby` `API `_ for a full list of features. The result file is by default saved to the `transient/transient_name/model/`, but the user can of course change this. +It is by default saved in a `json` format, which can be changed to 'hdf5' for more compression. Plotting lightcurves and corner plots ------------------------- @@ -21,7 +22,7 @@ Plotting a corner plot is as simple as result.plot_corner() The user can pass in different keyword arguments to change the look/format/what parameters are plotted. -See the :code:`bilby` `API `_ +See the :code:`bilby` `API `_. We can also plot the fit From 3703030fa7ddbbb4c22c1ec6721520abfcfce0e9 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 16:59:35 +0100 Subject: [PATCH 15/28] change frequencies to frequency everywhere. --- redback/sed.py | 48 +++++++++---------- redback/transient_models/kilonova_models.py | 4 +- .../magnetar_boosted_ejecta_models.py | 43 ++++++++--------- redback/utils.py | 18 +++---- 4 files changed, 56 insertions(+), 57 deletions(-) diff --git a/redback/sed.py b/redback/sed.py index 171bd2ec..3245454c 100644 --- a/redback/sed.py +++ b/redback/sed.py @@ -2,33 +2,33 @@ from redback.constants import * from redback.utils import nu_to_lambda, lambda_to_nu -def blackbody_to_flux_density(temperature, r_photosphere, dl, frequencies): +def blackbody_to_flux_density(temperature, r_photosphere, dl, frequency): """ A general blackbody_to_flux_density formula :param temperature: effective temperature in kelvin :param r_photosphere: photosphere radius in cm :param dl: luminosity_distance in cm - :param frequencies: frequencies to calculate in Hz - Must be same length as time array or a single number. In source frame + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number. In source frame :return: flux_density """ ## adding units back in to ensure dimensions are correct - frequencies = frequencies * uu.Hz + frequency = frequency * uu.Hz radius = r_photosphere * uu.cm dl = dl * uu.cm temperature = temperature * uu.K planck = cc.h.cgs speed_of_light = cc.c.cgs boltzmann_constant = cc.k_B.cgs - num = 2 * np.pi * planck * frequencies ** 3 * radius ** 2 + num = 2 * np.pi * planck * frequency ** 3 * radius ** 2 denom = dl ** 2 * speed_of_light ** 2 - frac = 1. / (np.expm1((planck * frequencies) / (boltzmann_constant * temperature))) + frac = 1. / (np.expm1((planck * frequency) / (boltzmann_constant * temperature))) flux_density = num / denom * frac return flux_density class CutoffBlackbody(object): def __init__(self, time, temperature, luminosity, r_photosphere, - frequencies, luminosity_distance, cutoff_wavelength, **kwargs): + frequency, luminosity_distance, cutoff_wavelength, **kwargs): """ Blackbody SED with a cutoff @@ -36,7 +36,7 @@ def __init__(self, time, temperature, luminosity, r_photosphere, :param luminosity: luminosity in cgs :param temperature: temperature in kelvin :param r_photosphere: photosphere radius in cm - :param frequencies: frequencies in Hz - must be a single number or same length as time array + :param frequency: frequency in Hz - must be a single number or same length as time array :param luminosity_distance: dl in cm :param kwargs: None """ @@ -44,7 +44,7 @@ def __init__(self, time, temperature, luminosity, r_photosphere, self.luminosity = luminosity self.temperature = temperature self.r_photosphere = r_photosphere - self.frequencies = frequencies + self.frequency = frequency self.luminosity_distance = luminosity_distance self.cutoff_wavelength = cutoff_wavelength self.reference = 'https://ui.adsabs.harvard.edu/abs/2017ApJ...850...55N/abstract' @@ -55,7 +55,7 @@ def __init__(self, time, temperature, luminosity, r_photosphere, def calculate_flux_density(self): # Mostly from Mosfit/SEDs cutoff_wavelength = self.cutoff_wavelength * angstrom_cgs - wavelength = nu_to_lambda(self.frequencies) + wavelength = nu_to_lambda(self.frequency) x_const = planck * speed_of_light * boltzmann_constant flux_const = 4 * np.pi * 2*np.pi * planck * speed_of_light**2 * angstrom_cgs @@ -107,19 +107,19 @@ def calculate_flux_density(self): class Blackbody(object): - def __init__(self, temperature, r_photosphere, frequencies, luminosity_distance, **kwargs): + def __init__(self, temperature, r_photosphere, frequency, luminosity_distance, **kwargs): """ Simple Blackbody SED :param temperature: effective temperature in kelvin :param r_photosphere: photosphere radius in cm - :param frequencies: frequencies to calculate in Hz - Must be same length as time array or a single number. In source frame + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number. In source frame :param luminosity_distance: luminosity_distance in cm :param kwargs: None """ self.temperature = temperature self.r_photosphere = r_photosphere - self.frequencies = frequencies + self.frequency = frequency self.luminosity_distance = luminosity_distance self.reference = 'It is a blackbody - Do you really need a reference for this?' @@ -127,17 +127,17 @@ def __init__(self, temperature, r_photosphere, frequencies, luminosity_distance, def calculate_flux_density(self): flux_density = blackbody_to_flux_density(temperature=self.temperature, r_photosphere=self.r_photosphere, - frequencies=self.frequencies, dl=self.luminosity_distance) + frequency=self.frequency, dl=self.luminosity_distance) return flux_density class Synchrotron(object): - def __init__(self, frequencies, luminosity_distance, + def __init__(self, frequency, luminosity_distance, pp,nu_max, source_radius=1e13, f0=1e-26, **kwargs): """ Synchrotron SED - :param frequencies: frequencies to calculate in Hz - Must be same length as time array or a single number. In source frame + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number. In source frame :param luminosity_distance: luminosity_distance in cm :param pp: synchrotron power law slope :param nu_max: max frequency @@ -145,7 +145,7 @@ def __init__(self, frequencies, luminosity_distance, :param f0: frequency normalization :param kwargs: None """ - self.frequencies = frequencies + self.frequency = frequency self.luminosity_distance = luminosity_distance self.pp = pp self.nu_max = nu_max @@ -158,16 +158,16 @@ def __init__(self, frequencies, luminosity_distance, def calculate_flux_density(self): fmax = self.f0 * self.source_radius**2 * self.nu_max ** 2.5 # for SSA - mask = self.frequencies < self.nu_max - sed = np.zeros(len(self.frequencies)) + mask = self.frequency < self.nu_max + sed = np.zeros(len(self.frequency)) # sed units are erg/s/Angstrom - need to turn them into flux density compatible units units = uu.erg / uu.s / uu.hz / uu.cm**2 sed[mask] = self.f0 * self.source_radius**2 * \ - (self.frequencies/self.nu_max)**2.5 * angstrom_cgs / speed_of_light * self.frequencies **2 - sed[~mask] = fmax * (self.frequencies/self.nu_max)**(-(self.pp - 1.)/2.) \ - * angstrom_cgs / speed_of_light * self.frequencies **2 + (self.frequency/self.nu_max)**2.5 * angstrom_cgs / speed_of_light * self.frequency **2 + sed[~mask] = fmax * (self.frequency/self.nu_max)**(-(self.pp - 1.)/2.) \ + * angstrom_cgs / speed_of_light * self.frequency **2 self.sed = sed sed = sed / (4*np.pi * self.luminosity_distance**2) * lambda_to_nu(1.) @@ -186,7 +186,7 @@ def __init__(self, time, luminosity, frequency, sed, luminosity_distance, line_w :param time: time in source frame :param luminosity: luminosity in cgs - :param frequency: frequencies to calculate in Hz - Must be same length as time array or a single number. In source frame + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number. In source frame :param sed: instantiated SED class object. :param luminosity_distance: luminosity_distance in cm :param line_wavelength: line wavelength in angstrom @@ -198,7 +198,7 @@ def __init__(self, time, luminosity, frequency, sed, luminosity_distance, line_w """ self.time = time self.luminosity = luminosity - self.frequencies = frequency + self.frequency = frequency self.SED = sed self.luminosity_distance = luminosity_distance self.line_wavelength = line_wavelength @@ -214,7 +214,7 @@ def __init__(self, time, luminosity, frequency, sed, luminosity_distance, line_w def calculate_flux_density(self): # Mostly from Mosfit/SEDs - wavelength = nu_to_lambda(self.frequencies) + wavelength = nu_to_lambda(self.frequency) amplitude = self.line_amplitude * np.exp(-0.5*((self.time - self.line_time)/self.line_duration)**2) seds = self.SED.sed * (1 - amplitude) diff --git a/redback/transient_models/kilonova_models.py b/redback/transient_models/kilonova_models.py index c03ad68c..b2c5a62c 100644 --- a/redback/transient_models/kilonova_models.py +++ b/redback/transient_models/kilonova_models.py @@ -267,13 +267,13 @@ def one_component_kilonova_model(time, redshift, mej, vej, kappa, **kwargs): temp_func = interp1d(time_temp, y=temperature) rad_func = interp1d(time_temp, y=r_photosphere) # convert to source frame time and frequency - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) temp = temp_func(time) photosphere = rad_func(time) flux_density = blackbody_to_flux_density(temperature=temp, r_photosphere=photosphere, - dl=dl, frequencies=frequency) + dl=dl, frequency=frequency) if kwargs['output_format'] == 'flux_density': return flux_density.to(uu.mJy).value diff --git a/redback/transient_models/magnetar_boosted_ejecta_models.py b/redback/transient_models/magnetar_boosted_ejecta_models.py index e42327c1..758aaaea 100644 --- a/redback/transient_models/magnetar_boosted_ejecta_models.py +++ b/redback/transient_models/magnetar_boosted_ejecta_models.py @@ -14,7 +14,6 @@ def metzger_magnetar_boosted_kilonova_model(time, redshift, mej, vej, beta, kapp """ :param time: observer frame time in days :param redshift: redshift - :param frequencies: frequencies to calculate - Must be same length as time array or a single number :param mej: ejecta mass in solar masses :param vej: minimum initial velocity :param beta: velocity power law slope (M=v^-beta) @@ -24,10 +23,10 @@ def metzger_magnetar_boosted_kilonova_model(time, redshift, mej, vej, beta, kapp :param nn: braking index :param thermalisation_efficiency: magnetar thermalisation efficiency :param kwargs: neutron_precursor_switch, pair_cascade_switch, ejecta_albedo, magnetar_heating, output_format - frequencies (frequencies to calculate - Must be same length as time array or a single number) + frequency (frequency to calculate - Must be same length as time array or a single number) :return: flux_density or magnitude """ - frequencies = kwargs['frequencies'] + frequency = kwargs['frequency'] time_temp = np.geomspace(1e-4, 1e7, 300) bolometric_luminosity, temperature, r_photosphere = _metzger_magnetar_boosted_kilonova_model(time_temp, mej, vej, beta, kappa_r, l0, tau_sd, nn, @@ -39,13 +38,13 @@ def metzger_magnetar_boosted_kilonova_model(time, redshift, mej, vej, beta, kapp rad_func = interp1d(time_temp, y=r_photosphere) # convert to source frame time and frequency time = time * 86400 - frequencies, time = calc_kcorrected_properties(frequencies=frequencies, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) temp = temp_func(time) photosphere = rad_func(time) flux_density = blackbody_to_flux_density(temperature=temp, r_photosphere=photosphere, - dl=dl, frequencies=frequencies) + dl=dl, frequency=frequency) if kwargs['output_format'] == 'flux_density': return flux_density.to(uu.mJy).value @@ -288,48 +287,48 @@ def _ejecta_dynamics_and_interaction(time, mej, beta, ejecta_radius, kappa, n_is return lorentz_factor, lbol_rest, comoving_temperature, radius, doppler_factor, tau -def _comoving_blackbody_to_flux_density(dl, frequencies, radius, temperature, doppler_factor): +def _comoving_blackbody_to_flux_density(dl, frequency, radius, temperature, doppler_factor): """ :param dl: luminosity distance in cm - :param frequencies: frequencies to calculate in Hz - Must be same length as time array or a single number + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number :param radius: ejecta radius in cm :param temperature: comoving temperature in K :param doppler_factor: doppler_factor :return: flux_density """ ## adding units back in to ensure dimensions are correct - frequencies = frequencies * uu.Hz + frequency = frequency * uu.Hz radius = radius * uu.cm dl = dl * uu.cm temperature = temperature * uu.K planck = cc.h.cgs speed_of_light = cc.c.cgs boltzmann_constant = cc.k_B.cgs - num = 2 * np.pi * planck * frequencies ** 3 * radius ** 2 + num = 2 * np.pi * planck * frequency ** 3 * radius ** 2 denom = dl ** 2 * speed_of_light ** 2 * doppler_factor ** 2 - frac = 1. / (np.exp((planck * frequencies) / (boltzmann_constant * temperature * doppler_factor)) - 1) + frac = 1. / (np.exp((planck * frequency) / (boltzmann_constant * temperature * doppler_factor)) - 1) flux_density = num / denom * frac return flux_density -def _comoving_blackbody_to_luminosity(frequencies, radius, temperature, doppler_factor): +def _comoving_blackbody_to_luminosity(frequency, radius, temperature, doppler_factor): """ - :param frequencies: frequencies to calculate in Hz - Must be same length as time array or a single number + :param frequency: frequency to calculate in Hz - Must be same length as time array or a single number :param radius: ejecta radius in cm :param temperature: comoving temperature in K :param doppler_factor: doppler_factor :return: luminosity """ ## adding units back in to ensure dimensions are correct - frequencies = frequencies * uu.Hz + frequency = frequency * uu.Hz radius = radius * uu.cm temperature = temperature * uu.K planck = cc.h.cgs speed_of_light = cc.c.cgs boltzmann_constant = cc.k_B.cgs - num = 8 * np.pi ** 2 * planck * frequencies ** 4 * radius ** 2 + num = 8 * np.pi ** 2 * planck * frequency ** 4 * radius ** 2 denom = speed_of_light ** 2 * doppler_factor ** 2 - frac = 1. / (np.exp((planck * frequencies) / (boltzmann_constant * temperature * doppler_factor)) - 1) + frac = 1. / (np.exp((planck * frequency) / (boltzmann_constant * temperature * doppler_factor)) - 1) luminosity = num / denom * frac return luminosity @@ -349,10 +348,10 @@ def mergernova(time, redshift, mej, beta, ejecta_radius, kappa, n_ism, l0, tau_s :param nn: braking index :param thermalisation_efficiency: magnetar thermalisation efficiency :param kwargs: output_format - whether to output flux density or AB magnitude - frequencies (frequencies to calculate - Must be same length as time array or a single number) + frequency (frequency to calculate - Must be same length as time array or a single number) :return: flux density or AB magnitude """ - frequencies = kwargs['frequencies'] + frequency = kwargs['frequency'] time_temp = np.geomspace(1e-4, 1e8, 1000, endpoint=True) dl = cosmo.luminosity_distance(redshift).cgs.value _, bolometric_luminosity, comoving_temperature, radius, doppler_factor, _ = _ejecta_dynamics_and_interaction( @@ -366,12 +365,12 @@ def mergernova(time, redshift, mej, beta, ejecta_radius, kappa, n_ism, l0, tau_s d_func = interp1d(time_temp, y=doppler_factor) # convert to source frame time and frequency time = time * 86400 - frequencies, time = calc_kcorrected_properties(frequencies=frequencies, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) temp = temp_func(time) rad = rad_func(time) df = d_func(time) - flux_density = _comoving_blackbody_to_flux_density(dl=dl, frequencies=frequencies, radius=rad, temperature=temp, + flux_density = _comoving_blackbody_to_flux_density(dl=dl, frequency=frequency, radius=rad, temperature=temp, doppler_factor=df) if kwargs['output_format'] == 'flux_density': return flux_density.to(uu.mJy).value @@ -413,7 +412,7 @@ def _trapped_magnetar_lum(time, mej, beta, ejecta_radius, kappa, n_ism, l0, tau_ df = d_func(time) optical_depth = tau_func(time) frequency = kwargs['frequency'] - trapped_ejecta_lum = _comoving_blackbody_to_luminosity(frequencies=frequency, radius=rad, + trapped_ejecta_lum = _comoving_blackbody_to_luminosity(frequency=frequency, radius=rad, temperature=temp, doppler_factor=df) lsd = _magnetar_only(time, l0=l0, tau=tau_sd, nn=nn) lum = np.exp(-optical_depth) * lsd + trapped_ejecta_lum @@ -439,8 +438,8 @@ def _trapped_magnetar_flux(time, redshift, mej, beta, ejecta_radius, kappa, n_is :param kwargs: 'photon_index' used to calculate k correction and convert from luminosity to flux :return: integrated flux """ - frequencies = kwargs['frequency'] - frequencies, time = calc_kcorrected_properties(frequencies=frequencies, redshift=redshift, time=time) + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) lum = _trapped_magnetar_lum(time, mej, beta, ejecta_radius, kappa, n_ism, l0, tau_sd, nn, thermalisation_efficiency, diff --git a/redback/utils.py b/redback/utils.py index 40f586af..898346d8 100644 --- a/redback/utils.py +++ b/redback/utils.py @@ -48,17 +48,17 @@ def nu_to_lambda(frequency): """ return 1.e10 * (speed_of_light_si / frequency) -def calc_kcorrected_properties(frequencies, redshift, time): +def calc_kcorrected_properties(frequency, redshift, time): """ Perform k-correction - :param frequencies: observer frame frequencies + :param frequency: observer frame frequency :param redshift: source redshift :param time: observer frame time - :return: k-corrected frequencies and source frame time + :return: k-corrected frequency and source frame time """ time = time / (1 + redshift) - frequencies = frequencies * (1 + redshift) - return frequencies, time + frequency = frequency * (1 + redshift) + return frequency, time def mjd_to_jd(mjd): @@ -146,7 +146,7 @@ def date_to_mjd(year, month, day): return mjd -def get_filter_frequencies(filter): +def get_filter_frequency(filter): pass @@ -205,16 +205,16 @@ def calc_flux_from_mag(magnitude, reference_flux, magnitude_system='AB'): return 1000 * flux # return in mJy -def bands_to_frequencies(bands): +def bands_to_frequency(bands): """ - Converts a list of bands into an array of frequencies in Hz + Converts a list of bands into an array of frequency in Hz ---------- bands: list[str] List of bands. Returns ---------- - array_like: An array of frequencies associated with the given bands. + array_like: An array of frequency associated with the given bands. """ if bands is None: bands = [] From 352d6761e025b8395b858c20ad054b6fe2118a4e Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 17:46:17 +0100 Subject: [PATCH 16/28] Documentation, energy injection function, and default kwargs --- redback/transient_models/afterglow_models.py | 404 ++++++++++++++++--- 1 file changed, 338 insertions(+), 66 deletions(-) diff --git a/redback/transient_models/afterglow_models.py b/redback/transient_models/afterglow_models.py index 6f744769..3fa3e42a 100644 --- a/redback/transient_models/afterglow_models.py +++ b/redback/transient_models/afterglow_models.py @@ -1,6 +1,6 @@ from astropy.cosmology import Planck18 as cosmo # noqa from inspect import isfunction -from redback.utils import logger, citation_wrapper +from redback.utils import logger, citation_wrapper, calc_ABmag_from_flux_density from redback.constants import day_to_s try: import afterglowpy as afterglow @@ -22,13 +22,41 @@ @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def cocoon(time, redshift, umax, umin, loge0, k, mej, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A cocoon afterglow model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param umax: initial outflow 4 velocity maximum + :param umin: minimum outflow 4 velocity + :param loge0: log10 fidicial energy in velocity distribution E(>u) = E0u^-k in erg + :param k: power law index of energy velocity distribution + :param mej: mass of material at umax in solar masses + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['cocoon'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] e0 = 10 ** loge0 n0 = 10 ** logn0 @@ -36,41 +64,83 @@ def cocoon(time, redshift, umax, umin, loge0, k, mej, logn0, p, logepse, logepsb epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'uMax': umax, 'Er': e0, 'uMin': umin, 'k': k, 'MFast_solar': mej, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def kn_afterglow(time, redshift, umax, umin, loge0, k, mej, logn0, p, logepse, logepsb, ksin, g0, **kwargs): - time = time * day_to_s - dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] - jettype = jettype_dict['cocoon'] - spectype = kwargs['spectype'] - frequency = kwargs['frequency'] - e0 = 10 ** loge0 - n0 = 10 ** logn0 - epse = 10 ** logepse - epsb = 10 ** logepsb - Z = {'jetType': jettype, 'specType': spectype, 'uMax': umax, 'Er': e0, - 'uMin': umin, 'k': k, 'MFast_solar': mej, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, - 'spread': spread, 'latRes': latres, 'tRes': tres} - flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + """ + A kilonova afterglow model from afterglowpy, similar to cocoon but with constraints. + + :param time: time in days in source frame + :param redshift: source redshift + :param umax: initial outflow 4 velocity maximum + :param umin: minimum outflow 4 velocity + :param loge0: log10 fidicial energy in velocity distribution E(>u) = E0u^-k in erg + :param k: power law index of energy velocity distribution + :param mej: mass of material at umax in solar masses + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ + + output = cocoon(time=time, redshift=redshift, umax=umax, umin=umin, loge0=loge0, + k=k, mej=mej, logn0=logn0,p=p,logepse=logepse,logepsb=logepsb, + ksin=ksin, g0=g0, **kwargs) + return output @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def cone_afterglow(time, redshift, thv, loge0, thw, thc, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A cone afterglow model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thw: wing truncation angle of jet thw = thw*thc + :param thc: half width of jet core in radians + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['cone'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] thw = thw * thc e0 = 10 ** loge0 @@ -79,20 +149,50 @@ def cone_afterglow(time, redshift, thv, loge0, thw, thc, logn0, p, logepse, loge epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres, 'thetaWing': thw} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def gaussiancore(time, redshift, thv, loge0, thc, thw, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A gaussiancore model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thw: wing truncation angle of jet thw = thw*thc + :param thc: half width of jet core in radians + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['gaussian_w_core'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] thw = thw * thc @@ -102,20 +202,50 @@ def gaussiancore(time, redshift, thv, loge0, thc, thw, logn0, p, logepse, logeps epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres, 'thetaWing': thw} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def gaussian(time, redshift, thv, loge0, thw, thc, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A gaussian structured jet model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thw: wing truncation angle of jet thw = thw*thc + :param thc: half width of jet core in radians + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['gaussian'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] thw = thw * thc e0 = 10 ** loge0 @@ -124,20 +254,51 @@ def gaussian(time, redshift, thv, loge0, thw, thc, logn0, p, logepse, logepsb, k epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres, 'thetaWing': thw} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def smoothpowerlaw(time, redshift, thv, loge0, thw, thc, beta, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A smoothpowerlaw structured jet model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thw: wing truncation angle of jet thw = thw*thc + :param thc: half width of jet core in radians + :param beta: index for power-law structure, theta^-b + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['smooth_power_law'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] thw = thw * thc e0 = 10 ** loge0 @@ -146,57 +307,136 @@ def smoothpowerlaw(time, redshift, thv, loge0, thw, thc, beta, logn0, p, logepse epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres, 'thetaWing': thw, 'b': beta} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def powerlawcore(time, redshift, thv, loge0, thw, thc, beta, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A power law with core structured jet model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thw: wing truncation angle of jet thw = thw*thc + :param thc: half width of jet core in radians + :param beta: index for power-law structure, theta^-b + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['powerlaw_w_core'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] thw = thw * thc e0 = 10 ** loge0 n0 = 10 ** logn0 epse = 10 ** logepse epsb = 10 ** logepsb + Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres, 'thetaWing': thw, 'b': beta} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') def tophat(time, redshift, thv, loge0, thc, logn0, p, logepse, logepsb, ksin, g0, **kwargs): + """ + A tophat jet model from afterglowpy + + :param time: time in days in source frame + :param redshift: source redshift + :param thv: viewing angle in radians + :param loge0: log10 on axis isotropic equivalent energy + :param thc: half width of jet core/jet opening angle in radians + :param beta: index for power-law structure, theta^-b + :param logn0: log10 number density of ISM in cm^-3 + :param p: electron distribution power law index. Must be greater than 2. + :param logepse: log10 fraction of thermal energy in electrons + :param logepsb: log10 fraction of thermal energy in magnetic field + :param ksin: fraction of electrons that get accelerated + :param g0: initial lorentz factor + :param kwargs: spread: whether jet can spread, defaults to False + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + l0, ts, q: energy injection parameters, defaults to 0 + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + :return: flux density or AB mag. + """ time = time * day_to_s dl = cosmo.luminosity_distance(redshift).cgs.value - spread = kwargs['spread'] - latres = kwargs['latres'] - tres = kwargs['tres'] + spread = kwargs.get('spread', False) + latres = kwargs.get('latres', 2) + tres = kwargs.get('tres', 100) + spectype = kwargs.get('spectype', 0) + l0 = kwargs.get('L0', 0) + q = kwargs.get('q', 0) + ts = kwargs.get('ts', 0) jettype = jettype_dict['tophat'] - spectype = kwargs['spectype'] frequency = kwargs['frequency'] - e0 = 10 ** loge0 n0 = 10 ** logn0 epse = 10 ** logepse epsb = 10 ** logepsb Z = {'jetType': jettype, 'specType': spectype, 'thetaObs': thv, 'E0': e0, 'thetaCore': thc, 'n0': n0, 'p': p, 'epsilon_e': epse, 'epsilon_B': epsb, - 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': 0, 'q': 0, 'ts': 0, 'g0': g0, + 'xi_N': ksin, 'd_L': dl, 'z': redshift, 'L0': l0, 'q': q, 'ts': ts, 'g0': g0, 'spread': spread, 'latRes': latres, 'tRes': tres} flux_density = afterglow.fluxDensity(time, frequency, **Z) - return flux_density + if kwargs['output_format'] == 'flux_density': + return flux_density + elif kwargs['output_format'] == 'magnitude': + return calc_ABmag_from_flux_density(flux_density).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') -def afterglow_models_no_jet_spread(time, **kwargs): +def afterglow_models_with_energy_injection(time, **kwargs): + """ + A base class for afterglowpy models with energy injection. + :param time: time in days in source frame + :param kwargs: all kwargs used by the specific jet model. + l0: Fiducial luminosity for energy injection + q: temporal powerlaw index for energy injection + ts: fiducial timescale for energy injection in days + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + base_model: A string to indicate the type of jet model to use. + :return: flux density or AB mag. + """ from redback.model_library import modules_dict # import model library in function to avoid circular dependency base_model = kwargs['base_model'] if isfunction(base_model): @@ -208,7 +448,39 @@ def afterglow_models_no_jet_spread(time, **kwargs): function = modules_dict['afterglow_models'][base_model] else: raise ValueError("Not a valid base model.") - kwargs['spread'] = False + kwargs['ts'] = kwargs['ts'] * day_to_s + kwargs['spread'] = True + output = function(time, **kwargs) + return output + +@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract') +def afterglow_models_with_jet_spread(time, **kwargs): + """ + A base class for afterglow models with jet spreading. Note, with these models you cannot sample in g0. + + :param time: time in days in source frame + :param kwargs: all kwargs used by the specific jet model. + latres: latitudinal resolution for structured jets, defaults to 2 + tres: time resolution of shock evolution, defaults to 100 + spectype: whether to have inverse compton, defaults to 0, i.e., no inverse compton. + change to 1 for including inverse compton emission. + output_format: Whether to output flux density or AB mag + base_model: A string to indicate the type of jet model to use. + :return: flux density or AB mag. + """ + from redback.model_library import modules_dict # import model library in function to avoid circular dependency + base_model = kwargs['base_model'] + if isfunction(base_model): + function = base_model + elif base_model not in jet_spreading_models: + logger.warning('{} is not implemented as a base model'.format(base_model)) + raise ValueError('Please choose a different base model') + elif isinstance(base_model, str): + function = modules_dict['afterglow_models'][base_model] + else: + raise ValueError("Not a valid base model.") + kwargs['spread'] = True kwargs.pop('g0') - flux_density = function(time, **kwargs) - return flux_density \ No newline at end of file + output = function(time, **kwargs) + return output + From e2acdbd242fd585410f1615be7bd68840d8547a6 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 17:50:47 +0100 Subject: [PATCH 17/28] change electron fraction to pair cascade fraction in pair cascades --- .../magnetar_boosted_ejecta_models.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/redback/transient_models/magnetar_boosted_ejecta_models.py b/redback/transient_models/magnetar_boosted_ejecta_models.py index 758aaaea..7e08fb0c 100644 --- a/redback/transient_models/magnetar_boosted_ejecta_models.py +++ b/redback/transient_models/magnetar_boosted_ejecta_models.py @@ -23,7 +23,8 @@ def metzger_magnetar_boosted_kilonova_model(time, redshift, mej, vej, beta, kapp :param nn: braking index :param thermalisation_efficiency: magnetar thermalisation efficiency :param kwargs: neutron_precursor_switch, pair_cascade_switch, ejecta_albedo, magnetar_heating, output_format - frequency (frequency to calculate - Must be same length as time array or a single number) + frequency (frequency to calculate - Must be same length as time array or a single number), + pair_cascade_fraction: fraction of magnetar spin down energy that turns into pair cascades :return: flux_density or magnitude """ frequency = kwargs['frequency'] @@ -63,15 +64,13 @@ def _metzger_magnetar_boosted_kilonova_model(time, mej, vej, beta, kappa_r, l0, :param tau_sd: magnetar spin down damping timescale :param nn: braking index :param thermalisation_efficiency: magnetar thermalisation efficiency - :param kwargs: neutron_precursor_switch, pair_cascade_switch, ejecta_albedo, magnetar_heating + :param kwargs: neutron_precursor_switch, pair_cascade_switch, ejecta_albedo, magnetar_heating, pair_cascade_fraction :return: bolometric_luminosity, temperature, photosphere_radius """ pair_cascade_switch = kwargs.get('pair_cascade_switch', True) neutron_precursor_switch = kwargs.get('neutron_precursor_switch', True) magnetar_heating = kwargs.get('magnetar_heating', 'first_layer') - if kwargs['pair_cascade_switch'] == True: - ejecta_albedo = kwargs.get('ejecta_albedo', 0.5) time = time tdays = time/86400 @@ -195,7 +194,9 @@ def _metzger_magnetar_boosted_kilonova_model(time, mej, vej, beta, kappa_r, l0, bolometric_luminosity = np.sum(lum_rad, axis=0) if pair_cascade_switch == True: - tlife_t = (0.6/(1 - ejecta_albedo))*(electron_fraction/0.1)**0.5 * (bolometric_luminosity/1.0e45)**0.5 \ + ejecta_albedo = kwargs.get('ejecta_albedo', 0.5) + pair_cascade_fraction = kwargs.get('pair_cascade_fraction', 0.01) + tlife_t = (0.6/(1 - ejecta_albedo))*(pair_cascade_fraction/0.1)**0.5 * (lsd/1.0e45)**0.5 \ * (v0/(0.3*speed_of_light))**(0.5) * (time/86400)**(-0.5) bolometric_luminosity = bolometric_luminosity / (1.0 + tlife_t) From 7f13d89499f6664c115c8946c811f6a47c416ac1 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 18:03:05 +0100 Subject: [PATCH 18/28] frequencies to frequency --- redback/transient_models/supernova_models.py | 24 +++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index 5a1cfdb0..ca791546 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -60,7 +60,7 @@ def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, """ frequency = kwargs['frequency'] # time = time * 86400 - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = exponential_powerlaw_bolometric(time=time, lbol_0=lbol_0, @@ -68,7 +68,7 @@ def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, interaction_process=interaction_process, **kwargs) photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl) + frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -133,13 +133,13 @@ def arnett(time, redshift, f_nickel, mej, interaction_process=ip.Diffusion, """ frequency = kwargs['frequency'] # time = time * 86400 - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = arnett_bolometric(time=time, f_nickel=f_nickel, mej=mej, interaction_process=interaction_process, **kwargs) photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl) + frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -203,7 +203,7 @@ def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interactio :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = basic_magnetar_powered_bolometric(time=time, p0=p0,bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, @@ -211,7 +211,7 @@ def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interactio photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl) + frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -263,13 +263,13 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu """ frequency = kwargs['frequency'] cutoff_wavelength = kwargs.get('cutoff_wavelength', 3000) - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = slsn_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, interaction_process=interaction_process) photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) + frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) flux_density = sed_1.flux_density @@ -300,7 +300,7 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol_mag = _basic_magnetar(time=time*86400, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) @@ -314,7 +314,7 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl) + frequency=frequency, luminosity_distance=dl) flux_density = sed_1.flux_density @@ -340,7 +340,9 @@ def type_1c(): pass @citation_wrapper('redback') -def homologous_expansion_supernova_model_bolometric(): +def homologous_expansion_supernova_model_bolometric(time, **kwargs): + v_ejecta = np.sqrt(10.0 * self._energy * FOE / + (3.0 * self._m_ejecta * M_SUN_CGS)) / KM_CGS pass @citation_wrapper('redback') From 7986f79c90542f4a9ad05d201f88862d93a7790c Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 18:04:26 +0100 Subject: [PATCH 19/28] frequencies to frequency --- redback/transient_models/tde_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/redback/transient_models/tde_models.py b/redback/transient_models/tde_models.py index c54b789d..35bf5ab3 100644 --- a/redback/transient_models/tde_models.py +++ b/redback/transient_models/tde_models.py @@ -61,13 +61,13 @@ def tde_analytical(time, redshift, l0, t_0, interaction_process=ip.Diffusion, """ frequency = kwargs['frequency'] cutoff_wavelength = kwargs.get('cutoff_wavelength', 3000) - frequency, time = calc_kcorrected_properties(frequencies=frequency, redshift=redshift, time=time) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value lbol = tde_analytical_bolometric(time=time, l0=l0, t_0=t_0, interaction_process=interaction_process) photo = photosphere(time=time, luminosity=lbol, **kwargs) sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequencies=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) + frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) flux_density = sed_1.flux_density From 3b78893fe037982511f5742cb4c1d7813a5cd537 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 18:33:04 +0100 Subject: [PATCH 20/28] homologous expansion and thin shells --- redback/transient_models/supernova_models.py | 131 ++++++++++++++++--- 1 file changed, 113 insertions(+), 18 deletions(-) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index ca791546..8099497a 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -4,9 +4,16 @@ import redback.sed as sed import redback.photosphere as photosphere from astropy.cosmology import Planck18 as cosmo # noqa -from redback.utils import calc_kcorrected_properties, citation_wrapper +from redback.utils import calc_kcorrected_properties, citation_wrapper, logger +from redback.constants import day_to_s, solar_mass, km_cgs +from inspect import isfunction import astropy.units as uu +homologous_expansion_models = ['exponential_powerlaw_bolometric', 'arnett_bolometric', + 'basic_magnetar_powered_bolometric','slsn_bolometric', + 'general_magnetar_slsn_bolometric','csm_interaction_bolometric', + 'type_1c_bolometric','type_1a_bolometric'] + def thermal_synchrotron(): """ From Margalit paper ... @@ -323,41 +330,129 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in elif kwargs['output_format'] == 'magnitude': return flux_density.to(uu.ABmag).value -@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract') -def csm_interaction(): +@citation_wrapper('redback') +def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, + **kwargs): + from redback.model_library import modules_dict # import model library in function to avoid circular dependency + base_model = kwargs['base_model'] + if isfunction(base_model): + function = base_model + elif base_model not in homologous_expansion_models: + logger.warning('{} is not implemented as a base model'.format(base_model)) + raise ValueError('Please choose a different base model') + elif isinstance(base_model, str): + function = modules_dict['supernova_models'][base_model] + else: + raise ValueError("Not a valid base model.") + v_ejecta = np.sqrt(10.0 * ek / (3.0 * mej * solar_mass)) / km_cgs + kwargs['vej'] = v_ejecta + kwargs['mej'] = mej + + lbol = function(time, interaction_process=interaction_process, **kwargs) + + return lbol + + +@citation_wrapper('redback') +def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, + **kwargs): + from redback.model_library import modules_dict # import model library in function to avoid circular dependency + base_model = kwargs['base_model'] + if isfunction(base_model): + function = base_model + elif base_model not in homologous_expansion_models: + logger.warning('{} is not implemented as a base model'.format(base_model)) + raise ValueError('Please choose a different base model') + elif isinstance(base_model, str): + function = modules_dict['supernova_models'][base_model] + else: + raise ValueError("Not a valid base model.") + v_ejecta = np.sqrt(2.0 * ek / (mej * solar_mass)) / km_cgs + kwargs['vej'] = v_ejecta + kwargs['mej'] = mej + + lbol = function(time, interaction_process=interaction_process, **kwargs) + return lbol + + +@citation_wrapper('redback') +def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, + photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + + lbol = homologous_expansion_supernova_model_bolometric(time=time, mej=mej, ek=ek, + interaction_process=interaction_process, **kwargs) + photo = photosphere(time=time, luminosity=lbol, **kwargs) + + sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency=frequency, luminosity_distance=dl) + + flux_density = sed_1.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value + +@citation_wrapper('redback') +def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, + photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + + lbol = thin_shell_supernova_model_bolometric(time=time, mej=mej, ek=ek, + interaction_process=interaction_process, **kwargs) + photo = photosphere(time=time, luminosity=lbol, **kwargs) + + sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency=frequency, luminosity_distance=dl) + + flux_density = sed_1.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value + +@citation_wrapper('redback') +def general_magnetar_slsn_bolometric(): pass -@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def csm_nickel(): +@citation_wrapper('redback') +def general_magnetar_slsn(): pass -@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def type_1a(): +@citation_wrapper('redback') +def csm_interaction_bolometric(): pass -@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def type_1c(): +@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract') +def csm_interaction(): pass -@citation_wrapper('redback') -def homologous_expansion_supernova_model_bolometric(time, **kwargs): - v_ejecta = np.sqrt(10.0 * self._energy * FOE / - (3.0 * self._m_ejecta * M_SUN_CGS)) / KM_CGS +@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') +def csm_nickel(): pass @citation_wrapper('redback') -def homologous_expansion_supernova_model(): +def csm_interaction_bolometric(): pass @citation_wrapper('redback') -def thin_shell_supernova_model_bolometric(): +def type_1a_bolometric(): pass -@citation_wrapper('redback') -def thin_shell_supernova_model(): +@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') +def type_1a(): pass @citation_wrapper('redback') -def general_magnetar_slsn(): +def type_1c_bolometric(): pass +@citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') +def type_1c(): + pass From 48cf523f7a69396a96d2a82b1fd99b0684a50500 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Tue, 8 Mar 2022 18:41:12 +0100 Subject: [PATCH 21/28] Docs for thin shell and homologous expansion --- redback/transient_models/supernova_models.py | 60 ++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index 8099497a..fc99c38f 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -333,6 +333,19 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in @citation_wrapper('redback') def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, **kwargs): + """ + Assumes homologous expansion to transform kinetic energy to ejecta velocity + + :param time: time in days in source frame + :param mej: ejecta mass in solar masses + :param ek: kinetic energy in ergs + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity + :param kwargs: Must be all the kwargs required by the specific interaction_process + e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor + 'base model' from homologous_expansion_models list + :return: bolometric_luminosity + """ from redback.model_library import modules_dict # import model library in function to avoid circular dependency base_model = kwargs['base_model'] if isfunction(base_model): @@ -356,6 +369,19 @@ def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_p @citation_wrapper('redback') def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip.Diffusion, **kwargs): + """ + Assumes thin shell ejecta to transform kinetic energy into ejecta velocity + + :param time: time in days in source frame + :param mej: ejecta mass in solar masses + :param ek: kinetic energy in ergs + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity + :param kwargs: Must be all the kwargs required by the specific interaction_process + e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor, + 'base model' from homologous_expansion_models list + :return: bolometric_luminosity + """ from redback.model_library import modules_dict # import model library in function to avoid circular dependency base_model = kwargs['base_model'] if isfunction(base_model): @@ -378,6 +404,23 @@ def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip. @citation_wrapper('redback') def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): + """ + Assumes homologous expansion to transform kinetic energy to ejecta velocity + + :param time: time in days in observer frame + :param redshift: source redshift + :param mej: ejecta mass in solar masses + :param ek: kinetic energy in ergs + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity + :param photosphere: Default is TemperatureFloor. + kwargs must have vej or relevant parameters if using different photosphere model + :param sed: Default is blackbody. + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + 'base model' from homologous_expansion_models list + :return: flux_density or magnitude depending on output_format kwarg + """ frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -399,6 +442,23 @@ def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_pr @citation_wrapper('redback') def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.Diffusion, photosphere=photosphere.TemperatureFloor, sed=sed.Blackbody, **kwargs): + """ + Assumes thin shell ejecta to transform kinetic energy into ejecta velocity + + :param time: time in days in observer frame + :param redshift: source redshift + :param mej: ejecta mass in solar masses + :param ek: kinetic energy in ergs + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity + :param photosphere: Default is TemperatureFloor. + kwargs must have vej or relevant parameters if using different photosphere model + :param sed: Default is blackbody. + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + 'base model' from homologous_expansion_models list + :return: flux_density or magnitude depending on output_format kwarg + """ frequency = kwargs['frequency'] frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value From a4894ef587ffd3f74bc60f58d294cc729da6552a Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 16:56:57 +0100 Subject: [PATCH 22/28] get csm properties utility function --- redback/utils.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/redback/utils.py b/redback/utils.py index 898346d8..b15aec62 100644 --- a/redback/utils.py +++ b/redback/utils.py @@ -32,6 +32,27 @@ def wrapper(f): return f return wrapper +def get_csm_properties(nn, eta): + csm_properties = namedtuple('csm_properties', ['AA', 'Bf', 'Br']) + filepath = f"{dirname}/tables/csm_table.txt" + ns, ss, bfs, brs, aas = np.loadtxt(filepath, delimiter=',', unpack=True) + bfs = np.reshape(bfs, (10, 30)).T + brs = np.reshape(brs, (10, 30)).T + aas = np.reshape(aas, (10, 30)).T + ns = np.unique(ns) + ss = np.unique(ss) + bf_func = RegularGridInterpolator((ss, ns), bfs) + br_func = RegularGridInterpolator((ss, ns), brs) + aa_func = RegularGridInterpolator((ss, ns), aas) + + Bf = bf_func([nn, eta])[0] + Br = br_func([nn, eta])[0] + AA = aa_func([nn, eta])[0] + + csm_properties.AA = AA + csm_properties.Bf = Bf + csm_properties.Br = Br + return csm_properties def lambda_to_nu(wavelength): """ From b27eb95b9d953aa6b491266b1e8956e786b6e1a8 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 16:57:23 +0100 Subject: [PATCH 23/28] cleaned up csm diffusion to use more derived parameters --- redback/interaction_processes.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/redback/interaction_processes.py b/redback/interaction_processes.py index 5a3583d0..b39d63d0 100644 --- a/redback/interaction_processes.py +++ b/redback/interaction_processes.py @@ -126,7 +126,7 @@ def convert_input_luminosity(self): return tau_diff, new_lums class CSMDiffusion(object): - def __init__(self, time, luminosity, kappa, csm_mass, mej, r0, eta, rho, **kwargs): + def __init__(self, time, luminosity, kappa, r_photosphere, mass_csm_threshold, csm_mass, **kwargs): """ :param time: source frame time in days :param luminosity: luminosity @@ -141,11 +141,9 @@ def __init__(self, time, luminosity, kappa, csm_mass, mej, r0, eta, rho, **kwarg self.time = time self.luminosity = luminosity self.kappa = kappa + self.r_photosphere = r_photosphere + self.mass_csm_threshold = mass_csm_threshold self.csm_mass = csm_mass * solar_mass - self.m_ejecta = mej * solar_mass - self.r0 = r0 * au_cgs - self.eta = eta - self.rho = rho self.reference = 'https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract' self.new_luminosity = [] @@ -156,18 +154,15 @@ def convert_input_luminosity(self): minimum_log_spacing = -3 # Derived parameters - # scaling constant for CSM density profile - qq = self.rho * self.r0 ** self.eta - # outer CSM shell radius - radius_csm = ((3.0 - self.eta) / (4.0 * np.pi * qq) * self.csm_mass + self.r0 ** (3.0 - self.eta)) ** (1.0 / (3.0 - self.eta)) + # photosphere radius - r_photosphere = abs((-2.0 * (1.0 - self.eta) / (3.0 * self.kappa * qq) + - radius_csm ** (1.0 - self.eta)) ** (1.0 / (1.0 - self.eta))) + r_photosphere = self.r_photosphere tau_diff = (self.kappa * self.csm_mass) / (13.8 * speed_of_light * r_photosphere) / day_to_s + # mass of the optically thick CSM (tau > 2/3). - mass_csm_threshold = np.abs(4.0 * np.pi * qq / (3.0 - self.eta) * ( - r_photosphere ** (3.0 - self.eta) - self.r0 ** (3.0 - self.eta))) + mass_csm_threshold = self.mass_csm_threshold + beta = 4. * np.pi ** 3. / 9. t0 = self.kappa * (mass_csm_threshold) / (beta * speed_of_light * r_photosphere) / day_to_s From a6bc883d416313eb3ff28f2ab961818f34da3a41 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 16:57:41 +0100 Subject: [PATCH 24/28] units of velocity fix --- redback/photosphere.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/redback/photosphere.py b/redback/photosphere.py index c937e5a7..d790e849 100644 --- a/redback/photosphere.py +++ b/redback/photosphere.py @@ -10,7 +10,7 @@ def __init__(self, time, luminosity, vej, temperature_floor, **kwargs): :param time: source frame time in days :param luminosity: luminosity - :param vej: ejecta velocity in cgs + :param vej: ejecta velocity in km/s :param temperature_floor: floor temperature in kelvin """ self.time = time @@ -117,7 +117,7 @@ def __init__(self, time, luminosity, mej, vej, kappa, envelope_slope=10, **kwarg :param time: time in source frame in days :param luminosity: luminosity :param mej: ejecta mass - :param vej: ejecta velocity in km + :param vej: ejecta velocity in km/s :param kappa: opacity :param envelope_slope: envelope slope, default = 10 """ From 6e9f70c378fc462f763176c6f7f2db4465a3b26b Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 16:57:53 +0100 Subject: [PATCH 25/28] csm and csmni models --- redback/transient_models/supernova_models.py | 257 ++++++++++++++++--- 1 file changed, 227 insertions(+), 30 deletions(-) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index fc99c38f..10cf5d9c 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -4,10 +4,11 @@ import redback.sed as sed import redback.photosphere as photosphere from astropy.cosmology import Planck18 as cosmo # noqa -from redback.utils import calc_kcorrected_properties, citation_wrapper, logger -from redback.constants import day_to_s, solar_mass, km_cgs +from redback.utils import calc_kcorrected_properties, citation_wrapper, logger, get_csm_properties +from redback.constants import day_to_s, solar_mass, km_cgs, au_cgs from inspect import isfunction import astropy.units as uu +from collections import namedtuple homologous_expansion_models = ['exponential_powerlaw_bolometric', 'arnett_bolometric', 'basic_magnetar_powered_bolometric','slsn_bolometric', @@ -32,7 +33,7 @@ def exponential_powerlaw_bolometric(time, lbol_0, alpha_1, alpha_2, tpeak_d, int :param alpha_2: second exponent :param tpeak_d: peak time in days :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, mej (solar masses), vej (km/s), temperature_floor :return: bolometric_luminosity @@ -57,7 +58,7 @@ def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, :param alpha_2: second exponent :param tpeak_d: peak time in days :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -108,7 +109,7 @@ def arnett_bolometric(time, f_nickel, mej, interaction_process=ip.Diffusion, **k :param f_nickel: fraction of nickel mass :param mej: total ejecta mass in solar masses :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor :return: bolometric_luminosity @@ -130,7 +131,7 @@ def arnett(time, redshift, f_nickel, mej, interaction_process=ip.Diffusion, :param f_nickel: fraction of nickel mass :param mej: total ejecta mass in solar masses :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -179,7 +180,7 @@ def basic_magnetar_powered_bolometric(time, p0, bp, mass_ns, theta_pb,interactio :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor :return: bolometric_luminosity @@ -201,7 +202,7 @@ def basic_magnetar_powered(time, redshift, p0, bp, mass_ns, theta_pb, interactio :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -238,7 +239,7 @@ def slsn_bolometric(time, p0, bp, mass_ns, theta_pb,interaction_process=ip.Diffu :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor :return: bolometric_luminosity @@ -260,7 +261,7 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is CutoffBlackbody. @@ -298,7 +299,7 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in :param mass_ns: mass of neutron star in solar masses :param theta_pb: angle between spin and magnetic field axes :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -340,7 +341,7 @@ def homologous_expansion_supernova_model_bolometric(time, mej, ek, interaction_p :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor 'base model' from homologous_expansion_models list @@ -376,7 +377,7 @@ def thin_shell_supernova_model_bolometric(time, mej, ek, interaction_process=ip. :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param kwargs: Must be all the kwargs required by the specific interaction_process e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor, 'base model' from homologous_expansion_models list @@ -412,7 +413,7 @@ def homologous_expansion_supernova_model(time, redshift, mej, ek, interaction_pr :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -450,7 +451,7 @@ def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.D :param mej: ejecta mass in solar masses :param ek: kinetic energy in ergs :param interaction_process: Default is Diffusion. - Can also be None in which case the output is just the raw engine luminosity + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. :param photosphere: Default is TemperatureFloor. kwargs must have vej or relevant parameters if using different photosphere model :param sed: Default is blackbody. @@ -477,29 +478,215 @@ def thin_shell_supernova_model(time, redshift, mej, ek, interaction_process=ip.D elif kwargs['output_format'] == 'magnitude': return flux_density.to(uu.ABmag).value -@citation_wrapper('redback') -def general_magnetar_slsn_bolometric(): - pass -@citation_wrapper('redback') -def general_magnetar_slsn(): - pass +def _csm_engine(time, mej, csm_mass, vej, eta, rho, kappa, r0, **kwargs): + """ + :param time: time in days in source frame + :param mej: ejecta mass in solar masses + :param csm_mass: csm mass in solar masses + :param vej: ejecta velocity in km/s + :param eta: csm density profile exponent + :param rho: csm density profile amplitude + :param kappa: opacity + :param r0: radius of csm shell in AU + :param kwargs: + efficiency: in converting between kinetic energy and luminosity, default 0.5 + delta: default 1, + nn: default 12, + :return: named tuple with 'lbol','r_photosphere' 'mass_csm_threshold' + """ + mej = mej * solar_mass + csm_mass = csm_mass * solar_mass + r0 = r0 * au_cgs + vej = vej * km_cgs + Esn = 3. * vej ** 2 * mej / 10. + ti = 1. + + delta = kwargs.get('delta', 1) + nn = kwargs.get('nn', 12) + efficiency = kwargs.get('efficiency', 0.5) + + csm_properties = get_csm_properties(nn, eta) + AA = csm_properties.AA + Bf = csm_properties.Bf + Br = csm_properties.Br + + # Derived parameters + # scaling constant for CSM density profile + qq = rho * r0 ** eta + # outer CSM shell radius + radius_csm = ((3.0 - eta) / (4.0 * np.pi * qq) * csm_mass + r0 ** (3.0 - eta)) ** ( + 1.0 / (3.0 - eta)) + # photosphere radius + r_photosphere = abs((-2.0 * (1.0 - eta) / (3.0 * kappa * qq) + + radius_csm ** (1.0 - eta)) ** (1.0 / (1.0 - eta))) + + # mass of the optically thick CSM (tau > 2/3). + mass_csm_threshold = np.abs(4.0 * np.pi * qq / (3.0 - eta) * ( + r_photosphere ** (3.0 - eta) - r0 ** (3.0 - eta))) + + # g**n is scaling parameter for ejecta density profile + g_n = (1.0 / (4.0 * np.pi * (nn - delta)) * ( + 2.0 * (5.0 - delta) * (nn - 5.0) * Esn) ** ((nn - 3.) / 2.0) / ( + (3.0 - delta) * (nn - 3.0) * mej) ** ((nn - 5.0) / 2.0)) + + # time at which shock breaks out of optically thick CSM - forward shock + t_FS = (abs((3.0 - eta) * qq ** ((3.0 - nn) / (nn - eta)) * ( + AA * g_n) ** ((eta - 3.0) / (nn - eta)) / + (4.0 * np.pi * Bf ** (3.0 - eta))) ** ( + (nn - eta) / ((nn - 3.0) * (3.0 - eta))) * (mass_csm_threshold) ** ( + (nn - eta) / ((nn - 3.0) * (3.0 - eta)))) + + # time at which reverse shock sweeps up all ejecta - reverse shock + t_RS = (vej / (Br * (AA * g_n / qq) ** ( + 1.0 / (nn - eta))) * + (1.0 - (3.0 - nn) * mej / + (4.0 * np.pi * vej ** + (3.0 - nn) * g_n)) ** (1.0 / (3.0 - nn))) ** ( + (nn - eta) / (eta - 3.0)) + + mask_1 = t_FS - time * day_to_s > 0 + mask_2 = t_RS - time * day_to_s > 0 + + lbol = efficiency * (2.0 * np.pi / (nn - eta) ** 3 * g_n ** ((5.0 - eta) / (nn - eta)) * qq ** + ((nn - 5.0) / (nn - eta)) * (nn - 3.0) ** 2 * (nn - 5.0) * Bf ** (5.0 - eta) * AA ** + ((5.0 - eta) / (nn - eta)) * (time * day_to_s + ti) ** + ((2.0 * nn + 6.0 * eta - nn * eta - 15.) / + (nn - eta)) + 2.0 * np.pi * (AA * g_n / qq) ** + ((5.0 - nn) / (nn - eta)) * Br ** (5.0 - nn) * g_n * ((3.0 - eta) / (nn - eta)) ** 3 * ( + time * day_to_s + ti) ** + ((2.0 * nn + 6.0 * eta - nn * eta - 15.0) / (nn - eta))) + + lbol[~mask_1] = 0 + lbol[~mask_2] = 0 + + csm_output = namedtuple('csm_output', ['lbol', 'r_photosphere', 'mass_csm_threshold']) + csm_output.lbol = lbol + csm_output.r_photosphere = r_photosphere + csm_output.mass_csm_threshold = mass_csm_threshold + return csm_output + @citation_wrapper('redback') -def csm_interaction_bolometric(): - pass +def csm_interaction_bolometric(time, mej, csm_mass, vej, eta, rho, kappa, r0, + interaction_process=ip.CSMDiffusion, **kwargs): + """ + :param time: time in days in source frame + :param mej: ejecta mass in solar masses + :param csm_mass: csm mass in solar masses + :param vej: ejecta velocity in km/s + :param eta: csm density profile exponent + :param rho: csm density profile amplitude + :param kappa: opacity + :param r0: radius of csm shell in AU + :param interaction_process: Default is CSMDiffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + :param kwargs: + efficiency: in converting between kinetic energy and luminosity, default 0.5 + delta: default 1, + nn: default 12, + If interaction process is different kwargs must include other keyword arguments that are required. + :return: bolometric_luminosity + """ + csm_output = _csm_engine(time=time, mej=mej, csm_mass=csm_mass, vej=vej, + eta=eta, rho=rho, kappa=kappa, r0=r0, **kwargs) + lbol = csm_output.lbol + r_photosphere = csm_output.r_photosphere + mass_csm_threshold = csm_output.mass_csm_threshold + + if interaction_process is not None: + interaction_class = interaction_process(time=time, luminosity=lbol, + kappa=kappa, r_photosphere=r_photosphere, + mass_csm_threshold=mass_csm_threshold, csm_mass=csm_mass, **kwargs) + lbol = interaction_class.new_luminosity + return lbol + @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2013ApJ...773...76C/abstract') -def csm_interaction(): - pass +def csm_interaction(time, redshift, mej, csm_mass, vej, eta, rho, kappa, r0, + interaction_process=ip.CSMDiffusion, photosphere=photosphere.TemperatureFloor, + sed=sed.Blackbody, **kwargs): + """ + :param time: time in days in observer frame + :param redshift: source redshift + :param mej: ejecta mass in solar masses + :param csm_mass: csm mass in solar masses + :param vej: ejecta velocity in km/s + :param eta: csm density profile exponent + :param rho: csm density profile amplitude + :param kappa: opacity + :param r0: radius of csm shell in AU + :param interaction_process: Default is CSMDiffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + :param photosphere: Default is TemperatureFloor. + kwargs must have vej or relevant parameters if using different photosphere model + :param sed: Default is blackbody. + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa_gamma, temperature_floor + 'base model' from homologous_expansion_models list + :return: flux_density or magnitude depending on output_format kwarg + """ + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + + lbol = csm_interaction_bolometric(time=time, mej=mej, csm_mass=csm_mass, vej=vej, eta=eta, + rho=rho, kappa=kappa, r0=r0, interaction_process=interaction_process, **kwargs) + + photo = photosphere(time=time, luminosity=lbol, vej=vej, **kwargs) + + sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency=frequency, luminosity_distance=dl) + + flux_density = sed_1.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def csm_nickel(): - pass +def csm_nickel(time, redshift, mej, f_nickel, csm_mass, ek, eta, rho, kappa, r0, **kwargs): + """ + Assumes csm and nickel engine with homologous expansion + + :param time: time in days in observer frame + :param redshift: source redshift + :param mej: ejecta mass in solar masses + :param csm_mass: csm mass in solar masses + :param ek: kinetic energy in ergs + :param eta: csm density profile exponent + :param rho: csm density profile amplitude + :param kappa: opacity + :param r0: radius of csm shell in AU + :param kwargs: kappa_gamma, temperature_floor, and any kwarg to + change any other input physics/parameters from default. + :return: flux_density or magnitude depending on output_format kwarg + """ + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + + vej = np.sqrt(2.0 * ek / (mej * solar_mass)) / km_cgs + kwargs['vej'] = vej + nickel_lbol = arnett_bolometric(time=time, f_nickel=f_nickel, + mej=mej, interaction_process=ip.Diffusion, **kwargs) + csm_lbol = csm_interaction_bolometric(time=time, mej=mej, csm_mass=csm_mass, eta=eta, + rho=rho, kappa=kappa, r0=r0, interaction_process=ip.CSMDiffusion, **kwargs) + lbol = nickel_lbol + csm_lbol + + photo = photosphere.TemperatureFloor(time=time, luminosity=lbol, vej=vej, **kwargs) + + sed_1 = sed.Blackbody(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency=frequency, luminosity_distance=dl) + + flux_density = sed_1.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value -@citation_wrapper('redback') -def csm_interaction_bolometric(): - pass @citation_wrapper('redback') def type_1a_bolometric(): @@ -516,3 +703,13 @@ def type_1c_bolometric(): @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') def type_1c(): pass + +@citation_wrapper('redback') +def general_magnetar_slsn_bolometric(): + pass + +@citation_wrapper('redback') +def general_magnetar_slsn(): + pass + + From 44b128f5ea7f2d18bdb8035a0fc74b9683161a4d Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 17:04:51 +0100 Subject: [PATCH 26/28] Warning to install from source --- docs/installation.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/installation.txt b/docs/installation.txt index 742093be..7fe4a76c 100644 --- a/docs/installation.txt +++ b/docs/installation.txt @@ -24,8 +24,10 @@ Installation This will install all requirements for running :code:`redback` for general transient fitting, including bilby and our default sampler `dynesty `_. Other samplers will need to be installed via pip or the appropriate means. Please look through :code:`redback` `fitting documentation `_ for what samplers are available. +Currently :code:`redback` is going significant development, and we can not guarantee that any `pip` or `conda` releases will be completely up to date. +Therefore, for the near future, we recommend installing :code:`redback` from source. -Install redback from source +Install :code:`redback` from source ------------------------- :code:`redback` is developed and tested with Python 3.7+. In the From 51e7401f8eddaef0821bbd00f1c446537785f8d5 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 18:19:59 +0100 Subject: [PATCH 27/28] CSM interpolation table --- redback/tables/csm_table.txt | 300 +++++++++++++++++++++++++++++++++++ 1 file changed, 300 insertions(+) create mode 100644 redback/tables/csm_table.txt diff --git a/redback/tables/csm_table.txt b/redback/tables/csm_table.txt new file mode 100644 index 00000000..d52fa553 --- /dev/null +++ b/redback/tables/csm_table.txt @@ -0,0 +1,300 @@ +0.0,6.0,1.255564,0.906057,2.422657 +0.0,6.275862068965517,1.224261,0.916353,1.970811 +0.0,6.551724137931035,1.203396,0.924596,1.612987 +0.0,6.827586206896552,1.188402,0.931348,1.344976 +0.0,7.103448275862069,1.176956,0.936982,1.133469 +0.0,7.379310344827586,1.168122,0.941747,0.973674 +0.0,7.655172413793103,1.161069,0.945837,0.843831 +0.0,7.931034482758621,1.155227,0.949398,0.737170 +0.0,8.206896551724139,1.150336,0.952522,0.650604 +0.0,8.482758620689655,1.146213,0.955352,0.552565 +0.0,8.758620689655173,1.142580,0.957675,0.516337 +0.0,9.03448275862069,1.139308,0.959901,0.468747 +0.0,9.310344827586206,1.136488,0.961901,0.424121 +0.0,9.586206896551724,1.134195,0.963680,0.384980 +0.0,9.862068965517242,1.132012,0.965301,0.351938 +0.0,10.137931034482758,1.130394,0.966895,0.306736 +0.0,10.413793103448276,1.128640,0.967815,0.291763 +0.0,10.689655172413794,1.127025,0.969136,0.270826 +0.0,10.96551724137931,1.125595,0.970348,0.251014 +0.0,11.241379310344827,1.124306,0.971468,0.232984 +0.0,11.517241379310345,1.123132,0.972509,0.216703 +0.0,11.793103448275861,1.122056,0.973479,0.202010 +0.0,12.068965517241379,1.121065,0.974386,0.188731 +0.0,12.344827586206897,1.120148,0.975236,0.176704 +0.0,12.620689655172413,1.119296,0.976035,0.165785 +0.0,12.89655172413793,1.118501,0.976787,0.155847 +0.0,13.172413793103448,1.117757,0.977496,0.146780 +0.0,13.448275862068964,1.117038,0.978166,0.138597 +0.0,13.724137931034482,1.116364,0.978799,0.131090 +0.0,14.0,1.115743,0.979399,0.124179 +0.2222222222222222,6.0,1.265533,0.910434,2.176640 +0.2222222222222222,6.275862068965517,1.233603,0.920151,1.759751 +0.2222222222222222,6.551724137931035,1.212489,0.927946,1.430740 +0.2222222222222222,6.827586206896552,1.197025,0.934345,1.192416 +0.2222222222222222,7.103448275862069,1.185267,0.939683,1.008465 +0.2222222222222222,7.379310344827586,1.175887,0.944227,0.867156 +0.2222222222222222,7.655172413793103,1.168560,0.948091,0.753864 +0.2222222222222222,7.931034482758621,1.162757,0.951510,0.655132 +0.2222222222222222,8.206896551724139,1.157844,0.954546,0.551150 +0.2222222222222222,8.482758620689655,1.153454,0.957318,0.485463 +0.2222222222222222,8.758620689655173,1.149662,0.959425,0.460698 +0.2222222222222222,9.03448275862069,1.146377,0.961538,0.414590 +0.2222222222222222,9.310344827586206,1.143533,0.963400,0.374990 +0.2222222222222222,9.586206896551724,1.141455,0.965137,0.337521 +0.2222222222222222,9.862068965517242,1.139097,0.966225,0.308853 +0.2222222222222222,10.137931034482758,1.137111,0.967746,0.283670 +0.2222222222222222,10.413793103448276,1.135358,0.969119,0.260480 +0.2222222222222222,10.689655172413794,1.133782,0.970377,0.239822 +0.2222222222222222,10.96551724137931,1.132349,0.971536,0.221480 +0.2222222222222222,11.241379310344827,1.131036,0.972611,0.205166 +0.2222222222222222,11.517241379310345,1.129828,0.973612,0.190592 +0.2222222222222222,11.793103448275861,1.128616,0.974546,0.178230 +0.2222222222222222,12.068965517241379,1.127618,0.975420,0.166206 +0.2222222222222222,12.344827586206897,1.126670,0.976241,0.155516 +0.2222222222222222,12.620689655172413,1.125755,0.977013,0.146015 +0.2222222222222222,12.89655172413793,1.124902,0.977739,0.137354 +0.2222222222222222,13.172413793103448,1.124104,0.978425,0.129439 +0.2222222222222222,13.448275862068964,1.123358,0.979073,0.122188 +0.2222222222222222,13.724137931034482,1.122658,0.979687,0.115527 +0.2222222222222222,14.0,1.122000,0.980268,0.109394 +0.4444444444444444,6.0,1.276282,0.915101,1.943555 +0.4444444444444444,6.275862068965517,1.243829,0.924205,1.556329 +0.4444444444444444,6.551724137931035,1.222099,0.931533,1.269973 +0.4444444444444444,6.827586206896552,1.206249,0.937540,1.059490 +0.4444444444444444,7.103448275862069,1.194383,0.942597,0.889087 +0.4444444444444444,7.379310344827586,1.184808,0.946907,0.754087 +0.4444444444444444,7.655172413793103,1.177047,0.950568,0.659833 +0.4444444444444444,7.931034482758621,1.170908,0.953829,0.555989 +0.4444444444444444,8.206896551724139,1.165972,0.956800,0.482497 +0.4444444444444444,8.482758620689655,1.161459,0.959059,0.452310 +0.4444444444444444,8.758620689655173,1.157606,0.961312,0.403738 +0.4444444444444444,9.03448275862069,1.154651,0.963263,0.361226 +0.4444444444444444,9.310344827586206,1.151578,0.964470,0.330717 +0.4444444444444444,9.586206896551724,1.149074,0.966252,0.299055 +0.4444444444444444,9.862068965517242,1.146864,0.967820,0.271406 +0.4444444444444444,10.137931034482758,1.144867,0.969238,0.247581 +0.4444444444444444,10.413793103448276,1.142935,0.970536,0.227861 +0.4444444444444444,10.689655172413794,1.141335,0.971731,0.209327 +0.4444444444444444,10.96551724137931,1.139790,0.972837,0.193557 +0.4444444444444444,11.241379310344827,1.138377,0.973865,0.179509 +0.4444444444444444,11.517241379310345,1.137080,0.974823,0.166929 +0.4444444444444444,11.793103448275861,1.135888,0.975719,0.155616 +0.4444444444444444,12.068965517241379,1.134788,0.976558,0.145400 +0.4444444444444444,12.344827586206897,1.133772,0.977347,0.136143 +0.4444444444444444,12.620689655172413,1.132831,0.978089,0.127725 +0.4444444444444444,12.89655172413793,1.131958,0.978789,0.120042 +0.4444444444444444,13.172413793103448,1.131149,0.979450,0.113004 +0.4444444444444444,13.448275862068964,1.130405,0.980076,0.106523 +0.4444444444444444,13.724137931034482,1.129746,0.980671,0.100475 +0.4444444444444444,14.0,1.128973,0.981236,0.095255 +0.6666666666666666,6.0,1.287803,0.920077,1.728850 +0.6666666666666666,6.275862068965517,1.254881,0.928547,1.372745 +0.6666666666666666,6.551724137931035,1.232803,0.935361,1.114627 +0.6666666666666666,6.827586206896552,1.216370,0.940990,0.929476 +0.6666666666666666,7.103448275862069,1.204149,0.945670,0.783812 +0.6666666666666666,7.379310344827586,1.194582,0.949746,0.666801 +0.6666666666666666,7.655172413793103,1.186942,0.953240,0.554187 +0.6666666666666666,7.931034482758621,1.180452,0.956416,0.478146 +0.6666666666666666,8.206896551724139,1.174975,0.958840,0.443915 +0.6666666666666666,8.482758620689655,1.170519,0.961205,0.392934 +0.6666666666666666,8.758620689655173,1.166721,0.962527,0.360783 +0.6666666666666666,9.03448275862069,1.163444,0.964665,0.317867 +0.6666666666666666,9.310344827586206,1.160284,0.966468,0.286406 +0.6666666666666666,9.586206896551724,1.157805,0.968074,0.257760 +0.6666666666666666,9.862068965517242,1.155402,0.969530,0.234593 +0.6666666666666666,10.137931034482758,1.153254,0.970862,0.214473 +0.6666666666666666,10.413793103448276,1.151329,0.972087,0.196818 +0.6666666666666666,10.689655172413794,1.149603,0.973220,0.181210 +0.6666666666666666,10.96551724137931,1.148061,0.974270,0.167297 +0.6666666666666666,11.241379310344827,1.146343,0.975248,0.155781 +0.6666666666666666,11.517241379310345,1.145354,0.976160,0.143899 +0.6666666666666666,11.793103448275861,1.144053,0.977013,0.134243 +0.6666666666666666,12.068965517241379,1.142830,0.977812,0.125549 +0.6666666666666666,12.344827586206897,1.141776,0.978563,0.117502 +0.6666666666666666,12.620689655172413,1.140862,0.979270,0.110079 +0.6666666666666666,12.89655172413793,1.139697,0.979935,0.103812 +0.6666666666666666,13.172413793103448,1.139276,0.980564,0.097180 +0.6666666666666666,13.448275862068964,1.136756,0.981158,0.092091 +0.6666666666666666,13.724137931034482,1.137323,0.981720,0.086976 +0.6666666666666666,14.0,1.136925,0.982252,0.082034 +0.8888888888888888,6.0,1.300164,0.925380,1.522374 +0.8888888888888888,6.275862068965517,1.267074,0.933202,1.194479 +0.8888888888888888,6.551724137931035,1.244420,0.939534,0.967964 +0.8888888888888888,6.827586206896552,1.227784,0.944676,0.800617 +0.8888888888888888,7.103448275862069,1.215355,0.949143,0.674534 +0.8888888888888888,7.379310344827586,1.205195,0.952838,0.559115 +0.8888888888888888,7.655172413793103,1.197427,0.956159,0.477563 +0.8888888888888888,7.931034482758621,1.190291,0.958868,0.435695 +0.8888888888888888,8.206896551724139,1.185505,0.961297,0.381248 +0.8888888888888888,8.482758620689655,1.180882,0.963023,0.341271 +0.8888888888888888,8.758620689655173,1.176893,0.965100,0.301193 +0.8888888888888888,9.03448275862069,1.173293,0.966923,0.269194 +0.8888888888888888,9.310344827586206,1.170007,0.968560,0.243148 +0.8888888888888888,9.586206896551724,1.166824,0.970046,0.220793 +0.8888888888888888,9.862068965517242,1.164965,0.971405,0.199601 +0.8888888888888888,10.137931034482758,1.162914,0.972654,0.182257 +0.8888888888888888,10.413793103448276,1.160943,0.973805,0.167212 +0.8888888888888888,10.689655172413794,1.158739,0.974870,0.154931 +0.8888888888888888,10.96551724137931,1.157478,0.975858,0.142472 +0.8888888888888888,11.241379310344827,1.155261,0.976776,0.133651 +0.8888888888888888,11.517241379310345,1.154650,0.977629,0.122614 +0.8888888888888888,11.793103448275861,1.153002,0.978423,0.115628 +0.8888888888888888,12.068965517241379,1.152222,0.979163,0.107028 +0.8888888888888888,12.344827586206897,1.151143,0.979864,0.100462 +0.8888888888888888,12.620689655172413,1.149810,0.980626,0.094690 +0.8888888888888888,12.89655172413793,1.149219,0.980821,0.088282 +0.8888888888888888,13.172413793103448,1.148362,0.981372,0.082967 +0.8888888888888888,13.448275862068964,1.147528,0.981945,0.077864 +0.8888888888888888,13.724137931034482,1.146481,0.982418,0.075373 +0.8888888888888888,14.0,1.146016,0.982891,0.070618 +1.1111111111111112,6.0,1.313240,0.931102,1.328203 +1.1111111111111112,6.275862068965517,1.280046,0.938211,1.035931 +1.1111111111111112,6.551724137931035,1.257607,0.944001,0.831155 +1.1111111111111112,6.827586206896552,1.240125,0.948644,0.671422 +1.1111111111111112,7.103448275862069,1.227565,0.952698,0.563969 +1.1111111111111112,7.379310344827586,1.217544,0.956193,0.475319 +1.1111111111111112,7.655172413793103,1.209437,0.959006,0.428631 +1.1111111111111112,7.931034482758621,1.202823,0.961434,0.366153 +1.1111111111111112,8.206896551724139,1.197155,0.963792,0.318823 +1.1111111111111112,8.482758620689655,1.192338,0.965850,0.281275 +1.1111111111111112,8.758620689655173,1.188195,0.967689,0.249724 +1.1111111111111112,9.03448275862069,1.184391,0.969349,0.224987 +1.1111111111111112,9.310344827586206,1.181124,0.970858,0.203849 +1.1111111111111112,9.586206896551724,1.178492,0.972235,0.184604 +1.1111111111111112,9.862068965517242,1.176000,0.973496,0.166565 +1.1111111111111112,10.137931034482758,1.173819,0.974654,0.151956 +1.1111111111111112,10.413793103448276,1.171755,0.975720,0.139862 +1.1111111111111112,10.689655172413794,1.169903,0.976702,0.128972 +1.1111111111111112,10.96551724137931,1.167979,0.977620,0.125786 +1.1111111111111112,11.241379310344827,1.166679,0.978587,0.110248 +1.1111111111111112,11.517241379310345,1.165251,0.979127,0.100821 +1.1111111111111112,11.793103448275861,1.163904,0.979915,0.092392 +1.1111111111111112,12.068965517241379,1.162700,0.980419,0.095047 +1.1111111111111112,12.344827586206897,1.161624,0.981071,0.082609 +1.1111111111111112,12.620689655172413,1.160467,0.981700,0.076763 +1.1111111111111112,12.89655172413793,1.159525,0.982167,0.077382 +1.1111111111111112,13.172413793103448,1.158685,0.982771,0.067752 +1.1111111111111112,13.448275862068964,1.157857,0.983171,0.065358 +1.1111111111111112,13.724137931034482,1.157051,0.983668,0.061026 +1.1111111111111112,14.0,1.156313,0.986211,0.048402 +1.3333333333333333,6.0,1.328001,0.937219,1.131763 +1.3333333333333333,6.275862068965517,1.294944,0.943562,0.878302 +1.3333333333333333,6.551724137931035,1.271693,0.948624,0.699505 +1.3333333333333333,6.827586206896552,1.253823,0.952877,0.579154 +1.3333333333333333,7.103448275862069,1.240536,0.956154,0.518465 +1.3333333333333333,7.379310344827586,1.230621,0.960426,0.382656 +1.3333333333333333,7.655172413793103,1.223061,0.962692,0.334033 +1.3333333333333333,7.931034482758621,1.216049,0.964965,0.293938 +1.3333333333333333,8.206896551724139,1.209961,0.967019,0.257453 +1.3333333333333333,8.482758620689655,1.205025,0.968871,0.232910 +1.3333333333333333,8.758620689655173,1.200938,0.970545,0.207735 +1.3333333333333333,9.03448275862069,1.197096,0.972064,0.182146 +1.3333333333333333,9.310344827586206,1.193982,0.973445,0.170255 +1.3333333333333333,9.586206896551724,1.190929,0.974704,0.160814 +1.3333333333333333,9.862068965517242,1.188577,0.975867,0.136135 +1.3333333333333333,10.137931034482758,1.186257,0.977095,0.123615 +1.3333333333333333,10.413793103448276,1.184155,0.978443,0.103013 +1.3333333333333333,10.689655172413794,1.182246,0.978462,0.107668 +1.3333333333333333,10.96551724137931,1.180506,0.979261,0.098613 +1.3333333333333333,11.241379310344827,1.178816,0.980368,0.090810 +1.3333333333333333,11.517241379310345,1.177449,0.980902,0.081622 +1.3333333333333333,11.793103448275861,1.176099,0.981528,0.076342 +1.3333333333333333,12.068965517241379,1.174852,0.982133,0.071046 +1.3333333333333333,12.344827586206897,1.173691,0.982773,0.065466 +1.3333333333333333,12.620689655172413,1.172613,0.982988,0.065101 +1.3333333333333333,12.89655172413793,1.171607,0.983761,0.058204 +1.3333333333333333,13.172413793103448,1.170654,0.984174,0.058137 +1.3333333333333333,13.448275862068964,1.169788,0.984164,0.066268 +1.3333333333333333,13.724137931034482,1.168960,0.985044,0.050503 +1.3333333333333333,14.0,1.168182,0.985819,0.041741 +1.5555555555555554,6.0,1.344585,0.943710,0.960476 +1.5555555555555554,6.275862068965517,1.309117,0.948967,0.803768 +1.5555555555555554,6.551724137931035,1.287028,0.953321,0.638809 +1.5555555555555554,6.827586206896552,1.269438,0.957846,0.485648 +1.5555555555555554,7.103448275862069,1.256959,0.962133,0.346817 +1.5555555555555554,7.379310344827586,1.246587,0.964436,0.305980 +1.5555555555555554,7.655172413793103,1.238157,0.966696,0.265271 +1.5555555555555554,7.931034482758621,1.230987,0.968753,0.230565 +1.5555555555555554,8.206896551724139,1.225268,0.970610,0.202992 +1.5555555555555554,8.482758620689655,1.219367,0.972283,0.194503 +1.5555555555555554,8.758620689655173,1.215850,0.973798,0.160196 +1.5555555555555554,9.03448275862069,1.212052,0.975229,0.154127 +1.5555555555555554,9.310344827586206,1.208714,0.975832,0.137845 +1.5555555555555554,9.586206896551724,1.205741,0.976899,0.124075 +1.5555555555555554,9.862068965517242,1.202999,0.977892,0.122915 +1.5555555555555554,10.137931034482758,1.200692,0.978813,0.099714 +1.5555555555555554,10.413793103448276,1.198531,0.980391,0.082680 +1.5555555555555554,10.689655172413794,1.196569,0.981007,0.078452 +1.5555555555555554,10.96551724137931,1.194777,0.981675,0.073346 +1.5555555555555554,11.241379310344827,1.192798,0.982318,0.074961 +1.5555555555555554,11.517241379310345,1.191628,0.982924,0.063598 +1.5555555555555554,11.793103448275861,1.190236,0.983490,0.059464 +1.5555555555555554,12.068965517241379,1.188946,0.984027,0.055996 +1.5555555555555554,12.344827586206897,1.187752,0.984709,0.052577 +1.5555555555555554,12.620689655172413,1.186638,0.984577,0.053147 +1.5555555555555554,12.89655172413793,1.185595,0.985750,0.041214 +1.5555555555555554,13.172413793103448,1.184435,0.985973,0.045235 +1.5555555555555554,13.448275862068964,1.183718,0.985862,0.044005 +1.5555555555555554,13.724137931034482,1.182863,0.986240,0.039809 +1.5555555555555554,14.0,1.182058,0.986612,0.040205 +1.7777777777777777,6.0,1.359841,0.950834,0.771436 +1.7777777777777777,6.275862068965517,1.327202,0.955663,0.585432 +1.7777777777777777,6.551724137931035,1.304325,0.959589,0.501989 +1.7777777777777777,6.827586206896552,1.287483,0.962928,0.382406 +1.7777777777777777,7.103448275862069,1.274438,0.966914,0.273171 +1.7777777777777777,7.379310344827586,1.264052,0.969184,0.233236 +1.7777777777777777,7.655172413793103,1.255580,0.971464,0.191511 +1.7777777777777777,7.931034482758621,1.248534,0.973258,0.166149 +1.7777777777777777,8.206896551724139,1.242580,0.973668,0.174614 +1.7777777777777777,8.482758620689655,1.237481,0.975029,0.159925 +1.7777777777777777,8.758620689655173,1.233065,0.976398,0.135567 +1.7777777777777777,9.03448275862069,1.229202,0.977556,0.118098 +1.7777777777777777,9.310344827586206,1.225794,0.978610,0.113834 +1.7777777777777777,9.586206896551724,1.222765,0.979545,0.095872 +1.7777777777777777,9.862068965517242,1.220055,0.981455,0.073715 +1.7777777777777777,10.137931034482758,1.217616,0.982157,0.068767 +1.7777777777777777,10.413793103448276,1.215315,0.983065,0.063647 +1.7777777777777777,10.689655172413794,1.213402,0.983707,0.054914 +1.7777777777777777,10.96551724137931,1.211570,0.983304,0.065254 +1.7777777777777777,11.241379310344827,1.209889,0.983891,0.056388 +1.7777777777777777,11.517241379310345,1.208344,0.984441,0.053210 +1.7777777777777777,11.793103448275861,1.206917,0.984948,0.049597 +1.7777777777777777,12.068965517241379,1.205595,0.986347,0.037645 +1.7777777777777777,12.344827586206897,1.204367,0.985883,0.043848 +1.7777777777777777,12.620689655172413,1.203224,0.987172,0.033026 +1.7777777777777777,12.89655172413793,1.202157,0.987481,0.032000 +1.7777777777777777,13.172413793103448,1.201159,0.987083,0.038950 +1.7777777777777777,13.448275862068964,1.200223,0.988299,0.026847 +1.7777777777777777,13.724137931034482,1.199343,0.988548,0.026177 +1.7777777777777777,14.0,1.198516,0.988103,0.030068 +2.0,6.0,1.376441,0.958424,0.640501 +2.0,6.275862068965517,1.342981,0.962421,0.540211 +2.0,6.551724137931035,1.323394,0.965730,0.385833 +2.0,6.827586206896552,1.306952,0.968462,0.315803 +2.0,7.103448275862069,1.294265,0.970789,0.247060 +2.0,7.379310344827586,1.284048,0.972838,0.206567 +2.0,7.655172413793103,1.275678,0.974592,0.173838 +2.0,7.931034482758621,1.268692,0.976126,0.149494 +2.0,8.206896551724139,1.262497,0.977495,0.150561 +2.0,8.482758620689655,1.257684,0.978716,0.116013 +2.0,8.758620689655173,1.253269,0.979808,0.103820 +2.0,9.03448275862069,1.249399,0.980788,0.093314 +2.0,9.310344827586206,1.245978,0.981686,0.083408 +2.0,9.586206896551724,1.242933,0.982475,0.072907 +2.0,9.862068965517242,1.239901,0.983234,0.076206 +2.0,10.137931034482758,1.237536,0.983919,0.066222 +2.0,10.413793103448276,1.235518,0.984558,0.056680 +2.0,10.689655172413794,1.233490,0.985107,0.055730 +2.0,10.96551724137931,1.231636,0.985665,0.045244 +2.0,11.241379310344827,1.229913,0.986170,0.046977 +2.0,11.517241379310345,1.228369,0.986657,0.041696 +2.0,11.793103448275861,1.226922,0.987117,0.037656 +2.0,12.068965517241379,1.225580,0.987487,0.036378 +2.0,12.344827586206897,1.224334,0.987875,0.031703 +2.0,12.620689655172413,1.223172,0.988243,0.030201 +2.0,12.89655172413793,1.222087,0.988575,0.028960 +2.0,13.172413793103448,1.221071,0.988899,0.027079 +2.0,13.448275862068964,1.220118,0.989208,0.025756 +2.0,13.724137931034482,1.219223,0.989204,0.034378 +2.0,14.0,1.218379,0.991537,0.022205 From f59d54ab12e01a527078067aed639a1636b23290 Mon Sep 17 00:00:00 2001 From: Nikhil Sarin Date: Wed, 9 Mar 2022 18:20:15 +0100 Subject: [PATCH 28/28] A bunch of new models --- redback/transient_models/supernova_models.py | 143 ++++++++++++++++--- 1 file changed, 122 insertions(+), 21 deletions(-) diff --git a/redback/transient_models/supernova_models.py b/redback/transient_models/supernova_models.py index 10cf5d9c..fd48c1f7 100644 --- a/redback/transient_models/supernova_models.py +++ b/redback/transient_models/supernova_models.py @@ -1,5 +1,6 @@ import numpy as np from redback.transient_models.phenomenological_models import exponential_powerlaw +from redback.transient_models.magnetar_models import _magnetar_only import redback.interaction_processes as ip import redback.sed as sed import redback.photosphere as photosphere @@ -67,7 +68,6 @@ def sn_exponential_powerlaw(time, redshift, lbol_0, alpha_1, alpha_2, tpeak_d, :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] - # time = time * 86400 frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -140,7 +140,6 @@ def arnett(time, redshift, f_nickel, mej, interaction_process=ip.Diffusion, :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] - # time = time * 86400 frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value @@ -185,7 +184,7 @@ def basic_magnetar_powered_bolometric(time, p0, bp, mass_ns, theta_pb,interactio e.g., for Diffusion: kappa, kappa_gamma, vej (km/s), temperature_floor :return: bolometric_luminosity """ - lbol = _basic_magnetar(time=time*86400, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) + lbol = _basic_magnetar(time=time*day_to_s, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) if interaction_process is not None: interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) lbol = interaction_class.new_luminosity @@ -267,6 +266,7 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu :param sed: Default is CutoffBlackbody. :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + and CutoffBlackbody: cutoff_wavelength, default is 3000 Angstrom :return: flux_density or magnitude depending on output_format kwarg """ frequency = kwargs['frequency'] @@ -276,8 +276,9 @@ def slsn(time, redshift, p0, bp, mass_ns, theta_pb, interaction_process=ip.Diffu lbol = slsn_bolometric(time=time, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb, interaction_process=interaction_process) photo = photosphere(time=time, luminosity=lbol, **kwargs) - sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, - frequency=frequency, luminosity_distance=dl, cutoff_wavelength=cutoff_wavelength) + sed_1 = sed(time=time, luminosity=lbol, temperature=photo.photosphere_temperature, + r_photosphere=photo.r_photosphere,frequency=frequency, luminosity_distance=dl, + cutoff_wavelength=cutoff_wavelength) flux_density = sed_1.flux_density @@ -311,7 +312,7 @@ def magnetar_nickel(time, redshift, f_nickel, mej, p0, bp, mass_ns, theta_pb, in frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) dl = cosmo.luminosity_distance(redshift).cgs.value - lbol_mag = _basic_magnetar(time=time*86400, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) + lbol_mag = _basic_magnetar(time=time*day_to_s, p0=p0, bp=bp, mass_ns=mass_ns, theta_pb=theta_pb) lbol_arnett = _nickelcobalt_engine(time=time, f_nickel=f_nickel, mej=mej) lbol = lbol_mag + lbol_arnett @@ -688,28 +689,128 @@ def csm_nickel(time, redshift, mej, f_nickel, csm_mass, ek, eta, rho, kappa, r0, return flux_density.to(uu.ABmag).value -@citation_wrapper('redback') -def type_1a_bolometric(): - pass - @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def type_1a(): - pass +def type_1a(time, redshift, f_nickel, mej, **kwargs): + """ + A nickel powered explosion with line absorption and cutoff blackbody SED for SNe 1A. + + :param time: time in days in observer frame + :param redshift: source redshift + :param f_nickel: fraction of nickel mass + :param mej: ejecta mass in solar masses + :param kwargs: kappa, kappa_gamma, vej (km/s), + temperature_floor (K), Cutoff_wavelength (default is 3000 Angstrom) + :return: flux_density or magnitude depending on output_format kwarg + """ + frequency = kwargs['frequency'] + cutoff_wavelength = kwargs.get('cutoff_wavelength', 3000) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + lbol = arnett_bolometric(time=time, f_nickel=f_nickel, mej=mej, + interaction_process=ip.Diffusion, **kwargs) + + photo = photosphere.TemperatureFloor(time=time, luminosity=lbol, **kwargs) + sed_1 = sed.CutoffBlackbody(time=time, luminosity=lbol, temperature=photo.photosphere_temperature, + r_photosphere=photo.r_photosphere,frequency=frequency, luminosity_distance=dl, + cutoff_wavelength=cutoff_wavelength) + sed_2 = sed.Line(time=time, luminosity=lbol, frequency=frequency, luminosity_distance=dl, + sed=sed_1, **kwargs) + + flux_density = sed_2.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value -@citation_wrapper('redback') -def type_1c_bolometric(): - pass @citation_wrapper('https://ui.adsabs.harvard.edu/abs/2018ApJS..236....6G/abstract') -def type_1c(): - pass +def type_1c(time, redshift, f_nickel, mej, **kwargs): + """ + A nickel powered explosion with synchrotron and blackbody SED's for SNe 1C. + + :param time: time in days in observer frame + :param redshift: source redshift + :param f_nickel: fraction of nickel mass + :param mej: ejecta mass in solar masses + :param kwargs: kappa, kappa_gamma, vej (km/s), temperature_floor (K), pp, nu_max (default is 1e9 Hz) + source_radius (default is 1e13), f0: synchrotron normalisation (default is 1e-26). + :return: flux_density or magnitude depending on output_format kwarg + """ + frequency = kwargs['frequency'] + pp = kwargs['pp'] + nu_max = kwargs.get('nu_max', 1e9) + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + lbol = arnett_bolometric(time=time, f_nickel=f_nickel, mej=mej, + interaction_process=ip.Diffusion, **kwargs) + + photo = photosphere.TemperatureFloor(time=time, luminosity=lbol, **kwargs) + sed_1 = sed.Blackbody(temperature=photo.photosphere_temperature, + r_photosphere=photo.r_photosphere,frequency=frequency, luminosity_distance=dl) + sed_2 = sed.Synchrotron(frequency=frequency, luminosity_distance=dl, pp=pp, nu_max=nu_max, **kwargs) + + flux_density = sed_1.flux_density + sed_2.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value @citation_wrapper('redback') -def general_magnetar_slsn_bolometric(): - pass +def general_magnetar_slsn_bolometric(time, l0, tsd, nn, interaction_process=ip.Diffusion, **kwargs): + """ + :param time: time in days in source frame + :param l0: magnetar energy normalisation in ergs + :param tsd: magnetar spin down damping timescale in source frame days + :param nn: braking index + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + :param kwargs: Must be all the kwargs required by the specific interaction_process, + e.g., for Diffusion: kappa, kappa_gamma, vej (km/s) + :return: + """ + lbol = _magnetar_only(time=time*day_to_s, l0=l0, tsd=tsd*day_to_s, nn=nn) + if interaction_process is not None: + interaction_class = interaction_process(time=time, luminosity=lbol, **kwargs) + lbol = interaction_class.new_luminosity + return lbol @citation_wrapper('redback') -def general_magnetar_slsn(): - pass +def general_magnetar_slsn(time, redshift, l0, tsd, nn, interaction_process = ip.Diffusion, + photosphere = photosphere.TemperatureFloor, sed = sed.Blackbody, ** kwargs): + """ + :param time: time in days in observer frame + :param redshift: source redshift + :param l0: magnetar energy normalisation in ergs + :param tsd: magnetar spin down damping timescale in source frame in days + :param nn: braking index + :param interaction_process: Default is Diffusion. + Can also be None in which case the output is just the raw engine luminosity, or another interaction process. + :param photosphere: Default is TemperatureFloor. + kwargs must have vej or relevant parameters if using different photosphere model + :param sed: Default is blackbody. + :param kwargs: Must be all the kwargs required by the specific interaction_process, photosphere, sed methods used + e.g., for Diffusion and TemperatureFloor: kappa, kappa_gamma, vej (km/s), temperature_floor + :return: flux_density or magnitude depending on output_format kwarg + """ + frequency = kwargs['frequency'] + frequency, time = calc_kcorrected_properties(frequency=frequency, redshift=redshift, time=time) + dl = cosmo.luminosity_distance(redshift).cgs.value + + lbol = general_magnetar_slsn_bolometric(time=time, l0=l0, tsd=tsd, nn=nn, + interaction_process = interaction_process, ** kwargs) + photo = photosphere(time=time, luminosity=lbol, **kwargs) + + sed_1 = sed(temperature=photo.photosphere_temperature, r_photosphere=photo.r_photosphere, + frequency = frequency, luminosity_distance = dl) + + flux_density = sed_1.flux_density + + if kwargs['output_format'] == 'flux_density': + return flux_density.to(uu.mJy).value + elif kwargs['output_format'] == 'magnitude': + return flux_density.to(uu.ABmag).value +