-
Notifications
You must be signed in to change notification settings - Fork 0
/
remnant_Record_thermal_spectra_parse.py
155 lines (92 loc) · 5.26 KB
/
remnant_Record_thermal_spectra_parse.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
# coding: utf-8
# Create thermal spectra for the FS and the RS. Last step (thought to be compared with real spectra of SNRs)
import numpy as np
from astropy.table import Table, vstack, Column
import os
from scipy import interpolate
import pyatomdb as pyat
#####################################################################################################################################
# IMPORTANT: READ THIS
# Definition of filename as an optional argument to run the script multiple times in
# the terminal with all the different values via a simple loop in a shell script:
# python crspectra_resp.py filename
# where filename is any subdirectory in the present location
# See https://docs.python.org/3/library/argparse.html
import argparse
pathw = os.getcwd() + '/'
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument(
'index', metavar='str', type=str, choices=[x[1] for x in os.walk(pathw)][0],
nargs='+', help='introduce any subdirectory in the present location')
args = parser.parse_args()
filename = args.index[0] # Subdirectory selected to create the synthetic spectra with PyAtomDB.
# e.g. DDTa.0.10, DDTa.0.50, DDTa.1.00, DDTa.5.00
# Note: usage of os.walk(pathw)
#for (path, dirs, files) in os.walk(path):
# print path
# print dirs
# print files
# So, if the script is to be run from the directory itself (not the upper level suggested at the beginning of the script),
# instead of [x[2] for x in os.walk(pathw)], just comment put this parsing section and redefine path_files:
#path_files = pathw
#####################################################################################################################################
path_files = pathw + filename + '/'
age_sk = np.genfromtxt(path_files + 'data_sk_vs_rad.dat', skip_footer=2)[:,2] # age recorded for each shock
explmodel = path_files.split(os.sep)[-2][6:]
ISMrho = path_files.split(os.sep)[-3][-3:]
t_RS = np.zeros((2), dtype='object')
t_FS = np.zeros((2), dtype='object')
# IMPORTANT: READ THIS
# PyAtomDB calculates emissivities in [ct cm^3 s^-1 bin^-1]
# for given energy bins, e.g., E = [1, 2, 3, 4, 5...] eV,
# and then retrieves a photon spectrum whose length is
# len(ebins) -1, corresponding to the mean energy for
# each bin, in this case E = [1.5, 2.5, 3.5, 4.5...] eV
# (see dummyfirst flag in make_ion_spectrum)
# By default, crhydronei calculates thermal RS/FS photon
# spectra in 10000 energy bins with a step of 1.2 eV, from
# 0.95 keV to 12.0938 keV. To summarize,
# Input: 10000 energy bins
# Output: 9999 mean energy values, 9999 spectral values
En_start = 9.50E-2
En_step = 1.20E-3
En_nelem = 10000
# Input energies
En_th_bins = np.asarray([En_start + En_step*x for x in range(En_nelem)]) # Energy bin limits from thermal spectra
# Ouput energies that will be saved in the final files
En_th = np.asarray([En_start + En_step/2. + En_step*x for x in range(En_nelem - 1)]) # Average values for each bin ("Ebin")
ebins = En_th_bins
Ebin = En_th
# These have the information for the last shock
thermal_photons_xspec_unabs_RS = np.genfromtxt(path_files + 'thermal_photons_xspec_unabs_RS.dat')[:,2]
thermal_photons_xspec_unabs_FS = np.genfromtxt(path_files + 'thermal_photons_xspec_unabs_FS.dat')[:,2]
# The units for these thermal spectra are [ph cm^-2 s^-1 bin^-1]
# So, given that the energy step is 1.2 eV,
# [ph cm^-2 s^-1 bin^-1] * 1 bin / (1.2E-3 keV) = [ph cm^-2 s^-1 keV^-1]
##############################################################################
# thermal.f90
#!! file= 'thermal_photons_xspec_unabs_FS(RS)_evo.dat'
# write(i_write_unit_2,"(A7, I5, 1p1e12.3, 2I5)") "#Epoch ", &
# i_sk_ct, time_yr_run, &
# i_skip_Doppler_1, &
# i_skip_Doppler_2
#!! file= 'thermal_photons_xspec_unabs_FS(RS)_evo.dat'
# write(i_write_unit_2,"(1p3e15.5)") &
# binlo, binhi, &
# totfluxk/bin_kev ![N/cm^2/s] per E bin
# enddo
#
# write(*,*) "Finished summing for shell", i_sk_ct, i_FS_RS
##############################################################################
t_RS[0] = Column(Ebin, name='Energy', unit='keV')
t_FS[0] = Column(Ebin, name='Energy', unit='keV')
t_RS[1] = Column(thermal_photons_xspec_unabs_RS / En_step, name = str(int(age_sk[-1])) + 'yr', unit = 'ph+1cm-2s-1keV-1')
t_FS[1] = Column(thermal_photons_xspec_unabs_FS / En_step, name = str(int(age_sk[-1])) + 'yr', unit = 'ph+1cm-2s-1keV-1')
# See https://astropy.readthedocs.io/en/v0.1/wcs/units.html for all the units (including counts and photons)
#####################################################################################################################################
# Finally, save a FITS table with the energy bins in the first column and the different photon spectra (ages) in the remaining ones
Table_end_RS = Table([x for x in t_RS])
Table_end_RS.write('spectra_ph_thermal_RS_' + explmodel + '_' + ISMrho + '_' + str(int(age_sk[-1])) + 'yr' + '.fits', overwrite='True')
Table_end_FS = Table([x for x in t_FS])
Table_end_FS.write('spectra_ph_thermal_FS_' + explmodel + '_' + ISMrho + '_' + str(int(age_sk[-1])) + 'yr' + '.fits', overwrite='True')