Analytical Calculation of Transmission Probability (TP)

Here we will show how to calculate the gamma-ray transmission probability analytically.

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

[1]:

# Imports:
from cosi_atmosphere.response.AtmosphericProfile import Atmosphere
from cosi_atmosphere.response.TPCalc import TPCalculator
import numpy as np
import os

Get atmosphere model

We’ll use values based on the 2016 COSI balloon flight:

[2]:

instance = Atmosphere()
date = np.array(['2016-06-13 12:00:00'], dtype="datetime64[h]")
lat = -5.66
lon = -107.38
alts = np.linspace(0, 200, 2001) # km; spacing is 0.1 km (100 m)
atm_model = instance.get_atm_profile("rep_atm_model.dat",date,lon,lat,alts)

Get mass attenuation coefficients

Define instance of TP class, and read in atmosphere and mass attenuation coefficient data. Here we will use the default mass attenaution data, calculated for an atomospheric composition with 78% N2 and 22% O2. Note that you can also pass a different attenuation data file from XCOM if needed.

[4]:

tp_instance = TPCalculator()
tp_instance.read_atm_model("rep_atm_model.dat")
tp_instance.read_mass_attenuation()

../../_images/tutorials_transmission_probability_analytical_TP_analytical_5_0.png

Let’s take a look at the atmospheric relative abundance as a function of altitude. As can be seen, the atmosphere is mostly N2 (~78%) and O2 (~22%), although the composition does begin to change above ~100 km.

[5]:

tp_instance.relative_abundance()

../../_images/tutorials_transmission_probability_analytical_TP_analytical_7_0.png

Calculate transmission probability

Now we can calculate the TP. Let’s use an altitude of 33.5 km, and an off-axis angle of 50 degrees:

[6]:

tp_instance.calc_tp(33.5,50,output_name="TP_33p6km_theta50_analytical.dat")

../../_images/tutorials_transmission_probability_analytical_TP_analytical_9_0.png

[6]:

0     0.000000
1     0.000000
2     0.000000
3     0.000000
4     0.000000
        ...
75    0.788183
76    0.788092
77    0.788092
78    0.788001
79    0.787911
Name: atten_tot_wCS[cm2/g], Length: 80, dtype: float64

Make cosima input file

Now let’s calculate the TP file for MEGAlib’s cosima:

[7]:

tp_instance.cosima_tp_file(33.5,"TP_cosima_33p5km")

Max theta for given altitude due to Earth occultation [deg]: 95.85962050942686

../../_images/tutorials_transmission_probability_analytical_TP_analytical_11_1.png

Let’s take a look at the first 50 lines of the TP file for cosima:

[8]:

%%capture
os.system("head -n 50 TP_cosima_33p5km.dat")

IP LIN

# Theta axis in degrees:
XA 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 90.001 180
# Energy axis in keV:
YA 49.99999999999999 55.709569209792114 62.07112203081239 69.1591093740999 77.05648380328007 85.85567035004112 95.65964818837426 106.58315582669209 118.75403392250328 132.3147214349539 147.42392262509864 164.2584644132395 183.01536583047076 203.91414378376078 227.19938211953857 253.14359405220867 282.05041145734094 314.25813835466175 350.14371016818893 390.1271050997625 434.6762592434215 484.31254296349834 539.6168626279089 601.2364591068078 669.8924826012209 746.388432451844 831.6195607012854 926.5833494621031 1032.3911847100012 1150.2813631236224 1281.6335841933885 1427.9850972043125 1591.0486920651067 1772.732744535005 1975.1635504427547 2200.7102102809863 2452.012355406897 2732.0110403361195 3043.983162670423 3391.5798134845986 3778.8690069970835 4210.383289600847 4691.172785435414 5226.864299306152 5823.727168647751 6488.746635214575 7229.705595185846 8055.27568442854 8975.118764912551 10000.0

