1
2
3 """
4 A tool for comparing continuum divided data and models in a given wavelength
5 range. Useful for fitting dust emission and absorption features.
6
7 Author: R. Lombaert
8
9 """
10
11 import os,types
12 from scipy.integrate import trapz
13
14 import cc.path
15 from cc.tools.numerical import Interpol
16 from cc.data import Data
17 from cc.modeling.codes import MCMax
18 from cc.plotting import Plotting2
19 from cc.tools.io import DataIO
20
22
23 '''
24 A class for dividing a spectrum by the continuum and comparing models and
25 data.
26
27 '''
28
29 - def __init__(self,star_grid=[],spec=[],franges=[2.6,2.85,3.3,3.7],plot=0,\
30 func='power',cfg=''):
31
32 '''
33 Initializing a ContinuumDivision instance.
34
35 @keyword star_grid: The parameter sets for which the fitting is done
36
37 (default: [])
38 @type star_grid: list[Star()]
39 @keyword spec: information concerning dust data available for star
40
41 (default: [])
42 @type spec: list[Sed()]
43 @keyword franges: The fitting ranges for this instance. 4 values in mic,
44 delimiting 2 parts of the spectrum blueward and
45 redward of the dust feature. Default is for 3.1
46 micron ice feature.
47 Can also give a list of lists, in which case the
48 franges are defined for every Star() and Sed()
49 included, such that
50 len(franges) == len(spec)+len(star_grid). The model
51 franges are listed first, then the sed franges.
52 Can also be given through the cfg file, and takes
53 priority in that case.
54
55 (default: [2.6,2.85,3.3,3.7])
56 @type franges: list[float] or list[list[float]]
57 @keyword plot: Show the continuum division and fitting plots.
58
59 (default: 0)
60 @type plot: bool
61 @keyword func: The function used for fitting the continuum. Can be a
62 list of functions as well, of len == len(star_grid) +
63 len(spec)
64 Can also be given through the cfg file, and takes
65 priority in that case.
66
67 (default: power)
68 @type func: string or list[string]
69 @keyword cfg: a configuration file for Plotting2.py.
70 @type cfg: string
71
72 '''
73
74 self.star_grid = star_grid
75 self.spec = spec
76 self.cont_division = dict()
77 self.eq_width = dict()
78 self.plot = plot
79 self.cfg = cfg
80 if cfg:
81 cfg_dict = DataIO.readDict(cfg,convert_lists=1,convert_floats=1)
82 else:
83 cfg_dict = dict()
84 if cfg_dict.has_key('franges'):
85 franges = cfg_dict['franges']
86 if type(franges[0]) is types.TupleType:
87 franges = [list(f) for f in franges]
88 if cfg_dict.has_key('func'):
89 func = cfg_dict['func']
90 if type(franges[0]) is types.ListType:
91 if not len(franges) == (len(spec) + len(star_grid)):
92 print 'Not enough sets of franges defined for all Star() and'+\
93 ' Sed() objects. Taking first set for all.'
94 self.franges = [sorted(franges[0])
95 for i in range(len(spec)+len(star_grid))]
96 else:
97 self.franges = [sorted(fr) for fr in franges]
98 else:
99 self.franges = [sorted(franges)
100 for i in range(len(spec)+len(star_grid))]
101
102 self.frmin = min([min(fr) for fr in self.franges])
103 self.frmax = max([max(fr) for fr in self.franges])
104 if type(func) is types.ListType:
105 if not len(func) == (len(spec) + len(star_grid)):
106 print 'Not enough functions defined for all Star() and'+\
107 ' Sed() objects. Taking first func for all.'
108 self.func = [func[0] for i in range(len(spec)+len(star_grid))]
109 else:
110 self.func = func
111 else:
112 self.func = [func
113 for i in range(len(spec)+len(star_grid))]
114
115
116
118
119 '''
120 Show the model calculations in a plot.
121
122 '''
123
124 x = []
125 y = []
126 keys = []
127 for i,sed in enumerate(self.spec):
128 x.append(self.cont_division['sws%i'%i]['w_feat'])
129 y.append(self.cont_division['sws%i'%i]['f_division'])
130 keys.append(sed.star_name_plots)
131 for i,star in enumerate(self.star_grid):
132 if star['LAST_MCMAX_MODEL']:
133 x.append(self.cont_division[star['LAST_MCMAX_MODEL']]\
134 ['w_feat'])
135 y.append(self.cont_division[star['LAST_MCMAX_MODEL']]\
136 ['f_division'])
137
138
139
140 keys.append(star['LAST_MCMAX_MODEL'].replace('_','\_'))
141
142
143
144
145
146 fn = Plotting2.plotCols(x=x,y=y,xmin=self.frmin*0.9,\
147 xmax=self.frmax*1.1,keytags=keys,\
148 key_location=(0.6,0.02),cfg=self.cfg)
149 if not fn is None:
150 print 'Your Continuum Division plot can be found at '
151 print fn
152
153
154
156
157 '''
158 Calculate the equivalent width after continuum division.
159
160 @param dtype: data type (only 'model','sws' for now)
161 @type dtype: string
162 @param frindex: The index in the franges list for this entry.
163 @type frindex: int
164
165 '''
166
167 fr1 = self.franges[frindex][0]
168 fr4 = self.franges[frindex][3]
169 w = self.cont_division[dtype]['w_feat']
170 w_feat = w[(w>fr1)*(w<fr4)]
171 f_feat = self.cont_division[dtype]['f_division'][(w>fr1)*(w<fr4)]
172 self.eq_width[dtype] = trapz(x=w_feat,y=(1-f_feat))
173
174
176
177 '''
178 Prepare models for dust feature continuum division.
179
180 The equivalent width of the feature is also calculated here.
181
182 '''
183
184 for i,s in enumerate(self.star_grid):
185 model_id = s['LAST_MCMAX_MODEL']
186 if not model_id: continue
187 dpath = os.path.join(cc.path.mout,'models',model_id)
188 fn_spec = 'spectrum{:04.1f}.dat'.format(s['RT_INCLINATION'])
189 w,f = MCMax.readModelSpectrum(dpath,s['RT_SPEC'],fn_spec)
190 self.divideContinuum(w,f,dtype=model_id,frindex=i)
191 self.calcEqWidth(dtype=model_id,frindex=i)
192
193
195
196 '''
197 Prepare data for dust feature continuum division.
198
199 For now only SWS is available.
200
201 The equivalent width of the feature is also calculated here.
202
203 '''
204
205 for i,sp in enumerate(self.spec):
206 sws_type = [k for k in sp.data.keys() if 'SWS' in k][0]
207 w = sp.data[sws_type][0]
208 f = sp.data[sws_type][1]
209 self.divideContinuum(w,f,dtype='sws%i'%i,\
210 frindex=len(self.star_grid)+i)
211 self.calcEqWidth(dtype='sws%i'%i,frindex=len(self.star_grid)+i)
212
213
215
216 '''
217 Divide flux by the continuum flux in a dust feature.
218
219 @param w: The wavelength grid
220 @type w: list/array
221 @param f: The flux grid
222 @type f: list/array
223 @param dtype: data type (only 'model','sws' for now)
224 @type dtype: string
225 @param frindex: The index in the franges list for this entry.
226 @type frindex: int
227
228 @keyword plot: Show a plot of the continuum division
229 @type plot: bool
230
231 '''
232
233 fr1,fr2 = self.franges[frindex][0],self.franges[frindex][1]
234 fr3,fr4 = self.franges[frindex][2],self.franges[frindex][3]
235 w_in = list(w[(w>fr1) * (w<fr2)]) + list(w[(w>fr3) * (w<fr4)])
236 f_in = list(f[(w>fr1) * (w<fr2)]) + list(f[(w>fr3) * (w<fr4)])
237 w_cont = w[(w>self.frmin*0.9)*(w<self.frmax*1.1)]
238 f_ori = f[(w>self.frmin*0.9)*(w<self.frmax*1.1)]
239 f_cont = Interpol.fitFunction(x_in=w_in,y_in=f_in,x_out=w_cont,\
240 func=self.func[frindex])
241
242 self.cont_division[dtype] = dict()
243 self.cont_division[dtype]['w_feat'] = w_cont
244 self.cont_division[dtype]['f_feat'] = f_ori
245 self.cont_division[dtype]['w_fitsel_feat'] = w_in
246 self.cont_division[dtype]['f_fitsel_feat'] = f_in
247 self.cont_division[dtype]['f_interp'] = f_cont
248 self.cont_division[dtype]['f_division'] = f_ori/f_cont
249
250 if self.plot:
251 x = [w_cont,w_cont]
252 y = [f_ori,f_cont]
253 Plotting2.plotCols(x=x,y=y,xmin=self.frmin*0.9,xmax=self.frmax*1.1,\
254 ylogscale=0,xlogscale=0)
255