radiation_fields.py
6.28 KB
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
# -*- coding: utf-8 -*-
"""
Created on Thu Jun 23 09:55:18 2022
@author: frede
"""
import numpy as np
from astropy.table import Table
'''
This program calculates the scaling factor of the radiation field G0 relative
to the data of an object or a set of objects radiating for given distances
A file containing my G0 values with the associated distance will then be created
We take the standard Habing value for the interstellar medium 5.6x10^-14 ergs cm-3
between 91.2 and 240nm as normalization factor of G0
The file to be used must contain in column 1 the wavelength in nm,
and in column 2 the intensity per wavelength per sr (in erg cm-2 s-1 nm-1 sr-1).
Each column should not have a name
Example of execution of the program for the star of 10 solar radius HD200775,
for a distance to the star at 20pc :
radiation_field(filename='HD200775_RF.txt', star_radius=10, parsec=20)
If I want to generate a hundred values of G0 :
radiation_field('HD200775_RF.txt', 10, 20, ISRF=False, RF_list=True)
The program will tell me that for 20pc my first value of G0 is not of the order of 1e6,
and will make me enter a new value of distance (in pc) until finding G0 of the order of 1e6
If I want to calculate G0 for the interstellar medium
with an approach described by Habing (1968),
(using Habing file with datas of wavelength and intensity per wavelength per sr) :
radiation_field('habing1968.txt', 1, 1, ISRF=True)
For the interstellar medium, the parameters star_radius, parsec and RF_list
are no longer important;
the part of the program executed will be independent of them
'''
''' Constants '''
h = 6.62607015e-34 #Planck constant in J s
c = 299792458 #Light speed in m s-1
''' Conversions '''
ev = 1.602176634e-19 # 1 ev = 1.602176634e-19 J and value of electron charge
def radiation_field(filename, star_radius , parsec, ISRF=False, RF_list=False):
"""
----------
filename : str,
name of the .txt containing intensity and wavelength data
parsec : float,
distance in parsec from the star
star_radius : float,
star radius, in unit of solar radius
ISRF : bool,
Interstellar Radiation Field; if true, then the filename is a file for ISRF
if false, the radiation field of a star is studied
RF_list: bool,
If false, the calculation of G0 is done for a given distance,
if true, it is done for a hundred distances whose input distance
will correspond to a G0 of the order of 1e6
returns the value of G0 for a given distance, the distance,
wavelengths from 91.2nm (13.6ev) to 2 microns, the corresponding intensity
at a distance d_0, the energy intensity, the energy, the values of ISRF and RF_list
-------
"""
''' association of file data '''
wave_intensity = Table.read(filename, format='ascii')
wavelength = wave_intensity['col1'] #in nm
wavelength_intensity = wave_intensity['col2'] #in erg cm-2 s-1 nm-1 sr-1
#=> intensity per nm per sr
''' useful parameters '''
radiation_field = []
distance = []
energy = ( ( (h/ev) * (c*1e2) ) / (wavelength*1e-7) )[::-1] #in ev
r = star_radius*7e10 #radius of the star in cm
d_0 = 3.086e18*parsec #1pc = 3.086e18cm
''' energy intensity in erg s-1 cm-2 sr-1 eV-1 '''
energy_intensity = wavelength_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)
''' indentation '''
i=0
if ISRF :
''' Radiation Field '''
intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
#far UV range, 91.2nm corresponding to 13.6ev
g_0 = np.trapz( 2*np.pi * wavelength_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c*1e2)
#g_0 dimensionless
radiation_field.append(g_0)
distance.append(0)
else:
if RF_list :
while i < 100:
''' distance '''
if i == 0 :
d = d_0
else:
d = 5*i * d_0
''' initialization '''
g_0 = 0
''' geometrical dillution '''
diluted_intensity = wavelength_intensity*(r/d)**2
#intensity diluted and extincted in erg cm-2 s-1 nm-1 sr-1
'''energy intensity in erg s-1 cm-2 sr-1 eV-1'''
energy_intensity = diluted_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)
''' Radiation Field '''
intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
g_0 = np.trapz(2 * np.pi * diluted_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c * 1e2)
distance.append(d)
radiation_field.append(g_0)
if (radiation_field[0] >= 2e6 or radiation_field[0] < 1e6):
radiation_field = []
parsec = float(input('radiation_field[0] = {}, set a value of parsec such that radiation_field[0] is proportional to 1e6 : '.format(g_0),))
distance = []
d_0 = 3.086e18*parsec
i = 0
else:
i = i + 1
else:
d = d_0
''' geometrical dillution '''
diluted_intensity = wavelength_intensity*(r/d)**2
'''energy intensity in erg s-1 cm-2 sr-1 eV-1'''
energy_intensity = diluted_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)
''' Radiation Field '''
intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
g_0 = np.trapz(2 * np.pi * diluted_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c * 1e2)
distance.append(d)
radiation_field.append(g_0)
''' creation of data files '''
document_1 = Table([wavelength , wavelength_intensity] , names =('col1','col2'))
document_1.meta['comments'] = ['wavelength (nm)' , 'intensity (erg s-1 cm-2 nm-1 sr-1)']
document_1.write('intensity_per_wavelength.txt', format = 'ascii', overwrite=True)
document_2 = Table([distance , radiation_field], names =('col1', 'col2'))
document_2.meta['comments'] = ['distance (cm)' , 'radiation field (dimensionless)']
document_2.write('radiation_field_per_distance.txt', format = 'ascii', overwrite=True)
return radiation_field[0], distance[0], wavelength, wavelength_intensity, energy_intensity, energy, ISRF, RF_list