Response

class cosi_atmosphere.response.Atmosphere

get_atm_profile(filename, dates, lons, lats, alts, version=2.1)

Get atmospher model from NRLMSIS: Naval Research Laboratory’s Mass Spectrometer Incoherent Scatter Radar model.

For pymsis homepage see: https://swxtrec.github.io/pymsis/index.html

Parameters:

• filename (str) – Name of output data file.

• dates (ArrayLike) – Dates and time of interest.

• lons (ArrayLike:) – Longitudes of interest in degrees.

• lats (ArrayLike) – Latitudes of interest in degrees.

• alts (ArrayLike) – Altitudes of interest in km.

• version (float, optional) – MSIS version number, one of (0,2.0,2.1). Default is 2.1.

Note

All inputs take arrays; however, for now use single values for dates, lons, and lats. The only array should be alts, since for this method we mainly want the altitude profile.

Returns:

atm_dict – Dictionary with atmospheric information. Keys are altitude[km], mass_density[kg/m3], N2[m-3], O2[m-3], O[m-3], He[m-3], H[m-3], Ar[m-3], N[m-3], anomalous_oxygen[m-3], NO[m-3], and Temperature[k].

Return type:

dict

class cosi_atmosphere.response.MakeMassModels(atmosphere_file, kwargs={})

Generates a mass model based on input atmospheric data.

Parameters:

• atmosphere_file (str) – Input file describing the atmosphere, calculated with Atmospheric_profile class, based on NRLMSIS.

• kwargs (dict, optional) – Pass any kwargs to pandas read_csv method.

get_cart_vectors(angle, altitude, beam_alt=200)

Get x position and direction vectors for off-axis beam.

Parameters:

• angle (float) – Incident angle of source in degrees.

• altitude (float) – Altitude of detecting plane in km.

• beam_alt (float, optional) – Altitude at which the beam is placed in km (default is 200 km).

plot_atmosphere(show_comparison=False)

Plot atmosphere model.

Parameters:

show_comparison (bool, optional) – If true, compare normalized atmospheric components to original calculation from AZ (default is False).

rectangular_model(watch_height)

Original rectangular mass model by AZ, from MEGAlib/resource/examples/advanced/Atmosphere.

This is also similar to what was done by Alex Lowell. The model consists of rectangular atmoshperic slabs. There is a very tiny surrounding sphere placed at the top of the atmoshpere, although this is not used for a unifrom beam source. The watched volume is the entire slab starting at the specified watch_height. The source is a narrow beam with radius = 1cm.

Parameters:

watch_height (float) – Altitude for watched volume in km. The actual volume will be a box starting at the specified value, with a height equal to the spacing specified in the atmospheric model.

spherical_model(watch_alt)

Model consists of concentric spherical shells, enclosed by a large surrounding sphere.

Parameters:

watch_alt (float) – Altitude for watched volume in km. The actual volume will be a shell with the inner radius starting at the specified value, and a height equal to the spacing specified in the atmospheric model.

class cosi_atmosphere.response.Simulate

run_sim(source_file, seed=0, verbosity=1)

Run atmospheric sims.

Parameters:

• source_file (str) – Cosima source file.

• seed (int, optional) – Option to pass seed to cosima (default is no seed).

• verbosity (int, optional) – Output level of cosima (default is 1).

parse_sim_file(sim_file, unique=True)

Parse sim file with StoreSimulationInfo = init-only.

Parameters:

• sim_file (str) – Cosima sim file.

• unique (bool) – A photon may cross the detecting volume numerous times. To only count the first pass, set unique=True (defualt). If false, will count all passes.

parse_sim_file_all_info(sim_file, unique=True)

Parse sim file with StoreSimulationInfo = all. This also returns the sequence of interactions for each measured photon.

Parameters:

• sim_file (str) – Cosima sim file.

• unique (bool) – A photon may cross the detecting volume numerous times. To only count the first pass, set unique=True (defualt). If false, will count all passes.

Note

This method can be used if you are interested in the sequence of interactions for each of the photons. It comes at the cost of a larger sim file, since everything needs to be recorded. The sequence of interactions is stored as a string in the parsed event file, encoded as follows:

S: start, C: Compton event, B: Bremsstrahlung event, A: Photo electric absorption, R: Rayliegh scattering, P: Pair production, X: Pair annihilation, I: Entered watched volume.

In the default mode, we only consider the first time a photon crosses the watched volume. If considering multiple crossings, the photon properties are recorded separately for each crossing. This includes the event sequence, which gives all interactions leading up to the time the photon entered the watched volume.