AP 0 0 0.22430069728542043
AP 0 1 0.24218057678547583
AP 0 2 0.25960463158194774
AP 0 3 0.27325375634265486
AP 0 4 0.28846149355284895
AP 0 5 0.3017370573984878
AP 0 6 0.31447377600552795
AP 0 7 0.32610394341598814
AP 0 8 0.33718088513948924
AP 0 9 0.3495227181710535
AP 0 10 0.3632738822000062
AP 0 11 0.37542014840684085
AP 0 12 0.3883142502072704
AP 0 13 0.4020438725318201
AP 0 14 0.4142621321167309
AP 0 15 0.4278756116761066
AP 0 16 0.44304363331008945
AP 0 17 0.4580620798823969
AP 0 18 0.47215621882114994
AP 0 19 0.48785978699436794
AP 0 20 0.502685060291009
AP 0 21 0.5183557741400744
AP 0 22 0.5336573269215051
AP 0 23 0.5496880294077212
AP 0 24 0.5639701218800571
AP 0 25 0.5798831062610221
AP 0 26 0.5961380415999531
AP 0 27 0.6114623496101949
AP 0 28 0.6276857900351837
AP 0 29 0.6425926310790486
AP 0 30 0.6584207962147799
AP 0 31 0.6733136402175969
AP 0 32 0.6875361920200677
AP 0 33 0.7012934126679449
AP 0 34 0.716621589384291
AP 0 35 0.7287216900085306
AP 0 36 0.7405662701411874
AP 0 37 0.7537633992744107
AP 0 38 0.767679914763288
AP 0 39 0.7778386695893907
AP 0 40 0.7891574666907925
AP 0 41 0.7997597764831746
AP 0 42 0.8092200203060063

Altitude variation

For fun, let’s make a short gif showing how the TP changes with altitude. First, calculate the TP for a range of altitudes:

[9]:

alt_list = np.arange(25,46,1).tolist()
image_list = []
for each in alt_list:
    print(each)
    this_name = "TP_%s" %str(each)
    image_list.append(this_name + ".png")
    tp_instance.cosima_tp_file(each,this_name,show_plot=False)

25
Max theta for given altitude due to Earth occultation [deg]: 95.06474157390615
26
Max theta for given altitude due to Earth occultation [deg]: 95.16470698248163
27
Max theta for given altitude due to Earth occultation [deg]: 95.2627486169415
28
Max theta for given altitude due to Earth occultation [deg]: 95.35897207377991
29
Max theta for given altitude due to Earth occultation [deg]: 95.45347360574536
30
Max theta for given altitude due to Earth occultation [deg]: 95.54634124035411
31
Max theta for given altitude due to Earth occultation [deg]: 95.6376557319047
32
Max theta for given altitude due to Earth occultation [deg]: 95.72749137632073
33
Max theta for given altitude due to Earth occultation [deg]: 95.81591671223384
34
Max theta for given altitude due to Earth occultation [deg]: 95.90299512713601
35
Max theta for given altitude due to Earth occultation [deg]: 95.98878538385713
36
Max theta for given altitude due to Earth occultation [deg]: 96.07334207980502
37
Max theta for given altitude due to Earth occultation [deg]: 96.15671604917632
38
Max theta for given altitude due to Earth occultation [deg]: 96.23895471656216
39
Max theta for given altitude due to Earth occultation [deg]: 96.3201024089411
40
Max theta for given altitude due to Earth occultation [deg]: 96.40020063189331
41
Max theta for given altitude due to Earth occultation [deg]: 96.47928831492634
42
Max theta for given altitude due to Earth occultation [deg]: 96.55740203003307
43
Max theta for given altitude due to Earth occultation [deg]: 96.63457618696705
44
Max theta for given altitude due to Earth occultation [deg]: 96.71084320819736
45
Max theta for given altitude due to Earth occultation [deg]: 96.78623368606944

I’ll make a gif following this nice example: https://www.codedex.io/projects/create-a-gif-with-python. Make sure you have the package installed.

[10]:

import imageio.v3 as iio

filenames = image_list
images = [ ]

for filename in filenames:
  images.append(iio.imread(filename))

iio.imwrite('TP_alt_variation.gif', images, duration = 500, loop = 0)

Now play the gif:

MyGif