forked from eblur/dust
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcmindex.py
136 lines (113 loc) · 4.59 KB
/
cmindex.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
import os.path
import numpy as np
import constants as c
from scipy.interpolate import interp1d
#------------- Index of Refraction object comes in handy --
#class CM(object): # Complex index of refraction
# def __init__( self, rp=1.0, ip=0.0 ):
# self.rp = rp # real part
# self.ip = ip # imaginary part
#-------------- Complex index of refraction calcs ---------
# ALL CM OBJECTS CONTAIN
# cmtype : string ('Drude', 'Graphite', or 'Silicate')
# rp : either a function or scipy.interp1d object that is callable
# : rp(E) where E is in [keV]
# ip : same as above, ip(E) where E is in [keV]
CMROOT = os.path.expanduser('~/dev/eblur/dust/data')
class CmDrude(object):
""" OBJECT cmindex.CmDrude
---------------------------------------
__init__(self, rho=3.0)
---------------------------------------
cmtype : 'Drude'
rho : grain density [g cm^-3]
---------------------------------------
FUNCTION
rp(E) : real part of complex index of refraction [E in keV]
ip(E) : imaginary part of complex index of refraction [always 0.0]
"""
def __init__(self, rho=3.0): # Returns a CM using the Drude approximation
self.cmtype = 'Drude'
self.rho = rho
def rp( self, E ):
mm1 = self.rho / ( 2.0*c.mp() ) * c.re()/(2.0*np.pi) * np.power( c.kev2lam()/E , 2 )
return mm1+1
def ip( self, E ):
if np.size(E) > 1:
return np.zeros( np.size(E) )
else:
return 0.0
class CmGraphite(object):
""" OBJECT cmindex.CmGraphite
---------------------------------------
__init__( self, size='big', orient='perp' )
---------------------------------------
cmtype : 'Graphite'
size : 'big' or 'small'
orient : 'perp' or 'para'
rp(E) : scipy.interp1d object
ip(E) : scipy.interp1d object [E in keV]
"""
def __init__( self, size='big', orient='perp' ):
# size : string ('big' or 'small')
# : 'big' gives results for 0.1 um sized graphite grains at 20 K [Draine (2003)]
# : 'small' gives results for 0.01 um sized grains at 20 K
# orient : string ('perp' or 'para')
# : 'perp' gives results for E-field perpendicular to c-axis
# : 'para' gives results for E-field parallel to c-axis
#
self.cmtype = 'Graphite'
self.size = size
self.orient = orient
D03vals = c.restore( os.path.join(CMROOT,'CM_D03.pysav')) # look up file
if size == 'big':
if orient == 'perp':
lamvals = D03vals['Cpe_010_lam']
revals = D03vals['Cpe_010_re']
imvals = D03vals['Cpe_010_im']
if orient == 'para':
lamvals = D03vals['Cpa_010_lam']
revals = D03vals['Cpa_010_re']
imvals = D03vals['Cpa_010_im']
if size == 'small':
if orient == 'perp':
lamvals = D03vals['Cpe_001_lam']
revals = D03vals['Cpe_001_re']
imvals = D03vals['Cpe_001_im']
if orient == 'para':
lamvals = D03vals['Cpa_001_lam']
revals = D03vals['Cpa_001_re']
imvals = D03vals['Cpa_001_im']
lamEvals = c.kev2lam() / c.micron2cm() / lamvals # keV
self.rp = interp1d( lamEvals, revals )
self.ip = interp1d( lamEvals, imvals )
class CmSilicate(object):
""" OBJECT cmindex.CmSilicate
---------------------------------------
__init__( self )
---------------------------------------
cmtype : 'Silicate'
rp(E) : scipy.interp1d object
ip(E) : scipy.interp1d object [E in keV]
"""
def __init__( self ):
self.cmtype = 'Silicate'
D03vals = c.restore( os.path.join(CMROOT,'CM_D03.pysav')) # look up file
lamvals = D03vals['Sil_lam']
revals = D03vals['Sil_re']
imvals = D03vals['Sil_im']
lamEvals = c.kev2lam() / c.micron2cm() / lamvals # keV
self.rp = interp1d( lamEvals, revals )
self.ip = interp1d( lamEvals, imvals )
#------------- A quick way to grab a single CM ------------
def getCM( E, model=CmDrude() ):
""" FUNCTION getCM( E, model=CmDrude() )
---------------------------------------
INPUTS
E : scalar or np.array [keV]
model : any Cm-type object
---------------------------------------
RETURNS
Complex index of refraction : scalar or np.array of dtype='complex'
"""
return model.rp(E) + 1j * model.ip(E)