For some interactions such as PAIR creation, two IA entries are generated, representing the electron and the positron. Thus, the sequence “PP” corresponds to one pair event.

plot_sequence(event_file, elow=10, ehigh=10000, num_ebins=30, show_plots=True)

Plots ditributions for all photon interactions.

Parameters:

• event_file (str) – Name of event file with measured photons. Must be parsed from parse_sim_file_all_info.

• elow (float, optional) – Lower energy bound in keV (defualt is 100 keV).

• ehigh (float, optional) – Upper energy bound in keV (default is 10000 keV.

• num_ebins (int, optional) – Number of energy bins to use (default is 30). Only log binning for now.

• show_plots (bool, optional) – Whether or not to show plots (default is True).

combine_sims(file_list, nsim, output_name)

Combine event files from multiple simulations.

Parameters:

• file_list (list .) – List of file names to combine. The files must be the output from either the parse_sim_file or parse_sim_file_all_info methods.

• nsim (int) – Number of simulated photons for each file. Must be the same for all files.

• output_name (str) – Prefix of output dat file.

class cosi_atmosphere.response.Process(theta, all_events_file='all_thrown_events.dat', measured_events_file='event_list.dat')

Analyze atmophere simulations.

Parameters:

• theta (float) – Off-axis angle of source in degrees.

• all_events_file (str, optional) – Event file with all thrown events (default is output from ParseSims method).

• measured_events_file (str, optional) – Event file with all measured events (default is output from ParseSims method).

dot_product(v1, v2)

Calculates dot product of two vectors, v1 and v2.

Parameters:

• v1 (array) – First vector.

• v2 (array) – Second vector.

Note

Input vector arrays must have shape (rows, cols) = (N,3), where cols is x,y,z coordinates, and rows is number of vectors.

Returns:

dp – Dot product of v1 and v2.

Return type:

float

length(v)

Calculates length of vector.

Parameters:

v (array) – Input vector.

Returns:

Vector length.

Return type:

float

angle(v1, v2)

Calculates angle between two vectors, in radians.

Parameters:

• v1 (array) – First vector.

• v2 (array) – Second vector.

Note

Input vector arrays must have shape (rows, cols) = (N,3), where cols is x,y,z coordinates, and rows is number of vectors.

Returns:

angle – Angle in radians.

Return type:

float

bin_sim(rsp_name, elow=10, ehigh=10000, num_ebins=24, rlow=1e-07, rhigh=1000, num_rbins=120)

Bins the simulations, and writes output files for both starting and measured photons.

Events are weigthed by the ratio of the cosine of incident angle to cosine of measured angle. This is a geometric correction factor to account for the effect of the projected area.

Parameters:

• rsp_name (str) – Name of output response file (use .hdf5 or .h5 extension).

• elow (float, optional) – Lower energy bound in keV (defualt is 10 keV).

• ehigh (float, optional) – Upper energy bound in keV (default is 10000 keV).

• num_ebins (int, optional) – Number of energy bins to use (default is 24). Only log binning for now.

• rlow (float, optional) – Lower radial bound in km (default is 1e-7).

• rhigh (float, optional) – Upper radial bound in km (default is 1000).

• num_rbins (int, optional) – Number of radial bins to use (default is 120).

Note

Response file contains 4 different histograms, each stored as it’s own group. The group names are: 1. starting photons 2. measured_photons_baseline 3. measured_photons 4. rsp_main

get_binning_info()

Get info from histogram.

load_response(name)

Load response files for starting and measured photons.

Parameters:

name (str) – Name of response file.

make_scattering_plots(starting_pos=True, measured_pos=True, spec_i=True, radial_dist=True, theta_prime=True, theta_prime_em=True, rad_em=True, rad_ei=True, show_baseline=True, rsp_file=None)

Visualize the photon scattering.

Parameters:

• starting_pos (bool, optional) – Scatter plot showing starting postions of all photons.

• measured_pos (bool, optional) – Scatter plot showing measured postions of all photons.

• spec_i (bool, optional) – Starting and measured photon spectrum.

• radial_dist (bool, optional) – Radial distribution of measured photons.

• theta_prime (bool, optional) – Incident angle distribution.

• theta_prime_em (bool, optional) – Incident angle versus measured energy.

• rad_em (bool, optional) – Radial distribution versus measured energy.

• rad_ei (bool, optional) – Radial distribution versus initial energy.

• show_baseline (bool, optional) – Option to show angular comparison to un-weighted histogram.

• rsp_file (str, optional) – Option to give name of response file to be loaded (default is None).

calc_tp(show_plot=True)

Calculate the transmission probability.

This gives the probability that a photon of a given energy will pass the atmosphere without being scattered, and thus cross the detecting volume. The transmitted photons do not scatter, and thus there is no corresponding energy dispersion. However, the transmission probability does not account for scattered photons which enter the detecting volume from off-axis angles. These photons will produce a significant energy dispersion.

Parameters:

show_plot (bool, optional) – Option to plot transmission probability and compare to original calculation.

Notes

The transmission probability is a function of energy and off-axis angle of the source.

The transmission probability can also be calculated using the get_total_edisp_matrix method, which is preferred.

plot_tp()

Plot the transmission probability.

get_tp_from_file()

Read original transmission probability from numpy array (for 33 km). This has been used for past COSI studies, including DC1, and it’s included here for a sanity check.

Returns:

tp_energy, tp_array – list of energies and array of TP values, respectively.

Return type:

list, array

get_total_edisp_matrix(show_sanity_checks=False, make_plots=True, tp_file=None)

Get the energy dispersion matrix.

The total energy dispersion is the sum of the beam photons (which don’t scatter) and the scattered photons. Here I calculate all three: beam, scattered, and total. Likewise, I calculate the transmission fraction for all three.

Parameters:

• show_sanity_checks (bool, optional) – Print a comparison of total energy dispersion to summed energy dispersion (beam + scattered), and also for transmission probability, to verify that they are the same (default is False).

• make_plots (bool, optional) – Show plots (default is True).

• tp_file (str, optional) – Option to overlay TP from analytical calculation in plot. Specify file name from TPCalc class. Default is None.

make_edisp_matrix(idm_watch, em_watch, ia_watch)

Make energy dispersion matrix.

Parameters:

• idm_watch (ArrayLike) – IDs of measured photons for watched region.

• em_watch (ArrayLike) – Energy histogram of photons for watched region.

• ia_watch (ArrayLike) – Angle histogram of photons for watched region.

plot_edisp_matrix(matrix, savefile)

Plot energy dispersion matrix.

plot_tp_from_edisp(tp_file=None)

Plot the transmission probability.

Parameters:

tp_file (str, optional) – Option to overlay TP from analytical calculation in plot. Specify file name from TPCalc class. Default is None.

PL(E, p)

General PL spectral model, normalized to 1e-3 at 100 keV.

Parameters:

• E (float) – Energy in keV.

• p (float) – Power of dN/dE.

Returns:

dn/de – Value at given energy.

Return type:

float

PL_interp(p)

Returns interpolated PL function for energies between 100 keV to 10 MeV.

Paramters

pfloat

Power of dN/dE

returns:

pl_interp – Interpolated function.

rtype:

scipy:interpolate:interp1d

ff_correction(model_flux, name, show_plots=True)

Calculate atmospheric correction factor by forward folding the atmospheric energy dispersion with the model counts.

The energy dispersion matrices must first be generated via the ‘get_total_edisp_matrix’ method.

Parameters:

• model_flux (scipy:interpolate:interp1d) – Interp1d object giving the model flux as a function of energy.

• name (str) – Name of saved files.

• show_plots (bool) – Option to plot correction factor and ratio.

class cosi_atmosphere.response.ProcessSpherical(r_alt, all_events_file='all_thrown_events.dat', measured_events_file='event_list.dat')

Analyze atmophere simulations from spherical mass model.

Parameters:

• r_alt (float) – Altitude of observations in km.

• all_events_file (str, optional) – Event file with all thrown events (default is output from ParseSims method).

• measured_events_file (str, optional) – Event file with all measured events (default is output from ParseSims method).

dot_product(v1, v2)

Calculates dot product of two vectors, v1 and v2.

Parameters:

• v1 (array) – First vector.

• v2 (array) – Second vector.

Note

Input vector arrays must have shape (rows, cols) = (N,3), where cols is x,y,z coordinates, and rows is number of vectors.

Returns:

dp – Dot product of v1 and v2.

Return type:

float

length(v)

Calculates length of vector.

Parameters:

v (array) – Input vector.

Returns:

Vector length.

Return type:

float

angle(v1, v2)

Calculates angle between two vectors, in radians.

Parameters:

• v1 (array) – First vector.

• v2 (array) – Second vector.

Note

Input vector arrays must have shape (rows, cols) = (N,3), where cols is x,y,z coordinates, and rows is number of vectors.

Returns:

angle – Angle in radians.

Return type:

float

line_point(o, u, d)

Finds point on a line segment.

Parameters:

• o (array-like) – Point on line. In our case it’s the starting photon position.

• u (array-like) – Direction of line segment. In our case it’s the starting photon direction.

• d (array-like) – Distance from specified point (o) in the specified direction (u).

Returns:

x – Point on line.

Return type:

array-like

find_projection(o, u, r)

Finds projection of line segment onto sphere.

Everything should be done in Cartesian coordinates. The closest point of intersection defines the normal vector to the surface, relative to the initial direction. The true incident angle (before scattering) is the angle between these two vectors.

Parameters:

• o (array-like) – Point on line. In our case it’s the starting photon position.

• u (array_like) – Direction of line segment. In our case it’s the starting photon direction.

• r (float) – Radius of watched sphere. Default is r_earth + 33.1 km.

Returns:

intersection – Point of intersection

Return type:

array

Note

All inputs must have self-consistent units.

bin_sim(name, elow=10, ehigh=10000, num_ebins=17, anglow=0, anghigh=100, num_angbins=26, nside=16, scheme='ring', weighted=True)

Bins the simulations, and writes output files for both starting and measured photons.

Events are weigthed by the ratio of the cosine of incident angle to cosine of measured angle. This is a geometric correction factor to account for the effect of the projected area.

Parameters:

• name (str) – Name of output response file (use .hdf5 or .h5 extension).

• elow (float, optional) – Lower energy bound in keV (defualt is 100 keV).

• ehigh (float, optional) – Upper energy bound in keV (default is 10000 keV.

• num_ebins (int, optional) – Number of energy bins to use (default is 17). Only log binning for now.

• anglow (float, optional) – Lower bound of angle in degrees (default is 0).

• anghigh (float, optional) – Upper bound of anlge in degrees (default is 100)

• num_angbins (int, optional) – Number of angular bins (default is 26).

• nside (int, optional) – nside for healpix binning (defualt is 16).

• scheme (str, optional) – Scheme for healpix binning (ring or nested). Default is ring.

• weighted (bool, optional) – Wether or nor to use a weighted histogram (default is True).

Note

Response file contains 5 different histograms, each stored as its own group. The group names are: 1. primary_rsp 2. starting_photons 3. starting_pe, ons_rsp 4. measured_photons 5. measured_photons_rsp

get_binning_info()

Get info from histogram.

make_scattering_plots(pos_init=True, pos_meas=True, pos_proj=True, pos_ang_ri=True, pos_ang_di=True, ri=True, pos_ang_rm=True, pos_ang_dm=True, rm=True, theta_dist=True, ang_dist=None, spec=True, rsp_file=None)

Visualize the simulated photons.

Parameters:

• pos_init (bool, optional) – 3d scatter plot of initial postions

• pos_meas (bool, optional) – 3d scatter plot of measured positions

• pos_proj (bool, optional) – 2d projection onto xy-axis of initial positions

• pos_ang_ri (bool, optional) – healpix map of initial positions

• pos_ang_di (bool, optional) – healpix map of initial directions

• ri (bool, optional) – initial radius (relative to both Earth center and surrounding sphere disk)

• pos_ang_rm (bool, optional) – healpix map of measured positions

• pos_ang_dm (bool, optional) – healpix map of measured directions

• rm (bool, optional) – measured radius with respect to Earth center

• theta_dist (bool, optional) – distributions of incident angle (initial and measured)

• ang_dist (list, optional) – Distribution of measured anlges for a given initial angle. Takes list with incident angle in degrees and energy in keV. Default is None.

• spec (bool, optional) – Spectrum of initial and measured photons (dN/dE).

• rsp_file (str, optional) – Name of response file to use.

load_response(name)

Load response file and corresponding normalization.

Parameters:

• name – Name of response file.

• str – Name of response file.

get_total_edisp_matrix(theta, make_plots=True, rsp_file=None, tp_file=None)

Get the energy dispersion matrix.

The total energy dispersion is the sum of the transmitted photons (which don’t scatter) and the scattered photons. Here I calculate all three: transmitted, scattered, and total. Likewise, I calculate the transmission probability for all three.

Parameters:

• theta (float) – Incident angle of source in degrees.

• make_plots (bool, optional) – Show plots (default is True).

• rsp_file (str, optional) – Name of response file to load.

• tp_file (str, optional) – Option to overlay TP from analytical calculation in plot. Specify file name from TPCalc class. Default is None.

get_theta_bin(theta, rsp_file=None)

Returns index of the theta bin.

Parameters:

• theta (float) – Incident angle in degrees.

• rsp_file (str, optional) – Prefix of response file to load.

Returns:

theta_bin – Index of theta bin.

Return type:

int

get_energy_bin(energy, rsp_file=None)

Returns index of the energy bin.

Parameters:

• energy (float) – Energy in keV.

• rsp_file (str, optional) – Prefix of response file to load.

Returns:

energy_index – Index of energy bin.

Return type:

int

normalize_response()

Normalizes the atmospheric response matrix by the total photons thrown in E_i and theta_i.

make_edisp_matrix(theta_bin, comp='total')

Make energy dispersion matrix.

Parameters:

• theta_bin (int) – index of incident angle.

• comp (str, optional) – Component to use for constructing matrix. Either transmitted, scattered, or total (default is total).

class cosi_atmosphere.response.ProcessSphericalRep2(r_alt, all_events_file='all_thrown_events.dat', measured_events_file='event_list.dat')

bin_sim(name, elow=10, ehigh=10000, num_ebins=17, anglow=0, anghigh=100, num_angbins=26, dtheta_low=0, dtheta_high=180, num_dtheta_bins=45)

Bins the simulations, and writes output files for both starting and measured photons.

Events are weigthed by the ratio of the cosine of incident angle to cosine of measured angle. This is a geometric correction factor to account for the effect of the projected area.

Parameters:

• name (str) – Name of output response file (use .hdf5 or .h5 extension).

• elow (float, optional) – Lower energy bound in keV (defualt is 100 keV).

• ehigh (float, optional) – Upper energy bound in keV (default is 10000 keV.

• num_ebins (int, optional) – Number of energy bins to use (default is 17). Only log binning for now.

• anglow (float, optional) – Lower bound of angle in degrees (default is 0).

• anghigh (float, optional) – Upper bound of anlge in degrees (default is 100)

• num_angbins (int, optional) – Number of angular bins (default is 26).

• dtheta_low (float, optional) – Lower bound of delta theta in degrees (default is 0).

• dtheta_high (float, optional) – Upper bound of delta theta in degrees (default is 180).

• num_dtheta_bins (int, optional) – Number of delta theta bins (default is 90).

Note

Response file contains 2 different histograms, each stored as its own group. The group names are: 1. starting_photons_rsp 2. measured_photons

get_binning_info()

Get info from histogram.

load_response(name)

Load response file and corresponding normalization.

Parameters:

• name – Name of response file.

• str – Name of response file.

make_scattering_plots(angle, energy, name, sig_y=0.1, sig_x=1.5, scale=100, rsp_file=None, dist_fig_kwargs={}, rainbow_fig_kwargs={})

Visualize the response.

Parameters:

• angle (float) – Off-axis angle of source in degrees.

• energy (float) – Energy band in keV.

• name (str) – Prefix of saved plots (do not include extension).

• sig_y (float, optional) – Standard deviation of Gaussian kernal in y-direction used for smoothing image (default is 0.1).

• sig_x (float, optional) – Standard deviation of Gaussian kernal in y-direction used for smoothing image (default is 0.1).

• scale (float, optional) – Factor of array max used for plotting contour. Default is 100.

• rsp_file (str, optional) – Name of response file to use.

• dist_fig_kwargs (dict, optional:) – Pass any kwargs to plt.gca().set(), for distribution plot.

• rainbow_fig_kwargs (dict, optional:) – Pass any kwargs to plt.gca().set(), for rainbow plot.

get_dtheta_bin(dtheta, rsp_file=None)

Returns index of the delta theta bin.

Parameters:

• dtheta (float) – delta theta in degrees.

• rsp_file (str, optional) – Prefix of response file to load.

Returns:

theta_bin – Index of theta bin.

Return type:

int

get_total_edisp_matrix(theta, dtheta_max, make_plots=True, rsp_file=None, tp_file=None)

Get the energy dispersion matrix.

The total energy dispersion is the sum of the transmitted photons (which don’t scatter) and the scattered photons. Here I calculate all three: transmitted, scattered, and total. Likewise, I calculate the transmission probability for all three.

Parameters:

• theta (float) – Incident angle of source in degrees.

• dtheta_max (float) – Max delta theta to use in degrees.

• make_plots (bool, optional) – Show plots (default is True).

• rsp_file (str, optional) – Name of response file to load.

• tp_file (str, optional) – Option to overlay TP from analytical calculation in plot. Specify file name from TPCalc class. Default is None.

make_edisp_matrix(theta_bin, dtheta_max, comp='total')

Make energy dispersion matrix.

Parameters:

• theta_bin (int) – index of incident angle.

• dtheta_max (int) – Bin of max delta theta to use.

• comp (str, optional) – Component to use for constructing matrix. Either transmitted, scattered, or total (default is total).

normalize_response()

Normalizes the atmospheric response matrix by the total photons thrown in E_i and theta_i.

class cosi_atmosphere.response.TPCalculator

Calculates transmission probability for gamma rays.

Mass attenuation coeffients are from NIST XCOM: https://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html

Atmospheric data is fom NRLMSIS: https://swx-trec.com/msis

A great review of the relationship between transmitted intensity, linear attenuation coefficient, and mass attenuation coefficient can be found here: https://www.nde-ed.org/Physics/X-Ray/attenuationCoef.xhtml

read_atm_model(atm_file)

Reads in atmospheric profile and makes interpolated function.

Parameters:

atm_file (str) – Name of atmosphere file.

Returns:

density – Interpolted function of altitude (in cm) vs. density (in g/cm3).

Return type:

scipy:interpolate:interp1d

Note

Atmosphere file is output from Atmosphere class.

relative_abundance()

Calculates relative abundance of each element of the atmosphere as a function of altitude.

read_mass_attenuation(nist_data='default', show_plot=True)

Reads in attenuation data from NIST XCOM.

NIST XCOM data can be downloaded from: https://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html

Data file must include all column options from XCOM, with column names as follows: energy[MeV] scat_coh[cm2/g] scat_incoh[cm2/g] photo_abs[cm2/g] pair_prod_NF[cm2/g] pair_prod_EF[cm2/g] atten_tot_wCS[cm2/g] atten_tot_woCS[cm2/g]

Parameters:

• nist_data (str, optional) – Name of NIST data file. Default file is calculated assuming an atmospheric composition of 78% N2 and 22% O2.

• show_plot (bool) – Option to plot attenuation data (defualt is True).

plot_attenuation()

Plots attenuation data from NIST XCOM.

coslaw(A, theta, x)

Uses law of cosines to determine the altitude for any distance along the line of sight, relative to Earth’s surface, for a given altitude of observation and off-axis angle of source.

Parameters:

• A (float) – Altitude of observations in cm.

• theta (float) – Off-axis angle of source in degrees.

• x (float) – distance along line-of-sight in cm.

Returns:

r – Side length of triangle in cm (opposite the angle pi - theta).

Return type:

float

los_integral(A, theta, dx=10000.0, alt_max=20000000.0)

Estimates line of sight integral with respect to radius of spherical reference frame.

Parameters:

• A (float) – Altitude of observations in cm.

• theta (float) – Off-axis angle of source in degrees.

• dx (float, optional) – x step size in cm (default is 1e4).

• alt_max (float, optional) – Top of atmosphere in cm (default is 200e5 cm).

Returns:

integral – Integral along line of sight, which as units of g/cm2

Return type:

float

calc_tp(A, theta, output_name=None, dx=10000.0, alt_max=200, show_plot=True)

Calculates transmission probability.

Parameters:

• A (float) – Altitude of observations in km.

• theta (float) – Off-axis angle of source in degrees.

• output_name (str, optional) – Name of output file (default is None, in which case no output file is written).

• dx (float, optional) – x step size in cm (default is 1e4, i.e. 100 m).

• alt_max (float, optional) – Top of atmosphere in km (default is 200).

• show_plot (bool, optional) – Option to plot transmission probability (default is True).

Returns:

tp – Transmission probability as a function of altitude.

Return type:

array

plot_tp()

Plots transmission probability.

cosima_tp_file(A, output_name, theta_list='default', energy_list='default', show_plot=True)

Makes input transmission probability file for MEGAlib’s cosima.

Parameters:

• A (float) – Altitude of observations in km.

• output_name (str) – Name output files. Output will be saved as both .npy file as well as .dat file for cosima input.

• theta_list (list, optional) – List of fff-axis angles to include in degrees. Default is 0-90 degrees with 5 degree spacing. TP is set to zero beyond 90 degrees.

• energy_list (list, optional) – List of energy values to use in keV. Default is 50 keV - 10 MeV, using log spacing.

• show_plot (bool, optional) – Option to plot 2D TP array.