You're Invited:Meet the Socket Team at RSAC and BSidesSF 2026, March 23–26.RSVP →
Socket
Book a DemoSign in
Socket
Book a DemoSign in

afterglowpy

Package Overview
Dependencies
Maintainers
1
Versions
15
Alerts
File Explorer

Advanced tools

License
Socket logo

Install Socket

Detect and block malicious and high-risk dependencies

Install

afterglowpy - pypi Package Compare versions

Comparing version
0.7.3
 to
0.8.0
+21
LICENSE.txt
MIT License
Copyright (c) 2018 Geoffrey Ryan
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
+95
-104
afterglowpy.egg-info/PKG-INFO
Metadata-Version: 2.1
Name: afterglowpy
Version: 0.7.3
Version: 0.8.0
Summary: GRB Afterglow Models
Home-page: https://github.com/geoffryan/afterglowpy
Author: Geoffrey Ryan
Author-email: gsryan@umd.edu
Author-email: gryan@perimeterinstitute.ca
License: MIT
Project-URL: Source Code, https://github.com/geoffryan/afterglowpy
Project-URL: Documentation, https://afterglowpy.readthedocs.io
Description: # Numeric GRB Afterglow models
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2019](https://arxiv.org/abs/1909.11691). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.
Documentation available at <https://afterglowpy.readthedocs.io/>
## Attribution
If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2019, arXiv:1910.11691" [arXiv link](https://arxiv.org/abs/1909.11691).
## Features
_afterglowpy_ computes synchrotron emission from the forward shock of a relativistic blast wave. It includes:
- Fully trans-relativistic shock evolution through a constant density medium.
- On-the-fly integration over the equal-observer-time slices of the shock surface.
- Approximate prescription for jet spreading.
- Arbitrary viewing angles.
- Angularly structured jets, ie. E(&theta;)
- Spherical velocity-stratified outflows, ie. E(u)
- Counter-jet emission.
It has limited support (these should be considered experimental) for:
- Initial energy injection
- Inverse comption spectra
- Spatially-resolved intensity maps
- Early coasting phase
It does not include (yet):
- External wind medium, ie. n &prop; r<sup>-2</sup>
- Synchrotron self-absorbtion
- Reverse shock emission
_afterglowpy_ has been calibrated to the BoxFit code ([van Eerten, van der Horst, & Macfadyen 2011](https://arxiv.org/abs/1110.5089), available at the [Afterglow Library](https://cosmo.nyu.edu/afterglowlibrary/boxfit2011.html)) and produces similar light curves for top hat jets (within 50% when same parameters are used) both on- and off-axis. Its jet models by default do not include an initial coasting phase, which may effect predictions for early observations.
## Installation/Building
_afterglowpy_ is available via `pip`:
```bash
$ pip install afterglowpy
```
If you are working on a local copy of this repo and would like to install from source, you can the run the following from the top level directory of the project.
```bash
$ pip install -e .
```
## Using
**This interface will be updated to be more sensible in the VERY near future**
In your python code, import the library with `import afterglowpy as grb`.
The main function of interest is`grb.fluxDensity(t, nu, jetType, specType, *pars, **kwargs)`. See `tests/plotLC.py` for a simple example.
`jetType` can be -1 (top hat), 0 (Gaussian), 1 (Power Law w/ core), 2 (Gaussian w/ core), 3 (Cocoon), or 4 (Smooth Power Law).
`specType` can be 0 (global cooling time, no inverse compton) or 1 (global cooling time, inverse compton).
For jet-like afterglows (`jetTypes` -2, -1, 0, 1, 2, and 4) `pars` has 14 positional arguments:
- `0 thetaV` viewing angle in radians
- `1 E0` on-axis isotropic equivalent energy in erg
- `2 thetaC` half-width of the jet core in radians (jetType specific)
- `3 thetaW` "wing" truncation angle of the jet, in radians
- `4 b` power for power-law structure, &theta;<sup>-b</sup>
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
For cocoon-like afterglows (`jetType` 3) `pars` has 14 positional arguments:
- `0 umax` Initial maximum outflow 4-velocity
- `1 umin` Minium outflow 4-velocity
- `2 Ei` Fiducial energy in velocity distribution, E(>u) = E<sub>i</sub> u<sup>-k</sup>.
- `3 k` Power-law index of energy velocity distribution
- `4 Mej` Mass of material at `umax' in solar masses
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
Keyword arguments are:
- `z` redshift (defaults to 0)
- `tRes` time resolution of shock-evolution scheme, number of sample points per decade in time
- `latRes` latitudinal resolution for structured jets, number of shells per `thetaC`
- `rtol` target relative tolerance of flux integration
- `spread` boolean (defaults to True), whether to allow the jet to spread.
Platform: UNKNOWN
Classifier: Programming Language :: Python :: 3

@@ -121,2 +19,95 @@ Classifier: Programming Language :: C

Description-Content-Type: text/markdown
License-File: LICENSE.txt
Requires-Dist: numpy>=1.13
Requires-Dist: scipy>=0.14
Provides-Extra: docs
Requires-Dist: numpydoc; extra == "docs"
# Numeric GRB Afterglow models
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2020](https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.
Documentation available at <https://afterglowpy.readthedocs.io/>
## Attribution
If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2020, Astrophysical Journal *896*, 166 (2020)" [arXiv link](https://arxiv.org/abs/1909.11691).
## Acknowledgements
This work is funded in part by the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158).
## Features
_afterglowpy_ computes synchrotron emission from the forward shock of a relativistic blast wave. It includes:
- Fully trans-relativistic shock evolution through a constant density medium.
- On-the-fly integration over the equal-observer-time slices of the shock surface.
- Approximate prescription for jet spreading.
- Arbitrary viewing angles.
- Angularly structured jets, ie. E(&theta;)
- Spherical velocity-stratified outflows, ie. E(u)
- Counter-jet emission.
- Deep Newtonian emission.
- Image moments suitable for astrometry: centroid position and image size.
It has limited support (these should be considered experimental) for:
- Initial energy injection
- Inverse comption spectra
- Early coasting phase
It does not include (yet):
- External wind medium, ie. n &prop; r<sup>-2</sup>
- Synchrotron self-absorbtion
- Reverse shock emission
_afterglowpy_ has been calibrated to the BoxFit code ([van Eerten, van der Horst, & Macfadyen 2011](https://arxiv.org/abs/1110.5089), available at the [Afterglow Library](https://cosmo.nyu.edu/afterglowlibrary/boxfit2011.html)) and produces similar light curves for top hat jets (within 50% when same parameters are used) both on- and off-axis. Its jet models by default do not include an initial coasting phase, which may effect predictions for early observations.
## Installation/Building
_afterglowpy_ is available via `pip`:
```bash
$ pip install afterglowpy
```
If you are working on a local copy of this repo and would like to install from source, you can the run the following from the top level directory of the project.
```bash
$ pip install -e .
```
## Using
In your python code, import the library with `import afterglowpy as grb`.
The main function of interest is`grb.fluxDensity(t, nu, **kwargs)`. See `examples/plotLightCurve.py` for a simple example.
For jet-like afterglows there are up to 13 required keyword arguments:
- `jetType` an integer code setting the jet structure. It can be `grb.jet.TopHat`, `grb.jet.Gaussian`, `grb.jet.PowerLawCore`, `grb.jet.GaussianCore`, `grb.jet.Spherical`, or `grb.jet.PowerLaw`.
- `specType` an integer code specifying flags for the emissivity function and spectrum. Can be `grb.jet.SimpleSpec` (basic spectrum with &nu;<sub>m</sub> and &nu;<sub>c</sub>), `grb.jet.DeepNewtonian`, `grb.jet.ICCooling` (simple inverse Compton effects on the cooling frequency, experimental).
- `thetaObs` viewing angle in radians
- `E0` on-axis isotropic equivalent energy in erg
- `thetaCore` half-width of the jet core in radians (jetType specific)
- `thetaWing` "wing" truncation angle of the jet, in radians
- `b` power for power-law structure, &theta;<sup>-b</sup>
- `n0` Number density of ISM, in cm<sup>-3</sup>
- `p` Electron distribution power-law index (p>2)
- `epsilon_e` Thermal energy fraction in electrons
- `epsilon_B` Thermal energy fraction in magnetic field
- `xi_N` Fraction of electrons that get accelerated
- `d_L` Luminosity distance in cm
Optional keyword arguments for all models are:
- `z` redshift (defaults to 0)
- `spread` boolean (defaults to True), whether to allow the jet to spread.
- `counterjet` boolean (defaults to False), whether to include the counterjet
- `moment` array (integer dtype, same shape as t and nu) which sky moment to compute.
- `L0` Fiducial luminosity for energy injection, in erg/s, default 0.0.
- `q` Temporal power-law index for energy injection, default 0.0.
- `ts` Fiducial time-scale for energy injection, in seconds, default 0.
- `tRes` time resolution of shock-evolution scheme, number of sample points per decade in time
- `latRes` latitudinal resolution for structured jets, number of shells per `thetaC`
- `rtol` target relative tolerance of flux integration
+1
-1
afterglowpy.egg-info/requires.txt

@@ -1,2 +0,2 @@

numpy>=1.10
numpy>=1.13
scipy>=0.14

@@ -3,0 +3,0 @@ 

+1
-0
afterglowpy.egg-info/SOURCES.txt

@@ -0,1 +1,2 @@

LICENSE.txt
MANIFEST.in

@@ -2,0 +3,0 @@ README.md

+13
-2
afterglowpy/cocoon.py

@@ -48,3 +48,6 @@ import math

em = jet.emissivity(nu, r, mu, te, u, us, n0, p, epsE,
rho0 = n0 * mp
Msw = rho0 * 4.0/3.0 * np.pi * r**3
em = jet.emissivity(nu, r, mu, te, u, us, rho0, Msw, p, epsE,
epsB, ksiN, specType)

@@ -82,2 +85,8 @@

# Environment Variables
envType = kwargs['envType'] if 'envType' in kwargs else jet.EnvISM
R0Env = kwargs['R0Env'] if 'R0Env' in kwargs else 1.0e18
kEnv = kwargs['kEnv'] if 'kEnv' in kwargs else 0.0
rho1Env = kwargs['rho1Env'] if 'rho1Env' in kwargs else 1.0
rho0 = mp * n0

@@ -104,3 +113,5 @@ Mej = MFast_solar * Msun

ar, au = shock.shockEvolRK4(ate, r0, uMax,
MFast_solar*Msun, rho0, Er, k, uMin, L0, q, ts)
MFast_solar*Msun,
envType, rho0, R0Env, kEnv, rho1Env,
Er, k, uMin, L0, q, ts)

@@ -107,0 +118,0 @@ P = np.zeros(t.shape)

+25
-16
afterglowpy/flux.py

@@ -49,4 +49,8 @@ from . import cocoon

specType : int
Code for type of spectrum. Options are: 0 broken power law
(Ryan+ 2020), 1 broken power law w/ inverse Compton cooling. Default: 0
Flags for type of spectrum/emissivity function. Spectrum flags are
available in ``afterglowpy.jet`` and include: ``jet.SimpleSpec``
broken power law with nu_m and nu_c (Ryan+ 2020, default),
``jet.ICCooling`` simple inverse-compton contribution to cooling,
``jet.DeepNewtonian`` better handling of late-time emission when
some electrons become non-relativistic (e.g. Sironi+ 2013).
thetaObs : float

@@ -114,2 +118,10 @@ Viewing angle in radians. Jet models only.

Whether to include counterjet emission. Defaults to False.
moment : array_like, optional
An integer array the same shape as the larger of `t` or `nu`. Selects
which image moment to compute. Observer's sky is on the x-y plane,
with the jet propagating in the x-direction. Moments are in terms of
proper length (in cm) Options are: `jet.MOM_0` Flux, default,
`jet.MOM_X` x^1 moment, `jet.MOM_Y` y^1 moment, `jet.MOM_Z` z^1 moment,
`jet.MOM_XX` x^2 moment, `jet.MOM_YY` y^2 moment, `jet.MOM_ZZ`, z^2
moment, `jet.MOM_XY` x^1 y^1, `jet.MOM_YZ`, `jet.MOM_XZ`.
tRes : int, optional

@@ -155,3 +167,5 @@ Time resolution, number of points per decade in t, for shock evolution.

# Check Arguments, will raise ValueError if args are bad
t, nu = checkTNu(t, nu)

@@ -165,2 +179,3 @@ jetType = argsDict['jetType']

# arguments are good, full steam ahead!

@@ -172,2 +187,3 @@ z = argsDict.pop('z') if 'z' in argsDict else 0.0

# Default spreading method

@@ -181,3 +197,2 @@ if 'spread' in argsDict:

# Intercept background luminosities
# This was a bad idea to add to this function, but is kept for

@@ -194,3 +209,3 @@ # backwards compatibility. Please don't use these.

if argsDict['jetType'] == jet.Spherical:
if jetType == jet.Spherical:
Fnu.flat[:] = cocoon.fluxDensity(tz.flat, nuz.flat, **argsDict)

@@ -223,3 +238,3 @@ else:

Fnu *= 1+z
return Fnu

@@ -489,7 +504,4 @@

for _, x in argsDict.items():
if not np.isfinite(x):
raise ValueError("All parameters must be finite")
return
jetType = argsDict['jetType']

@@ -681,6 +693,8 @@ specType = argsDict['specType']

'thetaWing', 'b', 'L0', 'q', 'ts', 'n0', 'p', 'epsilon_e',
'epsilon_B', 'xi_N', 'd_L', 'g0', 'LR', 'LO', 'LX', 'tAdd', 'z']
'epsilon_B', 'xi_N', 'd_L', 'g0', 'LR', 'LO', 'LX', 'tAdd', 'z',
'envType', 'R0Env', 'kEnv', 'rho1Env']
sphKeys = ['jetType', 'specType', 'uMax', 'uMin', 'Er',
'k', 'MFast_solar', 'L0', 'q', 'ts', 'n0', 'p', 'epsilon_e',
'epsilon_B', 'xi_N', 'd_L', 'g0', 'LR', 'LO', 'LX', 'tAdd', 'z']
'epsilon_B', 'xi_N', 'd_L', 'g0', 'LR', 'LO', 'LX', 'tAdd', 'z',
'envType', 'R0Env', 'kEnv', 'rho1Env']

@@ -695,6 +709,1 @@ jetType = args[0]

return argsDict
+34
-6
afterglowpy/integrate.c

@@ -350,8 +350,11 @@ #include <math.h>

double *eps, int verbose, int (*errf)(void *),
double *pfa, double *pfb)
double *pfa, double *pfb, Mesh9 *mout)
{
double I = m9_adapt(f, xa, xb, Nmax, cadreInitInterval,
cadreProcessInterval, cadreSplitInterval,
atol, rtol, args, Neval, eps, NULL, verbose, errf,
atol, rtol, args, Neval, eps, mout, verbose, errf,
pfa, pfb);
return I;

@@ -1468,3 +1471,6 @@ }

{
mesh9Free(&m);
if(mout == NULL)
mesh9Free(&m);
else
mout = &m;
return 0.0;

@@ -1476,3 +1482,6 @@ }

{
mesh9Free(&m);
if(mout == NULL)
mesh9Free(&m);
else
mout = &m;
return 0.0;

@@ -1514,3 +1523,6 @@ }

{
mesh9Free(&m);
if(mout == NULL)
mesh9Free(&m);
else
mout = &m;
return 0.0;

@@ -1533,3 +1545,6 @@ }

{
mesh9Free(&m);
if(mout == NULL)
mesh9Free(&m);
else
mout = &m;
return 0.0;

@@ -1576,2 +1591,15 @@ }

// DEBUG DELETE LATER
/*
char *meshStr;
mesh9Write(&m, &meshStr);
FILE *fp = fopen("afterglowpy_int_dump.txt", "w");
fprintf(fp, "%s\n", meshStr);
int idx;
for(idx=0; idx<m.N; idx++)
interval9Write(m.heap+idx, fp);
fclose(fp);
free(meshStr);
*/
if(mout == NULL)

@@ -1578,0 +1606,0 @@ mesh9Free(&m);

+1
-1
afterglowpy/integrate.h

@@ -50,3 +50,3 @@ #ifndef AFTERGLOWPY_INTEGRATE

double *eps, int verbose, int (*errf)(void *),
double *pfa, double *pfb);
double *pfa, double *pfb, Mesh9 *mout);
double gk49_adapt(double (*f)(double, void *), double xa, double xb, int Nmax,

@@ -53,0 +53,0 @@ double atol, double rtol, void *args, int *Neval,

+7
-4
afterglowpy/interval.c

@@ -483,4 +483,7 @@ #include <stdio.h>

m->N = 0;
free(m->heap);
m->heap = NULL;
if(m->heap != NULL)
{
free(m->heap);
m->heap = NULL;
}
}

@@ -626,5 +629,5 @@

i->a, i->b, i->I, i->err, i->refinement);
fprintf(stream, " [%.3le %.3le %.3le %.3le %.3le %.3le"
" %.3le %.3le %.3le]\n", i->fa, i->fll, i->fl, i->flr,
fprintf(stream, " [%.6le %.6le %.6le %.6le %.6le %.6le"
" %.6le %.6le %.6le]\n", i->fa, i->fll, i->fl, i->flr,
i->fm, i->frl, i->fr, i->frr, i->fb);
}
+210
-58
afterglowpy/jetmodule.c
#include <Python.h>
#define NPY_NO_DEPRECATED_API NPY_1_11_API_VERSION
#define NPY_NO_DEPRECATED_API NPY_1_13_API_VERSION
#include <numpy/arrayobject.h>
#include <time.h>
#include "offaxis_struct.h"
#include "shockEvolution.h"

@@ -19,5 +20,6 @@ #define PROFILE

static const double b_default = 0.0;
static const double L0_default = 0.0;
static const double q_default = 0.0;
static const double ts_default = 0.0;
static const double L0_inj_default = 0.0;
static const double q_inj_default = 0.0;
static const double ts_inj_default = 0.0;
static const double t0_inj_default = 1.0e3;
static const double n_0_default = 1.0;

@@ -36,2 +38,8 @@ static const double p_default = 2.2;

static const double g0_default = -1.0;
static const int envType_default = ENV_ISM;
static const double R0_env_default = 1e18;
static const double k_env_default = 0.0;
static const double rho1_env_default = 1.0;
static const double E_core_global_default = 0.0;

@@ -163,2 +171,3 @@ static const double theta_h_core_global_default = 0.0;

PyModule_AddIntConstant(module, "Spherical", _spherical);
PyModule_AddIntConstant(module, "TrapFixed", INT_TRAP_FIXED);

@@ -175,2 +184,3 @@ PyModule_AddIntConstant(module, "TrapAdapt", INT_TRAP_ADAPT);

PyModule_AddIntConstant(module, "GK1021Adapt", INT_GK1021_ADAPT);
PyModule_AddIntConstant(module, "GammaInf", GAMMA_INF);

@@ -180,2 +190,4 @@ PyModule_AddIntConstant(module, "GammaFlat", GAMMA_FLAT);

PyModule_AddIntConstant(module, "GammaStruct", GAMMA_STRUCT);
PyModule_AddIntConstant(module, "SimpleSpec", SIMPLE_SPEC);
PyModule_AddIntConstant(module, "ICCooling", IC_COOLING_FLAG);

@@ -185,2 +197,19 @@ PyModule_AddIntConstant(module, "EpsEBar", EPS_E_BAR_FLAG);

PyModule_AddIntConstant(module, "SSASharp", SSA_SHARP_FLAG);
PyModule_AddIntConstant(module, "NoCooling", NO_COOLING_FLAG);
PyModule_AddIntConstant(module, "DeepNewtonian", DEEP_NEWTONIAN_FLAG);
PyModule_AddIntConstant(module, "BulkBM", BULK_BM_FLAG);
PyModule_AddIntConstant(module, "EnvISM", ENV_ISM);
PyModule_AddIntConstant(module, "EnvWind", ENV_WIND);
PyModule_AddIntConstant(module, "EnvPL", ENV_PL);
PyModule_AddIntConstant(module, "EnvStep", ENV_STEP);
PyModule_AddIntConstant(module, "MOM_0", MOM_0);
PyModule_AddIntConstant(module, "MOM_X", MOM_X);
PyModule_AddIntConstant(module, "MOM_Y", MOM_Y);
PyModule_AddIntConstant(module, "MOM_Z", MOM_Z);
PyModule_AddIntConstant(module, "MOM_XX", MOM_XX);
PyModule_AddIntConstant(module, "MOM_YY", MOM_YY);
PyModule_AddIntConstant(module, "MOM_ZZ", MOM_ZZ);
PyModule_AddIntConstant(module, "MOM_XY", MOM_XY);
PyModule_AddIntConstant(module, "MOM_YZ", MOM_YZ);
PyModule_AddIntConstant(module, "MOM_XZ", MOM_XZ);

@@ -205,2 +234,3 @@ #if PY_MAJOR_VERSION >= 3

PyObject *mask_obj = NULL;
PyObject *moment_obj = NULL;

@@ -223,5 +253,6 @@ #ifdef PROFILE

double b = b_default;
double L0 = L0_default;
double q = q_default;
double ts = ts_default;
double L0_inj = L0_inj_default;
double q_inj = q_inj_default;
double t0_inj = t0_inj_default;
double ts_inj = ts_inj_default;
double n_0 = n_0_default;

@@ -237,2 +268,8 @@ double p = p_default;

double g0 = g0_default;
int envType = envType_default;
double R0_env = R0_env_default;
double k_env = k_env_default;
double rho1_env = rho1_env_default;
double E_core_global = E_core_global_default;

@@ -254,5 +291,9 @@ double theta_h_core_global = theta_h_core_global_default;

"thetaObs", "E0", "thetaCore", "thetaWing",
"b", "L0", "q", "ts", "n0", "p",
"b",
"L0", "q", "ts",
"n0", "p",
"epsilon_e", "epsilon_B", "xi_N", "d_L",
"g0",
"envType", "R0Env", "kEnv", "rho1Env",
"t0_inj",
"E0Global", "thetaCoreGlobal",

@@ -262,14 +303,18 @@ "tRes", "latRes", "intType", "rtolStruct",

"mask",
"spread", "counterjet", "gammaType",
"spread", "counterjet", "gammaType", "moment",
NULL};
//printf("About to parse args\n");
//Parse Arguments
if(!PyArg_ParseTupleAndKeywords(args, kwargs,
"OO|ii""ddddddddddddddd""dd""iiidddii""O""iii",
//"OO|ii ddddddddddddddd dd iiidddii O iii",
"OO|ii""ddddddddddddddd""iddd""d""dd""iiidddii""O""iii""O",
kwlist,
&t_obj, &nu_obj, &jet_type, &spec_type,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b, &L0,
&q, &ts, &n_0, &p, &epsilon_E, &epsilon_B, &ksi_N, &d_L,
&t_obj, &nu_obj,
&jet_type, &spec_type,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b,
&L0_inj, &q_inj, &ts_inj,
&n_0, &p, &epsilon_E, &epsilon_B, &ksi_N, &d_L,
&g0,
&envType, &R0_env, &k_env, &rho1_env,
&t0_inj,
&E_core_global, &theta_h_core_global,

@@ -279,8 +324,10 @@ &tRes, &latRes, &int_type, &rtol_struct, &rtol_phi,

&mask_obj,
&spread, &counterjet, &gamma_type))
&spread, &counterjet, &gamma_type, &moment_obj))
{
//PyErr_SetString(PyExc_RuntimeError, "Could not parse arguments.");
PyErr_SetString(PyExc_RuntimeError, "Could not parse arguments.");
return NULL;
}
//printf("Args loaded\n");
if(int_type < 0 || int_type >= INT_UNDEFINED)

@@ -293,2 +340,4 @@ {

//printf("Args parsed\n");
//Grab NUMPY arrays

@@ -298,3 +347,7 @@ PyArrayObject *t_arr;

PyArrayObject *mask_arr = NULL;
PyArrayObject *moment_arr = NULL;
int givenMask = mask_obj != NULL && mask_obj != Py_None;
int givenMoment = moment_obj != NULL && moment_obj != Py_None;
t_arr = (PyArrayObject *) PyArray_FROM_OTF(t_obj, NPY_DOUBLE,

@@ -304,8 +357,11 @@ NPY_ARRAY_IN_ARRAY);

NPY_ARRAY_IN_ARRAY);
if(mask_obj != NULL)
if(givenMask)
mask_arr = (PyArrayObject *) PyArray_FROM_OTF(mask_obj, NPY_DOUBLE,
NPY_ARRAY_IN_ARRAY);
if(givenMoment)
moment_arr = (PyArrayObject *) PyArray_FROM_OTF(moment_obj, NPY_INT64,
NPY_ARRAY_IN_ARRAY);
if(t_arr == NULL || nu_arr == NULL || (mask_obj != NULL
&& mask_arr == NULL))
if(t_arr == NULL || nu_arr == NULL || (givenMask && mask_arr == NULL)
|| (givenMoment && moment_arr == NULL))
{

@@ -316,4 +372,7 @@ PyErr_SetString(PyExc_RuntimeError, "Could not read input arrays.");

Py_XDECREF(mask_arr);
Py_XDECREF(moment_arr);
return NULL;
}
//printf("Arrays grabbed\n");

@@ -323,6 +382,10 @@ int t_ndim = (int) PyArray_NDIM(t_arr);

int mask_ndim = 0;
int moment_ndim = 0;
if(mask_obj != NULL)
mask_ndim = (int) PyArray_NDIM(mask_arr);
if(givenMoment)
moment_ndim = (int) PyArray_NDIM(moment_arr);
if(t_ndim != 1 || nu_ndim != 1 || (mask_obj != NULL && mask_ndim != 1))
if(t_ndim != 1 || nu_ndim != 1 || (givenMask && mask_ndim != 1)
|| (givenMoment && moment_ndim != 1))
{

@@ -332,6 +395,10 @@ PyErr_SetString(PyExc_RuntimeError, "Arrays must be 1-D.");

Py_DECREF(nu_arr);
if(mask_obj != NULL)
if(mask_arr != NULL)
Py_DECREF(mask_arr);
if(moment_arr != NULL)
Py_DECREF(moment_arr);
return NULL;
}
//printf("NDims checked\n");

@@ -341,6 +408,9 @@ int N = (int)PyArray_DIM(t_arr, 0);

int Nmask = 0;
if(mask_obj != NULL)
if(givenMask)
Nmask = (int)PyArray_DIM(mask_arr, 0);
int Nmoment = 0;
if(givenMoment)
Nmoment = (int)PyArray_DIM(moment_arr, 0);
if(N != Nnu)
if(N != Nnu || (givenMoment && (Nmoment != N)))
{

@@ -350,4 +420,6 @@ PyErr_SetString(PyExc_RuntimeError, "Arrays must be same size.");

Py_DECREF(nu_arr);
if(mask_obj != NULL)
if(mask_arr != NULL)
Py_DECREF(mask_arr);
if(moment_arr != NULL)
Py_DECREF(moment_arr);
return NULL;

@@ -362,4 +434,8 @@ }

Py_DECREF(mask_arr);
if(moment_arr != NULL)
Py_DECREF(moment_arr);
return NULL;
}
//printf("Dims checked\n");

@@ -369,6 +445,12 @@ double *t = (double *)PyArray_DATA(t_arr);

double *mask = NULL;
if(mask_obj != NULL)
if(mask_arr != NULL)
mask = (double *)PyArray_DATA(mask_arr);
int masklen = Nmask/9;
long *moment = NULL;
if(moment_arr != NULL)
moment = (long *)PyArray_DATA(moment_arr);
//printf("Data got\n");
//Allocate output array

@@ -384,2 +466,6 @@

Py_DECREF(nu_arr);
if(mask_arr != NULL)
Py_DECREF(mask_arr);
if(moment_arr != NULL)
Py_DECREF(moment_arr);
return NULL;

@@ -410,4 +496,5 @@ }

theta_h_wing, b,
L0, q, ts,
L0_inj, q_inj, t0_inj, ts_inj,
n_0, p, epsilon_E, epsilon_B, ksi_N, g0,
envType, R0_env, k_env, rho1_env,
E_core_global, theta_h_core_global, ta, tb,

@@ -420,4 +507,8 @@ tRes, latRes, int_type,

//printf("Ready to go!\n");
// Calculate the flux!
calc_flux_density(jet_type, spec_type, t, nu, Fnu, N, &fp);
calc_flux_density(jet_type, spec_type, t, nu, Fnu, moment, N, &fp);
//printf("Gone!\n");

@@ -431,4 +522,6 @@ if(fp.error)

//Free the parameters!
free_fluxParams(&fp);
//printf("Freed params!\n");

@@ -446,2 +539,6 @@ #ifdef PROFILE2

Py_DECREF(mask_arr);
if(moment_obj != NULL)
Py_DECREF(moment_arr);
//printf("Decref'ed inputs!\n");

@@ -451,2 +548,4 @@ //Build output

//printf("Built outputs!\n");
#ifdef PROFILE1

@@ -474,8 +573,9 @@ //Profile 1 and output

int spec_type = 0;
double nu, R, mu, te, u, us, n0, p, epse, epsB, xi_N;
double nu, R, mu, te, u, us, rho0, Msw, p, epse, epsB, xi_N;
//Parse Arguments
if(!PyArg_ParseTuple(args, "ddddddddddd|i", &nu, &R, &mu, &te,
&u, &us, &n0, &p, &epse, &epsB, &xi_N, &spec_type))
if(!PyArg_ParseTuple(args, "dddddddddddd|i", &nu, &R, &mu, &te,
&u, &us, &rho0, &Msw, &p, &epse, &epsB, &xi_N,
&spec_type))
{

@@ -487,4 +587,4 @@ //PyErr_SetString(PyExc_RuntimeError, "Could not parse arguments.");

// Calculate it!
double em = emissivity(nu, R, mu, te, u, us, n0, p, epse, epsB,
xi_N, spec_type);
double em = emissivity(nu, R, mu, te, u, us, rho0, Msw, p, epse, epsB,
xi_N, spec_type);

@@ -512,5 +612,6 @@ //Build output

double b = b_default;
double L0 = L0_default;
double q = q_default;
double ts = ts_default;
double L0_inj = L0_inj_default;
double q_inj = q_inj_default;
double t0_inj = t0_inj_default;
double ts_inj = ts_inj_default;
double n_0 = n_0_default;

@@ -526,2 +627,8 @@ double p = p_default;

double g0 = g0_default;
int envType = envType_default;
double R0_env = R0_env_default;
double k_env = k_env_default;
double rho1_env = rho1_env_default;
double E_core_global = E_core_global_default;

@@ -546,2 +653,4 @@ double theta_h_core_global = theta_h_core_global_default;

"g0",
"envType", "R0Env", "kEnv", "rho1Env",
"t0_inj",
"E0Global", "thetaCoreGlobal",

@@ -556,8 +665,11 @@ "tRes", "latRes", "intType", "rtolStruct",

if(!PyArg_ParseTupleAndKeywords(args, kwargs,
"OOOO|ii""ddddddddddddddd""dd""iiidddii""O""iii",
"OOOO|ii""ddddddddddddddd""iddd""d""dd""iiidddii""O""iii",
kwlist,
&theta_obj, &phi_obj, &t_obj, &nu_obj, &jet_type, &spec_type,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b, &L0,
&q, &ts, &n_0, &p, &epsilon_E, &epsilon_B, &ksi_N, &d_L,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b,
&L0_inj, &q_inj, &ts_inj,
&n_0, &p, &epsilon_E, &epsilon_B, &ksi_N, &d_L,
&g0,
&envType, &R0_env, &k_env, &rho1_env,
&t0_inj,
&E_core_global, &theta_h_core_global,

@@ -704,4 +816,5 @@ &tRes, &latRes, &int_type, &rtol_struct, &rtol_phi,

theta_h_wing, b,
L0, q, ts,
L0_inj, q_inj, t0_inj, ts_inj,
n_0, p, epsilon_E, epsilon_B, ksi_N, g0,
envType, R0_env, k_env, rho1_env,
E_core_global, theta_h_core_global, ta, tb,

@@ -754,5 +867,6 @@ tRes, latRes, int_type,

double b = b_default;
double L0 = L0_default;
double q = q_default;
double ts = ts_default;
double L0_inj = L0_inj_default;
double q_inj = q_inj_default;
double t0_inj = t0_inj_default;
double ts_inj = ts_inj_default;
double n_0 = n_0_default;

@@ -772,2 +886,8 @@ double p = p_default;

double g0 = g0_default;
int envType = envType_default;
double R0_env = R0_env_default;
double k_env = k_env_default;
double rho1_env = rho1_env_default;
double E_core_global = E_core_global_default;

@@ -788,2 +908,4 @@ double theta_h_core_global = theta_h_core_global_default;

"g0",
"envType", "R0Env", "kEnv", "rho1Env",
"t0_inj",
"E0Global", "thetaCoreGlobal",

@@ -798,8 +920,11 @@ "tRes", "latRes", "intType", "rtolStruct",

if(!PyArg_ParseTupleAndKeywords(args, kwargs,
"OOO|ii""ddddddddddddddd""dd""iiidddii""O""iii",
"OOO|ii""ddddddddddddddd""iddd""d""dd""iiidddii""O""iii",
kwlist,
&theta_obj, &phi_obj, &t_obj, &jet_type, &spec_type,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b, &L0,
&q, &ts, &n_0, &p, &epsilon_E, &epsilon_B, &ksi_N, &d_L,
&theta_obs, &E_iso_core, &theta_h_core, &theta_h_wing, &b,
&L0_inj, &q_inj, &ts_inj, &n_0, &p, &epsilon_E, &epsilon_B,
&ksi_N, &d_L,
&g0,
&envType, &R0_env, &k_env, &rho1_env,
&t0_inj,
&E_core_global, &theta_h_core_global,

@@ -943,4 +1068,5 @@ &tRes, &latRes, &int_type, &rtol_struct, &rtol_phi,

theta_h_wing, b,
L0, q, ts,
L0_inj, q_inj, t0_inj, ts_inj,
n_0, p, epsilon_E, epsilon_B, ksi_N, g0,
envType, R0_env, k_env, rho1_env,
E_core_global, theta_h_core_global, ta, tb,

@@ -1000,5 +1126,6 @@ tRes, latRes, int_type,

pars.n_0 = n0;
pars.L0 = L0;
pars.q = q;
pars.ts = ts;
pars.L0_inj = L0;
pars.q_inj = q;
pars.t0_inj = t0_inj_default;
pars.ts_inj = ts;
pars.tRes = tRes;

@@ -1088,5 +1215,6 @@ pars.E_tot = -1.0;

pars.n_0 = n0;
pars.L0 = L0;
pars.q = q;
pars.ts = ts;
pars.L0_inj = L0;
pars.q_inj = q;
pars.t0_inj = t0_inj_default;
pars.ts_inj = ts;
pars.tRes = tRes;

@@ -1108,3 +1236,3 @@ pars.E_tot = -1.0;

printf("set_jet_params\n");
//printf("set_jet_params\n");
set_jet_params(&pars, E0, thetah);

@@ -1117,3 +1245,3 @@ if(pars.error)

}
printf("done\n");
//printf("done\n");

@@ -1162,2 +1290,3 @@ //Allocate output arrays

double tobs, phi, theta_obs, theta_0;
int funcVer;
PyObject *t_obj = NULL;

@@ -1169,4 +1298,4 @@ PyObject *R_obj = NULL;

//Parse Arguments
if(!PyArg_ParseTuple(args, "OOOdddd", &t_obj, &R_obj, &thj_obj, &tobs,
&phi, &theta_obs, &theta_0))
if(!PyArg_ParseTuple(args, "OOOddddi", &t_obj, &R_obj, &thj_obj, &tobs,
&phi, &theta_obs, &theta_0, &funcVer))
{

@@ -1231,6 +1360,29 @@ //PyErr_SetString(PyExc_RuntimeError, "Could not parse arguments.");

double th = find_jet_edge(phi, cos(theta_obs), sin(theta_obs), theta_0,
mu, thj, N);
double *cth = (double *)malloc(N * sizeof(double));
double *sth = (double *)malloc(N * sizeof(double));
for(i=0; i<N; i++)
{
cth[i] = cos(thj[i]);
sth[i] = sin(thj[i]);
}
int idx_mu_neg1 = searchSorted(-1.0, mu, N);
if(idx_mu_neg1 > 0)
idx_mu_neg1--;
int idx_mu_pos1 = searchSorted(1.0, mu, N);
if(idx_mu_pos1 + 1 < N)
idx_mu_pos1++;
double th;
if(funcVer)
th = find_jet_edge(phi, cos(theta_obs), sin(theta_obs), theta_0,
mu, thj, N, idx_mu_neg1, idx_mu_pos1, cth, sth);
else
th = find_jet_edge_old(phi, cos(theta_obs), sin(theta_obs), theta_0,
mu, thj, N, idx_mu_neg1, idx_mu_pos1, cth, sth);
free(mu);
free(cth);
free(sth);

@@ -1237,0 +1389,0 @@ PyObject *ret = Py_BuildValue("d", th);

+73
-32
afterglowpy/offaxis_struct.h

@@ -10,2 +10,3 @@ #ifndef AFTERGLOWPY_STRUCT

#include "integrate.h"
#include "interval.h"

@@ -56,2 +57,3 @@ #define ERR_CHK_VOID(pars) if(pars->error){ return;}

#define SIMPLE_SPEC 0
#define IC_COOLING_FLAG 1

@@ -61,2 +63,6 @@ #define EPS_E_BAR_FLAG 2

#define SSA_SHARP_FLAG 8
#define NO_COOLING_FLAG 16
#define DEEP_NEWTONIAN_FLAG 32
#define BULK_BM_FLAG 64
#define FIXED_PL_FLAG 128

@@ -70,2 +76,5 @@ enum{INT_TRAP_FIXED, INT_TRAP_ADAPT, INT_SIMP_FIXED, INT_SIMP_ADAPT,

enum{MOM_0, MOM_X, MOM_Y, MOM_Z, MOM_XX, MOM_YY, MOM_ZZ,
MOM_XY, MOM_YZ, MOM_XZ};
struct fluxParams

@@ -76,2 +85,3 @@ {

double cp;
double sp;
double ct;

@@ -86,2 +96,3 @@ double st;

double d_L;
long moment;

@@ -108,11 +119,15 @@ double E_iso;

int envType;
double As;
double Rwind;
double R0_env;
double k_env;
double rho1_env;
double L0;
double q;
double ts;
double L0_inj;
double q_inj;
double t0_inj;
double ts_inj;
double current_theta_cone_hi;
double current_theta_cone_low;
double current_theta_b;
double current_theta_a;
double theta_obs_cur;

@@ -142,2 +157,3 @@ int tRes;

double t_NR;
int cur_entry;

@@ -149,2 +165,4 @@ double *t_table;

double *mu_table;
double *cth_table;
double *sth_table;
int table_entries;

@@ -157,4 +175,14 @@

double *mu_table_inner;
double *cth_table_inner;
double *sth_table_inner;
int table_entries_inner;
int idx_mu_neg1;
int idx_mu_pos1;
int idx_mu_neg1_inner;
int idx_mu_pos1_inner;
Mesh9 phi_mesh;
Mesh9 theta_mesh;
int spec_type;

@@ -190,16 +218,32 @@ int gamma_type;

double check_t_e(double t_e, double mu, double t_obs, double *mu_table, int N);
int searchSorted(double x, double *arr, int N);
int searchSorted(double x, const double *arr, int N);
double interpolateLin(int a, int b, double x, double *X, double *Y, int N);
double interpolateLog(int a, int b, double x, double *X, double *Y, int N);
double find_jet_edge_old(double phi, double cto, double sto, double theta0,
const double *a_mu, const double *a_thj, int N,
int idx_mu_neg1, int idx_mu_pos1,
const double *a_cthj, const double *a_sthj);
double find_jet_edge(double phi, double cto, double sto, double theta0,
double *a_mu, double *a_thj, int N);
const double *a_mu, const double *a_thj, int N,
int idx_mu_neg1, int idx_mu_pos1,
const double *a_cthj, const double *a_sthj);
double costheta_integrand(double a_theta, void* params); // inner integral
double phi_integrand(double a_phi, void* params); // outer integral
double emissivity(double nu, double R, double mu, double te,
double u, double us, double n0, double p, double epse,
double epsB, double ksiN, int specType); //emissivity of
double u, double us, double rho0, double Msw, double p,
double epse, double epsB, double ksiN,
int specType); //emissivity of
// a zone.
void calc_absorption_length(double R, double mu, double delta,
double betaS, double uS,
double *length_back, double *length_front);
double absorption_integral(double Rb, double dR, double taua, double taub,
int order);
double absorption_integral_core(double a, double b, int order);
double flux(struct fluxParams *pars, double atol); // determine flux for a given t_obs
double flux_cone(double t_obs, double nu_obs, double E_iso, double theta_h,
double flux_cone(double t_obs, double nu_obs, long moment,
double E_iso, double theta_h,
double theta_cone_low, double theta_cone_hi,

@@ -246,30 +290,23 @@ double atol, struct fluxParams *pars);

struct fluxParams *pars);
void lc_tophat(double *t, double *nu, double *F, int Nt,
void lc_tophat(double *t, double *nu, double *F, long *moment, int Nt,
double E_iso, double theta_h, struct fluxParams *pars);
void lc_cone(double *t, double *nu, double *F, int Nt, double E_iso,
double theta_h, double theta_wing, struct fluxParams *pars);
void lc_powerlawCore(double *t, double *nu, double *F, int Nt,
double E_iso_core, double theta_h_core,
double theta_h_wing, double beta,
double *theta_c_arr, double *E_iso_arr,
int res_cones, struct fluxParams *pars);
void lc_powerlaw(double *t, double *nu, double *F, int Nt,
double E_iso_core, double theta_h_core,
double theta_h_wing,
double *theta_c_arr, double *E_iso_arr,
int res_cones, struct fluxParams *pars);
void lc_Gaussian(double *t, double *nu, double *F, int Nt,
void lc_cone(double *t, double *nu, double *F, long *moment, int Nt,
double E_iso, double theta_h, double theta_wing,
struct fluxParams *pars);
void lc_struct(double *t, double *nu, double *F, long *moment, int Nt,
double E_iso_core,
double theta_h_core, double theta_h_wing,
double *theta_c_arr, double *E_iso_arr,
int res_cones, struct fluxParams *pars);
void lc_GaussianCore(double *t, double *nu, double *F, int Nt,
double E_iso_core,
int res_cones, double (*f_E)(double,void *),
struct fluxParams *pars);
void lc_structCore(double *t, double *nu, double *F, long *moment, int Nt,
double E_iso_core,
double theta_h_core, double theta_h_wing,
double *theta_c_arr, double *E_iso_arr,
int res_cones, struct fluxParams *pars);
int res_cones, double (*f_E)(double,void *),
struct fluxParams *pars);
void calc_flux_density(int jet_type, int spec_type,
double *t, double *nu, double *Fnu, int N,
struct fluxParams *fp);
double *t, double *nu, double *Fnu, long *moment,
int N, struct fluxParams *fp);
void calc_intensity(int jet_type, int spec_type, double *theta, double *phi,

@@ -286,3 +323,4 @@ double *t, double *nu, double *Inu, int N,

double E_iso_core, double theta_core, double theta_wing,
double b, double L0, double q, double ts,
double b,
double L0_inj, double q_inj, double t0_inj, double ts_inj,
double n_0,

@@ -294,2 +332,3 @@ double p,

double g0,
int envType, double R0_env, double k_env, double rho1_env,
double E_core_global,

@@ -304,4 +343,6 @@ double theta_core_global,

int spread, int counterjet, int gamma_type);
void set_jet_params(struct fluxParams *pars, double E_iso, double theta_h);
void set_obs_params(struct fluxParams *pars, double t_obs, double nu_obs,
void set_obs_params(struct fluxParams *pars,
double t_obs, double nu_obs, long moment,
double theta_obs_cur, double current_theta_cone_hi,

@@ -308,0 +349,0 @@ double current_theta_cone_low);

+344
-48
afterglowpy/shockEvolution.c

@@ -73,2 +73,79 @@ #include <stdio.h>

double envDensityPar(double R, struct shockParams *par)
{
return envDensity(R, par->envType, par->rho0_env, par->R0_env, par->k_env,
par->rho1_env);
}
double envMassPar(double R, struct shockParams *par)
{
return envMass(R, par->envType, par->rho0_env, par->R0_env, par->k_env,
par->rho1_env);
}
double envRadiusPar(double M, struct shockParams *par)
{
return envRadius(M, par->envType, par->rho0_env, par->R0_env, par->k_env,
par->rho1_env);
}
double envDensity(double R, int envType, double rho0, double R0, double k,
double rho1)
{
if(envType == ENV_ISM)
return rho0;
else if(envType == ENV_WIND)
return rho0 * R0*R0/(R*R);
else if(envType == ENV_PL)
return rho0 * pow(R/R0, -k);
else if(envType == ENV_STEP)
return R < R0 ? rho0 : rho1;
else
return 0;
}
double envMass(double R, int envType, double rho0, double R0, double k,
double rho1)
{
if(envType == ENV_ISM)
return 4.0/3.0 * M_PI * rho0 * R*R*R;
else if(envType == ENV_WIND)
return 4.0*M_PI * rho0 * R0*R0 * R;
else if(envType == ENV_PL)
return 4.0*M_PI / (3.0-k) * rho0 * pow(R/R0, 3-k) * R0*R0*R0;
else if (envType == ENV_STEP)
{
if(R < R0)
return 4.0/3.0 * M_PI * rho0 * R*R*R;
else
return 4.0/3.0 * M_PI * (rho0 * R0*R0*R0
+ rho1 * (R-R0)*(R-R0)*(R-R0));
}
else
return 0;
}
double envRadius(double M, int envType, double rho0, double R0, double k,
double rho1)
{
if(envType == ENV_ISM)
return pow(0.75 * M / (M_PI * rho0), 1.0/3.0);
else if(envType == ENV_WIND)
return M / (4.0*M_PI * rho0 * R0*R0);
else if(envType == ENV_PL)
return R0 * pow((3.0-k) * M / (4.0*M_PI * rho0 * R0*R0*R0), 1.0/(3-k));
else if(envType == ENV_STEP)
{
double M0 = 4.0/3.0 * M_PI * rho0 * R0*R0*R0;
if(M <= M0)
return pow(0.75 * M / (M_PI * rho0), 1.0/3.0);
else
return R0 + pow(0.75 * (M-M0) / (M_PI * rho1), 1.0/3.0);
}
else
return 0;
}
void shockInitDecel(double t0, double *R0, double *u0, void *argv)

@@ -176,10 +253,130 @@ {

{
double *args = (double *)argv;
double E0 = args[0];
double Mej = args[1];
double rho0 = args[2];
double L0 = args[6];
double q = args[7];
double ts = args[8];
struct shockParams *par = (struct shockParams *)argv;
if(par->envType == ENV_ISM)
{
shockInitFindISM(t0, R0, u0, tRes, argv);
return;
}
double E0 = par->E0;
double Mej = par->Mej;
double c2 = v_light * v_light;
double t00, R00, u00;
if(Mej > 0.0)
{
double gm1_coast = E0 / (Mej * c2);
double g_coast = gm1_coast + 1;
double u2_coast = gm1_coast * (gm1_coast + 2);
double be2_coast = u2_coast / (g_coast * g_coast);
double Msw = 3.0 * E0 / ((4*u2_coast + 3) * be2_coast * c2);
//Rough deceleration radius
double Rd = envRadiusPar(Msw, par);
//shock 4-velocity
double u_coast = sqrt(u2_coast);
double uS = shockVel(u_coast);
double beS = uS / sqrt(uS*uS + 1);
double td = Rd / (beS * v_light);
if(t0 < 0.001*td)
{
*R0 = beS * v_light * t0;
*u0 = u_coast;
return;
}
t00 = 0.001*td;
R00 = beS * v_light * t00;
u00 = u_coast;
}
else
{
//Setting ref ultra-relativistic velocity
double u_UR = 1000;
double g2_UR = u_UR*u_UR + 1;
double be2_UR = u_UR*u_UR / g2_UR;
//swept mass at u_UR;
double M_UR = 3.0 * E0 / ((4*u_UR*u_UR + 3) * be2_UR * c2);
//Rough UR radius
double R_UR = envRadiusPar(M_UR, par);
double rho_UR = envDensityPar(R_UR, par);
// = 3-k for a power law rho ~ R^-k
double volIdx = 4*M_PI * R_UR*R_UR*R_UR * rho_UR / M_UR;
double t_UR = R_UR/v_light * (1.0 + 1.0/(4*(1+volIdx) * g2_UR));
if(t0 < t_UR)
{
double R_apprx = v_light * t0;
double M0 = envMassPar(R_apprx, par);
double g2 = 0.75 * E0 / (M0 * c2);
*R0 = R_apprx * (1 - 1.0/(4 * (1+volIdx) * g2));
*u0 = sqrt(g2 - 1.0);
return;
}
t00 = t_UR;
R00 = R_UR;
u00 = u_UR;
}
double dt, x[2], x0[2], k1[2], k2[2], k3[2], k4[2];
x[0] = R00;
x[1] = u00;
double t = t00;
double fac = pow(10, 1.0/tRes);
int j;
while(t < t0)
{
dt = (fac-1)*t;
if(fac*t >= t0)
dt = t0-t;
x0[0] = x[0];
x0[1] = x[1];
Rudot2D(t, x, argv, k1);
for(j=0; j<2; j++)
x[j] = x0[j] + 0.5*dt*k1[j];
Rudot2D(t, x, argv, k2);
for(j=0; j<2; j++)
x[j] = x0[j] + 0.5*dt*k2[j];
Rudot2D(t, x, argv, k3);
for(j=0; j<2; j++)
x[j] = x0[j] + dt*k3[j];
Rudot2D(t, x, argv, k4);
for(j=0; j<2; j++)
x[j] = x0[j] + dt*(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6.0;
t *= fac;
}
*R0 = x[0];
*u0 = x[1];
}
void shockInitFindISM(double t0, double *R0, double *u0, double tRes,
void *argv)
{
struct shockParams *par = (struct shockParams *)argv;
double E0 = par->E0;
double Mej = par->Mej;
double rho0 = par->rho0_env;
double L0 = par->L0_inj;
double q = par->q_inj;
double ts = par->ts_inj;
//printf("shockInitFind: t0=%.6lg E0=%.6lg Mej=%.6lg rho0=%.6lg\n",

@@ -285,15 +482,15 @@ // t0, E0, Mej, rho0);

x0[1] = x[1];
Rudot2D(t, x, args, k1);
Rudot2D(t, x, argv, k1);
for(j=0; j<2; j++)
x[j] = x0[j] + 0.5*dt*k1[j];
Rudot2D(t, x, args, k2);
Rudot2D(t, x, argv, k2);
for(j=0; j<2; j++)
x[j] = x0[j] + 0.5*dt*k2[j];
Rudot2D(t, x, args, k3);
Rudot2D(t, x, argv, k3);
for(j=0; j<2; j++)
x[j] = x0[j] + dt*k3[j];
Rudot2D(t, x, args, k4);
Rudot2D(t, x, argv, k4);

@@ -311,12 +508,11 @@ for(j=0; j<2; j++)

{
double *args = (double *)argv;
struct shockParams *par = (struct shockParams *)argv;
double Mej = args[1];
double rho0 = args[2];
double Einj = args[3];
double k = args[4];
double umin = args[5];
double L0 = args[6];
double q = args[7];
double ts = args[8];
double Mej = par->Mej;
double Einj = par->E0_refresh;
double k = par->k_refresh;
double umin = par->umin_refresh;
double L0 = par->L0_inj;
double q = par->q_inj;
double ts = par->ts_inj;

@@ -345,8 +541,18 @@ double R = x[0];

double num = -16*M_PI/3.0 * rho0*R*R * be*u*u * v_light
+ dEdt/(v_light*v_light);
double denom = be*Mej + 8*M_PI*rho0*R*R*R*u*(2*u*u+1)*(2*u*u+3)/(9*g*g*g*g)
- dEdu/(v_light*v_light);
double dudt = num/denom;
double rho = envDensityPar(R, par);
double M = envMassPar(R, par);
double dR_Esh = 4*M_PI/3.0 * (4*u*u+3)*be*be * rho * R*R * v_light*v_light;
double du_Esh = 2.0 * (2*u*u+1)*(2*u*u+3)*u * M * v_light*v_light
/ (3.0 * g*g*g*g);
double du_Eej = be * Mej * v_light*v_light;
double dudt = (-dR_Esh * dRdt + dEdt) / (du_Eej + du_Esh - dEdu);
//double num = -16*M_PI/3.0 * rho0*R*R * be*u*u * v_light
// + dEdt/(v_light*v_light);
//double denom = be*Mej
// + 8*M_PI*rho0*R*R*R*u*(2*u*u+1)*(2*u*u+3)/(9*g*g*g*g)
// - dEdu/(v_light*v_light);
//double dudt = num/denom;
xdot[0] = dRdt;

@@ -356,17 +562,18 @@ xdot[1] = dudt;

void RuThdot3D(double t, double *x, void *argv, double *xdot, int spread)
void RuThdot3D(double t, double *x, void *argv, double *xdot)
{
double *args = (double *)argv;
struct shockParams *par = (struct shockParams *)argv;
double Mej = args[1];
double rho0 = args[2];
double Einj = args[3];
double k = args[4];
double umin = args[5];
double L0 = args[6];
double q = args[7];
double ts = args[8];
double thC = args[9];
double th0 = args[10];
double thCg = args[11];
double Mej = par->Mej;
//double rho0 = par->rho0_env;
double Einj = par->E0_refresh;
double k = par->k_refresh;
double umin = par->umin_refresh;
double L0 = par->L0_inj;
double q = par->q_inj;
double ts = par->ts_inj;
double thC = par->thetaCore;
double th0 = par->theta0;
double thCg = par->thetaCoreGlobal;
int spread = par->spread;

@@ -387,2 +594,6 @@ double R = x[0];

double rho = envDensityPar(R, par);
double M = envMassPar(R, par);
double env3mk = 4*M_PI*rho*R*R*R / M;
double dThdt = 0.0;

@@ -500,4 +711,9 @@ //TODO: verify spreading procedure 190401

{
double Q0 = 2.0;
double Q = sqrt(2.0)*3.0;
double envk = 3 - env3mk;
double Q = sqrt(2.0)*env3mk;
double Q0;
if(envk < 2.0)
Q0 = 0.5 * (2.0 * (2.0 - envk) + 1.4 * envk);
else
Q0 = 1.4 * (3.0 - envk);

@@ -515,4 +731,9 @@ if(th < 0.5*M_PI && Q0*u*thC < 1)

{
double Q0 = 2.0;
double Q = sqrt(2.0)*3.0;
double envk = 3 - env3mk;
double Q = sqrt(2.0)*env3mk;
double Q0;
if(envk < 2.0)
Q0 = 0.5 * (2.0 * (2.0 - envk) + 1.4 * envk);
else
Q0 = 1.4 * (3.0 - envk);

@@ -528,2 +749,37 @@ if(th < 0.5*M_PI && Q0*u*thCg < 1)

}
else if(spread == 9)
{
//Fernandez & Lamb 2021 spreading.
if(th < 0.5*M_PI && u*thC < 100)
{
double gs = 1 / sqrt((1-bes)*(1+bes));
double fac = 1.0;
if (th0 < thC)
fac *= tan(0.5*th0)/tan(0.5*thC); //th0/thC;
dThdt = fac / (gs*gs * th) * dRdt / R;
}
}
else if(spread == 11)
{
double envk = 3 - env3mk;
double Q = sqrt(2.0)*env3mk;
double Q0;
if(envk < 2.0)
Q0 = 0.5 * (2.0 * (2.0 - envk) + 1.4 * envk);
else
Q0 = 1.4 * (3.0 - envk);
if(th < 0.5*M_PI && Q0*u*thC < 1)
{
double bew = 0.5*sqrt((2*u*u+3)/(4*u*u+3))*bes/g;
double fac = u*thC*Q < 1.0 ? 1.0 : Q*(1-Q0*u*thC) / (Q-Q0);
if (th0 < thC)
fac *= tan(0.5*th0)/tan(0.5*thC); //th0/thC;
double gs = 1.0 / sqrt((1-bes)*(1+bes));
double dThdt_lum = fac * v_light / (R * g);
dThdt = 4.0 * fac * bew * v_light / R;
if(dThdt > dThdt_lum)
dThdt = dThdt_lum;
}
}
}

@@ -543,2 +799,12 @@

double c2 = v_light*v_light;
double dR_Esh = 4*M_PI/3.0 * (4*u*u+3)*be*be * rho * R*R * om * c2;
double du_Esh = 2.0*(2*u*u+1)*(2*u*u+3)*u * M * om * c2 / (3.0 * g*g*g*g);
double dTh_Esh = (4*u*u+3)*be*be * M * 2*sinth*costh * c2 / 3.0;
double du_Eej = be * Mej * c2;
double dudt = (-dR_Esh * dRdt - dTh_Esh * dThdt + dEdt)
/ (du_Eej + du_Esh - dEdu);
/*
double num = -16*M_PI/3.0 * om*rho0*R*R * be*u*u * v_light

@@ -551,2 +817,3 @@ -8*M_PI/9.0*rho0*R*R*R*(4*u*u+3)*be*be*sinth*costh*dThdt

double dudt = num/denom;
*/

@@ -595,4 +862,3 @@ xdot[0] = dRdt;

void shockEvolveSpreadRK4(double *t, double *R, double *u, double *th, int N,
double R0, double u0, double th0, void *args,
int spread)
double R0, double u0, double th0, void *args)
{

@@ -613,15 +879,15 @@ int i,j;

x0[2] = th[i];
RuThdot3D(t[i], x0, args, k1, spread);
RuThdot3D(t[i], x0, args, k1);
for(j=0; j<3; j++)
x[j] = x0[j] + 0.5*dt*k1[j];
RuThdot3D(t[i], x, args, k2, spread);
RuThdot3D(t[i], x, args, k2);
for(j=0; j<3; j++)
x[j] = x0[j] + 0.5*dt*k2[j];
RuThdot3D(t[i], x, args, k3, spread);
RuThdot3D(t[i], x, args, k3);
for(j=0; j<3; j++)
x[j] = x0[j] + dt*k3[j];
RuThdot3D(t[i], x, args, k4, spread);
RuThdot3D(t[i], x, args, k4);

@@ -639,1 +905,31 @@ for(j=0; j<3; j++)

}
void setup_shockParams(struct shockParams *pars, int spread,
double E0, double Mej,
int envType, double rho0_env, double R0_env,
double k_env, double rho1_env,
double L0_inj, double q_inj, double t0_inj,
double ts_inj,
double E0_refresh, double k_refresh, double umin_refresh,
double thetaCore, double theta0,
double thetaCoreGlobal)
{
pars->spread = spread;
pars->E0 = E0;
pars->Mej = Mej;
pars->envType = envType;
pars->rho0_env = rho0_env;
pars->R0_env = R0_env;
pars->k_env = k_env;
pars->rho1_env = rho1_env;
pars->L0_inj = L0_inj;
pars->q_inj = q_inj;
pars->t0_inj = t0_inj;
pars->ts_inj = ts_inj;
pars->E0_refresh = E0_refresh;
pars->k_refresh = k_refresh;
pars->umin_refresh = umin_refresh;
pars->thetaCore = thetaCore;
pars->theta0 = theta0;
pars->thetaCoreGlobal = thetaCoreGlobal;
}
+47
-3
afterglowpy/shockEvolution.h
#ifndef AFTERGLOWPY_SHOCK
#define AFTERGLOWPY_SHOCK
enum{ENV_ISM, ENV_WIND, ENV_PL, ENV_STEP};
static const double T0_inj = 1.0e3;
struct shockParams
{
int spread;
double E0;
double Mej;
double L0_inj;
double q_inj;
double t0_inj;
double ts_inj;
double E0_refresh;
double k_refresh;
double umin_refresh;
double thetaCore;
double theta0;
double thetaCoreGlobal;
int envType;
double rho0_env;
double R0_env;
double k_env;
double rho1_env;
};
void setup_shockParams(struct shockParams *pars, int spread,
double E0, double Mej,
int envType, double env_rho0, double R0_env,
double k_env, double rho1_env,
double L0_inj, double q_inj, double t0_inj,
double ts_inj,
double E0_refresh, double k_refresh, double umin_refresh,
double thetaCore, double theta0,
double thetaCoreGlobal);
double shockVel(double u);
double E_inj(double te, double L0, double q, double ts);
double L_inj(double te, double L0, double q, double ts);
double envDensityPar(double R, struct shockParams *pars);
double envMassPar(double R, struct shockParams *pars);
double envRadiusPar(double M, struct shockParams *pars);
double envDensity(double R, int envType, double rho0, double R0, double k,
double rho1);
double envMass(double R, int envType, double rho0, double R0, double k,
double rho1);
double envRadius(double M, int envType, double rho0, double R0, double k,
double rho1);
void shockInitDecel(double t0, double *R0, double *u0, void *argv);
void shockInitFind(double t0, double *R0, double *u0, double tRes, void *argv);
void shockInitFindISM(double t0, double *R0, double *u0, double tRes,
void *argv);
void Rudot2D(double t, double *x, void *argv, double *xdot);
void RuThdot3D(double t, double *x, void *argv, double *xdot, int spread);
void RuThdot3D(double t, double *x, void *argv, double *xdot);
void shockEvolveRK4(double *t, double *R, double *u, int N, double R0,
double u0, void *args);
void shockEvolveSpreadRK4(double *t, double *R, double *u, double *th, int N,
double R0, double u0, double th0, void *args,
int spread);
double R0, double u0, double th0, void *args);
#endif
+37
-11
afterglowpy/shockmodule.c

@@ -103,7 +103,12 @@ #include <Python.h>

double R0, u0;
double Mej, rho0, Einj, k, umin, L0, q, ts;
double Mej, rho0_env, R0_env, k_env, rho1_env, Einj, k, umin, L0, q, ts;
int envType;
//Parse Arguments
if(!PyArg_ParseTuple(args, "Odddddddddd", &t_obj, &R0, &u0, &Mej, &rho0,
&Einj, &k, &umin, &L0, &q, &ts))
if(!PyArg_ParseTuple(args, "Odd""d""idddd""ddd""ddd",
&t_obj, &R0, &u0,
&Mej,
&envType, &rho0_env, &R0_env, &k_env, &rho1_env,
&Einj, &k, &umin,
&L0, &q, &ts))
{

@@ -155,4 +160,10 @@ PyErr_SetString(PyExc_RuntimeError, "Could not parse arguments.");

// Evolve the shock!
double shockArgs[9] = {u0, Mej, rho0, Einj, k, umin, L0, q, ts};
shockEvolveRK4(t, R, u, N, R0, u0, shockArgs);
struct shockParams shock_pars;
setup_shockParams(&shock_pars, 0, 0.0, Mej,
envType, rho0_env, R0_env, k_env, rho1_env,
L0, q, 1.0e3, ts,
Einj, k, umin,
0.0, 0.0, 0.0);
//double shockArgs[9] = {0.0, Mej, rho0, Einj, k, umin, L0, q, ts};
shockEvolveRK4(t, R, u, N, R0, u0, &shock_pars);

@@ -173,8 +184,15 @@ // Clean up!

double R0, u0, th0;
double Mej, rho0, Einj, k, umin, L0, q, ts, thC;
double Mej, rho0_env, R0_env, k_env, rho1_env;
double Einj, k, umin, L0, q, ts, thC;
int envType;
int spread;
//Parse Arguments
if(!PyArg_ParseTuple(args, "Oddddddddddddi", &t_obj, &R0, &u0, &th0,
&Mej, &rho0, &Einj, &k, &umin, &L0, &q, &ts,
if(!PyArg_ParseTuple(args, "Oddd""d""idddd""ddd""ddd""di",
&t_obj, &R0, &u0, &th0,
&Mej,
&envType, &rho0_env, &R0_env, &k_env,
&rho1_env,
&Einj, &k, &umin,
&L0, &q, &ts,
&thC, &spread))

@@ -230,6 +248,14 @@ {

// Evolve the shock!
double shockArgs[12] = {u0, Mej, rho0, Einj, k, umin, L0, q, ts, thC, th0,
thC};
shockEvolveSpreadRK4(t, R, u, th, N, R0, u0, th0, shockArgs, spread);
//double shockArgs[12] = {0.0, Mej, rho0, Einj, k, umin, L0, q, ts, thC, th0,
// thC};
struct shockParams shock_pars;
setup_shockParams(&shock_pars, spread, 0.0, Mej,
envType, rho0_env, R0_env, k_env, rho1_env,
L0, q, 1.0e3, ts,
Einj, k, umin,
thC, th0, thC);
shockEvolveSpreadRK4(t, R, u, th, N, R0, u0, th0, &shock_pars);
// Clean up!

@@ -236,0 +262,0 @@ Py_DECREF(t_arr);

+2
-2
afterglowpy/version.py
#!/usr/bin/env python3
"""Version info"""
__short_version__ = '0.7'
__version__ = '0.7.3'
__short_version__ = '0.8'
__version__ = '0.8.0'
+95
-104
PKG-INFO
Metadata-Version: 2.1
Name: afterglowpy
Version: 0.7.3
Version: 0.8.0
Summary: GRB Afterglow Models
Home-page: https://github.com/geoffryan/afterglowpy
Author: Geoffrey Ryan
Author-email: gsryan@umd.edu
Author-email: gryan@perimeterinstitute.ca
License: MIT
Project-URL: Source Code, https://github.com/geoffryan/afterglowpy
Project-URL: Documentation, https://afterglowpy.readthedocs.io
Description: # Numeric GRB Afterglow models
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2019](https://arxiv.org/abs/1909.11691). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.
Documentation available at <https://afterglowpy.readthedocs.io/>
## Attribution
If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2019, arXiv:1910.11691" [arXiv link](https://arxiv.org/abs/1909.11691).
## Features
_afterglowpy_ computes synchrotron emission from the forward shock of a relativistic blast wave. It includes:
- Fully trans-relativistic shock evolution through a constant density medium.
- On-the-fly integration over the equal-observer-time slices of the shock surface.
- Approximate prescription for jet spreading.
- Arbitrary viewing angles.
- Angularly structured jets, ie. E(&theta;)
- Spherical velocity-stratified outflows, ie. E(u)
- Counter-jet emission.
It has limited support (these should be considered experimental) for:
- Initial energy injection
- Inverse comption spectra
- Spatially-resolved intensity maps
- Early coasting phase
It does not include (yet):
- External wind medium, ie. n &prop; r<sup>-2</sup>
- Synchrotron self-absorbtion
- Reverse shock emission
_afterglowpy_ has been calibrated to the BoxFit code ([van Eerten, van der Horst, & Macfadyen 2011](https://arxiv.org/abs/1110.5089), available at the [Afterglow Library](https://cosmo.nyu.edu/afterglowlibrary/boxfit2011.html)) and produces similar light curves for top hat jets (within 50% when same parameters are used) both on- and off-axis. Its jet models by default do not include an initial coasting phase, which may effect predictions for early observations.
## Installation/Building
_afterglowpy_ is available via `pip`:
```bash
$ pip install afterglowpy
```
If you are working on a local copy of this repo and would like to install from source, you can the run the following from the top level directory of the project.
```bash
$ pip install -e .
```
## Using
**This interface will be updated to be more sensible in the VERY near future**
In your python code, import the library with `import afterglowpy as grb`.
The main function of interest is`grb.fluxDensity(t, nu, jetType, specType, *pars, **kwargs)`. See `tests/plotLC.py` for a simple example.
`jetType` can be -1 (top hat), 0 (Gaussian), 1 (Power Law w/ core), 2 (Gaussian w/ core), 3 (Cocoon), or 4 (Smooth Power Law).
`specType` can be 0 (global cooling time, no inverse compton) or 1 (global cooling time, inverse compton).
For jet-like afterglows (`jetTypes` -2, -1, 0, 1, 2, and 4) `pars` has 14 positional arguments:
- `0 thetaV` viewing angle in radians
- `1 E0` on-axis isotropic equivalent energy in erg
- `2 thetaC` half-width of the jet core in radians (jetType specific)
- `3 thetaW` "wing" truncation angle of the jet, in radians
- `4 b` power for power-law structure, &theta;<sup>-b</sup>
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
For cocoon-like afterglows (`jetType` 3) `pars` has 14 positional arguments:
- `0 umax` Initial maximum outflow 4-velocity
- `1 umin` Minium outflow 4-velocity
- `2 Ei` Fiducial energy in velocity distribution, E(>u) = E<sub>i</sub> u<sup>-k</sup>.
- `3 k` Power-law index of energy velocity distribution
- `4 Mej` Mass of material at `umax' in solar masses
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
Keyword arguments are:
- `z` redshift (defaults to 0)
- `tRes` time resolution of shock-evolution scheme, number of sample points per decade in time
- `latRes` latitudinal resolution for structured jets, number of shells per `thetaC`
- `rtol` target relative tolerance of flux integration
- `spread` boolean (defaults to True), whether to allow the jet to spread.
Platform: UNKNOWN
Classifier: Programming Language :: Python :: 3

@@ -121,2 +19,95 @@ Classifier: Programming Language :: C

Description-Content-Type: text/markdown
License-File: LICENSE.txt
Requires-Dist: numpy>=1.13
Requires-Dist: scipy>=0.14
Provides-Extra: docs
Requires-Dist: numpydoc; extra == "docs"
# Numeric GRB Afterglow models
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2020](https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.
Documentation available at <https://afterglowpy.readthedocs.io/>
## Attribution
If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2020, Astrophysical Journal *896*, 166 (2020)" [arXiv link](https://arxiv.org/abs/1909.11691).
## Acknowledgements
This work is funded in part by the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158).
## Features
_afterglowpy_ computes synchrotron emission from the forward shock of a relativistic blast wave. It includes:
- Fully trans-relativistic shock evolution through a constant density medium.
- On-the-fly integration over the equal-observer-time slices of the shock surface.
- Approximate prescription for jet spreading.
- Arbitrary viewing angles.
- Angularly structured jets, ie. E(&theta;)
- Spherical velocity-stratified outflows, ie. E(u)
- Counter-jet emission.
- Deep Newtonian emission.
- Image moments suitable for astrometry: centroid position and image size.
It has limited support (these should be considered experimental) for:
- Initial energy injection
- Inverse comption spectra
- Early coasting phase
It does not include (yet):
- External wind medium, ie. n &prop; r<sup>-2</sup>
- Synchrotron self-absorbtion
- Reverse shock emission
_afterglowpy_ has been calibrated to the BoxFit code ([van Eerten, van der Horst, & Macfadyen 2011](https://arxiv.org/abs/1110.5089), available at the [Afterglow Library](https://cosmo.nyu.edu/afterglowlibrary/boxfit2011.html)) and produces similar light curves for top hat jets (within 50% when same parameters are used) both on- and off-axis. Its jet models by default do not include an initial coasting phase, which may effect predictions for early observations.
## Installation/Building
_afterglowpy_ is available via `pip`:
```bash
$ pip install afterglowpy
```
If you are working on a local copy of this repo and would like to install from source, you can the run the following from the top level directory of the project.
```bash
$ pip install -e .
```
## Using
In your python code, import the library with `import afterglowpy as grb`.
The main function of interest is`grb.fluxDensity(t, nu, **kwargs)`. See `examples/plotLightCurve.py` for a simple example.
For jet-like afterglows there are up to 13 required keyword arguments:
- `jetType` an integer code setting the jet structure. It can be `grb.jet.TopHat`, `grb.jet.Gaussian`, `grb.jet.PowerLawCore`, `grb.jet.GaussianCore`, `grb.jet.Spherical`, or `grb.jet.PowerLaw`.
- `specType` an integer code specifying flags for the emissivity function and spectrum. Can be `grb.jet.SimpleSpec` (basic spectrum with &nu;<sub>m</sub> and &nu;<sub>c</sub>), `grb.jet.DeepNewtonian`, `grb.jet.ICCooling` (simple inverse Compton effects on the cooling frequency, experimental).
- `thetaObs` viewing angle in radians
- `E0` on-axis isotropic equivalent energy in erg
- `thetaCore` half-width of the jet core in radians (jetType specific)
- `thetaWing` "wing" truncation angle of the jet, in radians
- `b` power for power-law structure, &theta;<sup>-b</sup>
- `n0` Number density of ISM, in cm<sup>-3</sup>
- `p` Electron distribution power-law index (p>2)
- `epsilon_e` Thermal energy fraction in electrons
- `epsilon_B` Thermal energy fraction in magnetic field
- `xi_N` Fraction of electrons that get accelerated
- `d_L` Luminosity distance in cm
Optional keyword arguments for all models are:
- `z` redshift (defaults to 0)
- `spread` boolean (defaults to True), whether to allow the jet to spread.
- `counterjet` boolean (defaults to False), whether to include the counterjet
- `moment` array (integer dtype, same shape as t and nu) which sky moment to compute.
- `L0` Fiducial luminosity for energy injection, in erg/s, default 0.0.
- `q` Temporal power-law index for energy injection, default 0.0.
- `ts` Fiducial time-scale for energy injection, in seconds, default 0.
- `tRes` time resolution of shock-evolution scheme, number of sample points per decade in time
- `latRes` latitudinal resolution for structured jets, number of shells per `thetaC`
- `rtol` target relative tolerance of flux integration
+30
-42
README.md
# Numeric GRB Afterglow models
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2019](https://arxiv.org/abs/1909.11691). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.
A Python 3 module to calculate GRB afterglow light curves and spectra. Details of the methods can be found in [Ryan et al 2020](https://ui.adsabs.harvard.edu/abs/2020ApJ...896..166R/abstract). Builds on [van Eerten & MacFadyen 2010](https://arxiv.org/abs/1006.5125) and [van Eerten 2018](https://arxiv.org/abs/1801.01848). This code is under active development.

@@ -9,4 +9,8 @@ Documentation available at <https://afterglowpy.readthedocs.io/>

If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2019, arXiv:1910.11691" [arXiv link](https://arxiv.org/abs/1909.11691).
If you use this code in a publication, please refer to the package by name and cite "Ryan, G., van Eerten, H., Piro, L. and Troja, E., 2020, Astrophysical Journal *896*, 166 (2020)" [arXiv link](https://arxiv.org/abs/1909.11691).
## Acknowledgements
This work is funded in part by the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158).
## Features

@@ -22,2 +26,4 @@

- Counter-jet emission.
- Deep Newtonian emission.
- Image moments suitable for astrometry: centroid position and image size.

@@ -27,3 +33,2 @@ It has limited support (these should be considered experimental) for:

- Inverse comption spectra
- Spatially-resolved intensity maps
- Early coasting phase

@@ -53,52 +58,35 @@

**This interface will be updated to be more sensible in the VERY near future**
In your python code, import the library with `import afterglowpy as grb`.
The main function of interest is`grb.fluxDensity(t, nu, jetType, specType, *pars, **kwargs)`. See `tests/plotLC.py` for a simple example.
The main function of interest is`grb.fluxDensity(t, nu, **kwargs)`. See `examples/plotLightCurve.py` for a simple example.
`jetType` can be -1 (top hat), 0 (Gaussian), 1 (Power Law w/ core), 2 (Gaussian w/ core), 3 (Cocoon), or 4 (Smooth Power Law).
For jet-like afterglows there are up to 13 required keyword arguments:
`specType` can be 0 (global cooling time, no inverse compton) or 1 (global cooling time, inverse compton).
- `jetType` an integer code setting the jet structure. It can be `grb.jet.TopHat`, `grb.jet.Gaussian`, `grb.jet.PowerLawCore`, `grb.jet.GaussianCore`, `grb.jet.Spherical`, or `grb.jet.PowerLaw`.
- `specType` an integer code specifying flags for the emissivity function and spectrum. Can be `grb.jet.SimpleSpec` (basic spectrum with &nu;<sub>m</sub> and &nu;<sub>c</sub>), `grb.jet.DeepNewtonian`, `grb.jet.ICCooling` (simple inverse Compton effects on the cooling frequency, experimental).
- `thetaObs` viewing angle in radians
- `E0` on-axis isotropic equivalent energy in erg
- `thetaCore` half-width of the jet core in radians (jetType specific)
- `thetaWing` "wing" truncation angle of the jet, in radians
- `b` power for power-law structure, &theta;<sup>-b</sup>
- `n0` Number density of ISM, in cm<sup>-3</sup>
- `p` Electron distribution power-law index (p>2)
- `epsilon_e` Thermal energy fraction in electrons
- `epsilon_B` Thermal energy fraction in magnetic field
- `xi_N` Fraction of electrons that get accelerated
- `d_L` Luminosity distance in cm
For jet-like afterglows (`jetTypes` -2, -1, 0, 1, 2, and 4) `pars` has 14 positional arguments:
- `0 thetaV` viewing angle in radians
- `1 E0` on-axis isotropic equivalent energy in erg
- `2 thetaC` half-width of the jet core in radians (jetType specific)
- `3 thetaW` "wing" truncation angle of the jet, in radians
- `4 b` power for power-law structure, &theta;<sup>-b</sup>
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
For cocoon-like afterglows (`jetType` 3) `pars` has 14 positional arguments:
- `0 umax` Initial maximum outflow 4-velocity
- `1 umin` Minium outflow 4-velocity
- `2 Ei` Fiducial energy in velocity distribution, E(>u) = E<sub>i</sub> u<sup>-k</sup>.
- `3 k` Power-law index of energy velocity distribution
- `4 Mej` Mass of material at `umax' in solar masses
- `5 L0` Fiducial luminosity for energy injection, in erg/s, typically 0.
- `6 q` Temporal power-law index for energy injection, typically 0.
- `7 ts` Fiducial time-scale for energy injection, in seconds, typically 0.
- `8 n0` Number density of ISM, in cm<sup>-3</sup>
- `9 p` Electron distribution power-law index (p>2)
- `10 epsilon_e` Thermal energy fraction in electrons
- `11 epsilon_B` Thermal energy fraction in magnetic field
- `12 xi_N` Fraction of electrons that get accelerated
- `13 d_L` Luminosity distance in cm
Keyword arguments are:
Optional keyword arguments for all models are:
- `z` redshift (defaults to 0)
- `spread` boolean (defaults to True), whether to allow the jet to spread.
- `counterjet` boolean (defaults to False), whether to include the counterjet
- `moment` array (integer dtype, same shape as t and nu) which sky moment to compute.
- `L0` Fiducial luminosity for energy injection, in erg/s, default 0.0.
- `q` Temporal power-law index for energy injection, default 0.0.
- `ts` Fiducial time-scale for energy injection, in seconds, default 0.
- `tRes` time resolution of shock-evolution scheme, number of sample points per decade in time
- `latRes` latitudinal resolution for structured jets, number of shells per `thetaC`
- `rtol` target relative tolerance of flux integration
- `spread` boolean (defaults to True), whether to allow the jet to spread.
+2
-2
setup.py

@@ -43,3 +43,3 @@ #!/usr/bin/env python3

author="Geoffrey Ryan",
author_email="gsryan@umd.edu",
author_email="gryan@perimeterinstitute.ca",
description='GRB Afterglow Models',

@@ -60,3 +60,3 @@ long_description=long_description,

"Topic :: Scientific/Engineering :: Astronomy"],
install_requires=['numpy>=1.10', 'scipy>=0.14'],
install_requires=['numpy>=1.13', 'scipy>=0.14'],
extras_require={

@@ -63,0 +63,0 @@ 'docs': ['numpydoc']

afterglowpy/offaxis_struct_funcs.c

Sorry, the diff of this file is too big to display