calculus_habing.py
2.67 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
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 14 18:29:52 2022
@author: frede
"""
import numpy as np
from astropy.table import Table
import matplotlib.pyplot as plt
''' 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
''' calculus of wavelength '''
wave_intensity = Table.read('./converted/new_habing1968.txt', format = 'ascii')
wavelength = wave_intensity['col1']
density_range = (wavelength >= h/ev * c *1e9/13.6) & (wavelength <= h/ev * c * 1e9/6)
habing_RF = np.trapz( ( (-25/6) * wavelength[density_range]**2/1e6 +\
(25/2) * wavelength[density_range]/1e4 -\
(13/300) ) * 1e-14, wavelength[density_range])
frac = wavelength/100
wavelength_intensity = []
Draine_list = []
Mathis_list = []
for i in frac:
# test = ( (c*1e2)/(2*np.pi) ) * ( (-25/6) * wavelength[i]**2/1e6 +\
# (25/2) * wavelength[i]/1e4 -\
# (13/300) ) * 1e-14
# wavelength_intensity.append(test)
if i >= 0.912:
Draine = 4e-14 * 1.71 * (1/i**5) * (31.016*i**2 - 49.913*i + 19.897)
Draine_list.append(Draine)
for i in wavelength:
if i == 91.2 :
Mathis_list.append(0)
if i > 91.2 and i <= 110:
Mathis = 1.287e-9 * (i*1e-3)**(4.4172)
Mathis_list.append(Mathis)
if i > 110 and i <= 134:
Mathis = 6.825e-13 * (i*1e-3)
Mathis_list.append(Mathis)
if i > 134 and i <= 246:
Mathis = 2.373e-14 * (i*1e-3)**(-0.6678)
Mathis_list.append(Mathis)
Habing = [6.6e-14, 7e-14, 4e-14]
plt.figure(figsize=(10.90,8.93))
plt.axis()
plt.semilogy((h/ev *c*1e9)/wavelength, Draine_list,ls='dashdot', c='k',label='Draine (1978)')
plt.semilogy((h/ev *c*1e9)/wavelength[0:len(Mathis_list)], Mathis_list,c='k',label='MMP83')
plt.semilogy([(h/ev *c*1e9)/220, (h/ev *c*1e9)/140, (h/ev *c*1e9)/100], Habing,ls='--',c='k',marker="s",label='Habing (1968)')
plt.ylim([1e-14,2e-13])
plt.xlim([5,14])
plt.xlabel('Energy (eV)')
plt.xticks([6, 8, 10, 12, 13.6])
plt.ylabel('νu$_ν$ (erg cm$^{-3}$)')
plt.text(5.8, 1.6e-14, 'Local Interstellar UV', size=25)
plt.legend()
# pl = plt.plot([1, 4, 3, 7, 8], color = 'blue')
# plt.legend(pl, ['blue'])
# ax2 = plt.gca().twinx()
# pl2 = ax2.plot([80, 70, 30, 40, 10], color = 'red')
# plt.legend(pl2, ['red'])
# plt.show()
# document = Table([wavelength , wavelength_intensity] , names =('col1','col2'))
# document.meta['comments'] = ['wavelength (nm)' , 'intensity (erg s-1 cm-2 nm-1 sr-1)']
# document.write('habing1968.txt', format = 'ascii', overwrite=True)
print(habing_RF)