1
2
3 """
4 A toolbox for reading dust opacities in every shape and form.
5
6 Author: R. Lombaert
7
8 """
9
10 import os
11 import numpy as np
12 from numpy import array
13 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d
14 from scipy.interpolate import interp1d
15 from astropy import units as u
16 import cc.path
17 from cc.tools.io import DataIO
18
19
21
22 """
23 An interface for reading opacities of multiple dust species.
24
25 Does not inherit from the Reader object, because the files are structured
26 too differently from common input/output data.
27
28 """
29
31
32 """
33 Initiating an instance of the KappaReader.
34
35 """
36
37 self.lspecies = DataIO.getInputData(path=cc.path.usr,\
38 keyword='SPECIES_SHORT',\
39 filename='Dust.dat')
40 self.lfilenames = DataIO.getInputData(path=cc.path.usr,\
41 keyword='PART_FILE',\
42 filename='Dust.dat')
43 self.lspec_dens = DataIO.getInputData(path=cc.path.usr,\
44 keyword='SPEC_DENS',\
45 filename='Dust.dat')
46 self.kappas = dict()
47 self.qext_a = dict()
48 self.waves = dict()
49 self.fns = dict()
50 self.spec_dens = dict()
51
52
53
55
56 """
57 Read kappas (cm2/g) and Q_ext/a (cm-1) for a dust species from the
58 MCMax INPUT files.
59
60 This also reads the absorption and scattering kappas separately.
61
62 @param species: The dust species (from Dust.dat)
63 @type species: string
64
65 """
66
67 if self.waves.has_key(species):
68 return
69 try:
70 ispecies = self.lspecies.index(species)
71 except ValueError:
72 print 'Species not found in Dust.dat.'
73 return
74 fn = os.path.join(cc.path.mopac,self.lfilenames[ispecies])
75 sd = self.lspec_dens[ispecies]
76 if fn[-9:] == '.particle':
77 part_file = DataIO.readFile(filename=fn,delimiter=' ')
78 wav = array([float(q[0])
79 for q in part_file if len(q) == 4])
80 kappa = [array([float(q[1])
81 for q in part_file if len(q) == 4]),
82 array([float(q[2])
83 for q in part_file if len(q) == 4]),
84 array([float(q[3])
85 for q in part_file if len(q) == 4])]
86 else:
87 part_file = DataIO.readCols(filename=fn)
88 wav = part_file[0]
89 kappa = part_file[1:]
90 self.spec_dens[species] = sd
91 self.fns[species] = fn
92 self.waves[species] = wav
93 self.kappas[species] = kappa
94 self.qext_a[species] = array(kappa) * 4/3. * sd
95
96
97
99
100 """
101 Return the wavelength grid for a given species.
102
103 @param species: The dust species (from Dust.dat)
104 @type species: string
105
106 @keyword unit: The unit of the wavelength. Can be given as u.Unit()
107 object or as a string representation of those objects.
108 Can range from length, to frequency, and energy
109
110 (default: micron)
111 @type unit: str/u.Unit()
112
113 @return: wavelength (given unit)
114 @rtype: array
115
116
117 """
118
119
120
121 if isinstance(unit,str) and unit.lower() in ['cm-1','cm^-1']:
122 unit = 1./u.cm
123 elif isinstance(unit,str):
124 unit = getattr(u,unit)
125
126 self.readKappas(species)
127 if not self.waves.has_key(species):
128 return np.empty(0)
129
130 wav = self.waves[species]*u.micron
131
132 if (isinstance(unit,u.Quantity) and unit.unit.is_equivalent(u.K)) \
133 or (isinstance(unit,u.UnitBase) and unit.is_equivalent(u.K)):
134 wav = wav.to(u.erg,equivalencies=u.spectral())
135 return wav.to(unit,equivalencies=u.temperature_energy()).value
136 else:
137 return wav.to(unit,equivalencies=u.spectral()).value
138
139
140
142
143 """
144 Return the kappas for given species.
145
146 The index determines if you want extinction, absorption or scattering.
147
148 @param species: The dust species (from Dust.dat)
149 @type species: string
150
151 @keyword index: The index of the kappas in the .opacity/.particle file.
152 0: extinction, 1: absorption, 2: scattering
153
154 (default: 0)
155 @type index: int
156
157 @return: kappas (cm2/g)
158 @rtype: array
159
160 """
161
162 index = int(index)
163 self.readKappas(species)
164 if self.kappas.has_key(species):
165 return self.kappas[species][index]
166 else:
167 return np.empty(0)
168
169
170
171 - def getExtEff(self,species,index=0):
172
173 """
174 Return the extinction efficiencies per grain size for given species.
175
176 The index determines if you want extinction, scattering or absorption.
177
178 @param species: The dust species (from Dust.dat)
179 @type species: string
180
181 @keyword index: The index of the kappas in the .opacity/.particle file.
182 0: extinction, 1: absorption, 2: scattering
183
184 (default: 0)
185 @type index: int
186
187 @return: q_ext/a [micron,cm-1]
188 @rtype: array
189
190 """
191
192 self.readKappas(species)
193 if self.qext_a.has_key(species):
194 return self.qext_a[species][index]
195 else:
196 return np.empty(0)
197
198
199
200 - def interpolate(self,species,index=0,unit='micron',*args,**kwargs):
201
202 """
203 Create an interpolation object for the mass extinction/absorption/
204 scattering coefficients.
205
206 Additional arguments can be passed to the interpolation object.
207
208 The unit of the wavelength for the interpolation can be chosen.
209
210 @param species: The dust species (from Dust.dat)
211 @type species: string
212
213 @keyword index: The index of the kappas in the .opacity/.particle file.
214 0: extinction, 1: absorption, 2: scattering
215
216 (default: 0)
217 @type index: int
218 @keyword unit: The unit of the wavelength. Can be given as u.Unit()
219 object or as a string representation of those objects.
220 Can range from length, to frequency, and energy
221
222 (default: micron)
223 @type unit: str/u.Unit()
224
225 @return: The interpolator for the mass extinction/absorption/scattering
226 coefficients.
227 @rtype: spline1d
228
229 """
230
231 return spline1d(x=self.getWavelength(species,unit=unit),\
232 y=self.getKappas(species,index),*args,**kwargs)
233