diff --git a/.gitignore b/.gitignore index cb6bb4a..33fec8a 100644 --- a/.gitignore +++ b/.gitignore @@ -111,3 +111,9 @@ venv.bak/ /docs/notebooks/*.rsp /docs/notebooks/*.fits /docs/notebooks/*.ipynb + +# hdf5 +*.h5 + +# Dask worker space +dask-worker-space/ \ No newline at end of file diff --git a/cosmogrb/grb/grb.py b/cosmogrb/grb/grb.py index 88ccc84..01a89fc 100644 --- a/cosmogrb/grb/grb.py +++ b/cosmogrb/grb/grb.py @@ -134,6 +134,26 @@ def display_energy_dependent_light_curve( ].display_energy_dependent_light_curve( time=time, energy=energy, ax=ax, cmap=cmap, **kwargs ) + + def display_time_dependent_spectrum( + self, time, energy, ax=None, cmap="viridis", **kwargs + ): + """FIXME! briefly describe function + + :param time: + :param energy: + :param ax: + :param cmap: + :returns: + :rtype: + + """ + + list(self._lightcurves.values())[ + 0 + ].display_time_dependent_spectrum( + time=time, energy=energy, ax=ax, cmap=cmap, **kwargs + ) def display_energy_integrated_light_curve(self, time, ax=None, **kwargs): """FIXME! briefly describe function diff --git a/cosmogrb/io/grb_save.py b/cosmogrb/io/grb_save.py index 798b01e..a90a528 100644 --- a/cosmogrb/io/grb_save.py +++ b/cosmogrb/io/grb_save.py @@ -18,6 +18,7 @@ def __init__( dec, duration, z, + time_adjustment, lightcurves, responses, source_params, @@ -72,7 +73,7 @@ def __init__( lightcurve=lightcurves[key], response=responses[key] ) - # + super(GRBSave, self).__init__(data) @@ -91,6 +92,7 @@ def __init__( self._dec = dec self._z = z self._duration = duration + self._time_adjustment = time_adjustment def __setitem__(self, key, value): @@ -135,6 +137,10 @@ def duration(self): def T0(self): return self._T0 + @property + def time_adjustment(self): + return self._time_adjustment + @property def extra_info(self): return self._extra_info @@ -247,6 +253,7 @@ def from_file(cls, file_name): ra=ra, dec=dec, z=z, + time_adjustment = time_adjustment, duration=duration, lightcurves=lightcurves, responses=responses, @@ -256,7 +263,7 @@ def from_file(cls, file_name): def __repr__(self): - return self._output().to_string() + return self._output().to_string(float_format=lambda x: "{:.3e}".format(x)) def info(self): diff --git a/cosmogrb/lightcurve/light_curve_storage.py b/cosmogrb/lightcurve/light_curve_storage.py index 1187909..f4ad070 100644 --- a/cosmogrb/lightcurve/light_curve_storage.py +++ b/cosmogrb/lightcurve/light_curve_storage.py @@ -473,7 +473,7 @@ def _bin_spectrum(self, tmin, tmax, times, pha): """ - idx = np.ones_like(times, dtype=bool) + idx = np.zeros_like(times, dtype=bool) # filter times diff --git a/cosmogrb/lightcurve/lightcurve.py b/cosmogrb/lightcurve/lightcurve.py index 6680de7..0361ead 100644 --- a/cosmogrb/lightcurve/lightcurve.py +++ b/cosmogrb/lightcurve/lightcurve.py @@ -220,6 +220,24 @@ def display_energy_dependent_light_curve( time=time, energy=energy, ax=ax, cmap=cmap, **kwargs ) + def display_time_dependent_spectrum( + self, time, energy, ax=None, cmap="viridis", **kwargs + ): + """FIXME! briefly describe function + + :param time: + :param energy: + :param ax: + :param cmap: + :returns: + :rtype: + + """ + + self._source.display_time_dependent_spectrum( + time=time, energy=energy, ax=ax, cmap=cmap, **kwargs + ) + def display_energy_integrated_light_curve(self, time, ax=None, **kwargs): """FIXME! briefly describe function diff --git a/cosmogrb/response/response.py b/cosmogrb/response/response.py index cc7cc51..35dbcf2 100644 --- a/cosmogrb/response/response.py +++ b/cosmogrb/response/response.py @@ -260,7 +260,7 @@ def to_fits( instrument_name, ) - fits_file.writeto(filename, clobber=overwrite) + fits_file.writeto(filename, overwrite=overwrite) __all__ = ["Response"] diff --git a/cosmogrb/sampler/source.py b/cosmogrb/sampler/source.py index 0d042eb..fa0f507 100644 --- a/cosmogrb/sampler/source.py +++ b/cosmogrb/sampler/source.py @@ -77,6 +77,24 @@ def display_energy_dependent_light_curve( time, energy, ax, cmap, **kwargs ) + def display_time_dependent_spectrum( + self, time, energy, ax=None, cmap="viridis", **kwargs + ): + """FIXME! briefly describe function + + :param time: + :param energy: + :param ax: + :param cmap: + :returns: + :rtype: + + """ + + self._source_function.display_time_dependent_spectrum( + time, energy, ax, cmap, **kwargs + ) + def set_response(self, response): """ called if there is no response upon creation diff --git a/cosmogrb/sampler/source_function.py b/cosmogrb/sampler/source_function.py index 13814f0..7271119 100644 --- a/cosmogrb/sampler/source_function.py +++ b/cosmogrb/sampler/source_function.py @@ -1,6 +1,10 @@ import abc +import numpy as np + +import matplotlib as mpl import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable from cosmogrb.utils.array_to_cmap import array_to_cmap @@ -97,9 +101,11 @@ def sample_energy(self, times): pass - def display_energy_integrated_light_curve(self, time, ax=None, **kwargs): + def display_energy_integrated_light_curve(self, time,ax=None, **kwargs): """ plot the latent light curve integrated over energy + in the observer frame + set F_peak is given in burst frame -> peak is in observer frame at F/(1+z^2) :param time: :param ax: @@ -116,20 +122,30 @@ def display_energy_integrated_light_curve(self, time, ax=None, **kwargs): fig = ax.get_figure() - y = [self.energy_integrated_evolution(t) for t in time] + n_energies = 75 + + energy_grid = np.power( + 10.0, np.linspace(np.log10(self._emin*(1+self._z)), np.log10(self._emax*(1+self._z)), n_energies) + ) + + y = np.zeros(len(time)) + for i,t in enumerate(time): + # flux at index (time, energy) + energy_slice = self.evolution(energy_grid, t) + y[i] = np.trapz(energy_slice[0, :]*energy_grid, energy_grid) ax.plot(time, y, **kwargs) - ax.set_xlabel("time") - ax.set_ylabel("flux") + ax.set_xlabel("Time [s]") + ax.set_ylabel(r"Flux [keV/s/cm$^2$]") return fig def display_energy_dependent_light_curve( - self, time, energy, ax=None, cmap="viridis", **kwargs + self, time, energy, ax=None, cmap="viridis", uselog=True, **kwargs ): """ - plot the latent light curve integrated over energy + plot the latent light curve for different energies :param time: :param ax: @@ -146,20 +162,85 @@ def display_energy_dependent_light_curve( fig = ax.get_figure() - # index (time, flux) + # flux at index (time, energy) grid = self.evolution(energy, time) - _, colors = array_to_cmap(energy, cmap=cmap, use_log=True) + _, colors = array_to_cmap(energy, cmap=cmap, use_log=uselog) + + cmap = mpl.colormaps[cmap] + + if uselog: + norm = mpl.colors.LogNorm(vmin=min(energy), vmax=max(energy)) + else: + norm = mpl.colors.Normalize(vmin=min(energy), vmax=max(energy)) + + #Add an axis for colorbar on right + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + + fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), + cax=cax, label='Energy [keV]') for i, lc in enumerate(grid.T): ax.plot(time, lc, color=colors[i], **kwargs) - ax.set_xlabel("time") + ax.set_xlabel("Time [s]") ax.set_ylabel("flux") return fig + def display_time_dependent_spectrum( + self, time, energy, ax=None, cmap="viridis", uselog=False, **kwargs + ): + """ + plot the latent nu *F_nu spectrum at different times + + :param time: + :param ax: + :returns: + :rtype: + + """ + + if ax is None: + + fig, ax = plt.subplots() + + else: + + fig = ax.get_figure() + + # index (time, flux) + grid = self.evolution(energy, time) + + _, colors = array_to_cmap(time, cmap=cmap, use_log=uselog) + + cmap = mpl.colormaps[cmap] + + if uselog: + norm = mpl.colors.LogNorm(vmin=min(time), vmax=max(time)) + else: + mpl.colors.Normalize(vmin=min(time), vmax=max(time)) + + #Add an axis for colorbar on right + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + + fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), + cax=cax, label='Time [s]') + + for i, lc in enumerate(grid): + + ax.plot(energy, energy**2*lc, color=colors[i], **kwargs) + + ax.set_xlabel("Energy [keV]") + ax.set_ylabel(r"$\nu F_{\nu}$ [keV$^2$/s/cm$^2$/keV]") + ax.set_xscale('log') + ax.set_yscale('log') + + return fig + @property def index(self): return self._index diff --git a/docs/notebooks/grb.md b/docs/notebooks/grb.md index b6498f0..5e0d70f 100644 --- a/docs/notebooks/grb.md +++ b/docs/notebooks/grb.md @@ -5,10 +5,10 @@ jupyter: text_representation: extension: .md format_name: markdown - format_version: '1.2' - jupytext_version: 1.8.0 + format_version: '1.3' + jupytext_version: 1.14.1 kernelspec: - display_name: python3 + display_name: Python 3 (ipykernel) language: python name: python3 --- @@ -16,14 +16,15 @@ jupyter: ```python import matplotlib.pyplot as plt import numpy as np -%matplotlib inline -#from jupyterthemes import jtplot -#plt.style.use('mike') -#jtplot.style(context='talk', fscale=1, grid=False) - +#print(plt.style.available) import cosmogrb ``` +```python +%matplotlib widget +plt.style.use('seaborn-v0_8-pastel') +``` + # GRBs This section describes how to handle the low level simulation of GRBs. AS the code currently is built for simulating GRBs as observed by Fermi-GBM, we will focus our attention there. As the code expands, I will update the docs. @@ -45,7 +46,7 @@ grb = cosmogrb.gbm.GBMGRB_CPL( ra=312.0, dec=-62.0, z=1.0, - peak_flux=5e-7, + peak_flux=1e-6, alpha=-0.66, ep=500.0, tau=2.0, @@ -66,8 +67,6 @@ grb.info() time = np.linspace(0, 20, 500) grb.display_energy_integrated_light_curve(time, color="#A363DE"); - - ``` ```python hidden=true @@ -76,6 +75,14 @@ energy = np.logspace(1, 3, 1000) grb.display_energy_dependent_light_curve(time, energy, cmap='PRGn', lw=.25, alpha=.5) ``` +```python +fig,ax = plt.subplots(figsize=(7,5)) +time = np.linspace(0, 5, 15) +energy = np.logspace(1, 3, 100) +grb.display_time_dependent_spectrum(time,energy,ax=ax) +ax.set_ylim(1e-3,9e2) +``` + ## Simulate the GRB Now we can create all the light curves from the GRB. Since are not currently running a Dask server, we tell the GRB to process serially, i.e., computing each light curve one at a time. @@ -139,11 +146,13 @@ for k,v in grb_reload.items(): axes[3,2].set_visible(False) axes[3,3].set_visible(False) plt.tight_layout() +plt.show() ``` And we can look at the generated count spectra/ ```python +%matplotlib widget fig, axes = plt.subplots(4,4,sharex=False,sharey=False,figsize=(10,10)) row=0 col = 0 @@ -153,9 +162,9 @@ for k, v in grb_reload.items(): lightcurve = v['lightcurve'] - lightcurve.display_count_spectrum(tmin=0, tmax=5, ax=ax,color='#25C68C') - lightcurve.display_count_spectrum_source(tmin=0, tmax=5, ax=ax,color="#A363DE") - lightcurve.display_count_spectrum_background(tmin=0, tmax=5, ax=ax, color="#2C342E") + lightcurve.display_count_spectrum(tmin=0, tmax=5,ax=ax,color='#25C68C',label='All') + lightcurve.display_count_spectrum_source(tmin=0, tmax=5, ax=ax,color="#A363DE",label='Source') + lightcurve.display_count_spectrum_background(tmin=0, tmax=5, ax=ax, color="#2C342E",label='Bkg') ax.set_title(k,size=8) if col < 3: @@ -163,7 +172,8 @@ for k, v in grb_reload.items(): else: row+=1 col=0 - + +axes[0,0].legend() axes[3,2].set_visible(False) axes[3,3].set_visible(False) plt.tight_layout() @@ -175,7 +185,7 @@ In the case of GBM, we can convert the saved HDF5 files into TTE files for analy ```python cosmogrb.grbsave_to_gbm_fits("test_grb.h5") -!ls SynthGRB_* +#!ls SynthGRB_* ``` ```python