Source Code for Module ComboCode.cc.modeling.tools.DustOpacity

  1  # -*- coding: utf-8 -*- 
  2   
  3  """ 
  4  Tools for making and managing opacity files. 
  5   
  6  Author: R. Lombaert 
  7   
  8  """ 
  9   
 10  from scipy import array,argmin 
 11  from glob import glob 
 12  from math import pi 
 13  import os 
 14  import types 
 15   
 16  import cc.path 
 17  from cc.tools.io import DataIO 
 18  from cc.tools.numerical import Interpol 
 19   
 20   
 21 -def massFractionGSD(acut,amin=0.01,amax=100.,slope=-3.5): 
 22       
 23      ''' 
 24      Calculate the mass fraction of a subset of grain sizes from a grain size 
 25      distribution. 
 26       
 27      A minimum of two subsets is created, so two masses are returned. Can be  
 28      anything higher than that, depending on how many cuts are requested. 
 29       
 30      @param acut: The grain size(s) at which the cuts are made. Can be a float  
 31                   or a list.  
 32      @type acut: list(float) 
 33       
 34      @keyword amin: The minimum grain size in the distribution 
 35                       
 36                     (default: 0.01) 
 37      @type amin: float 
 38      @keyword amax: The maximum grain size in the distribution 
 39                       
 40                     (default: 100) 
 41      @type amax: float 
 42      @keyword slope: The slope of the a-dependence of the distribution. Default 
 43                      is the MRN value of -3.5 
 44                       
 45                      (default: -3.5) 
 46      @type slope: float 
 47       
 48      @return: The mass fraction of the first slice, the second slice, etc,  
 49               depending on how many cuts are requested. Minimum of two. 
 50      @rtype: [float,float,...]          
 51       
 52      ''' 
 53       
 54      if not type(acut) is types.ListType: 
 55          acut = [float(acut)] 
 56      acut = sorted(acut) 
 57      #-- Slope is negative typically. If slope == -4, this doesn't work 
 58      if slope == -4.: return np.empty(0) 
 59      #-- rho(a) = rho_spec * V(a) * n(a) = Cst * Int(a**3 * a**-3.5 da) ==> Cst * [a**0.5]^amax_amin 
 60      p = (3.+slope)+1. 
 61      subsets = [] 
 62      #-- first subset: 
 63      subsets.append(acut[0]**p-amin**p)     
 64      #-- additional subsets if len(acut) > 1 
 65      for i in xrange(len(acut[1:])): 
 66          subsets.append(acut[i+1]**p-acut[i]**p) 
 67      #-- final subset 
 68      subsets.append(amax**p-acut[-1]**p) 
 69      #-- Determine the fractions with respect to the total number density. 
 70      #   Constant does not matter because of it being ratios. 
 71      fractions = array(subsets)/sum(subsets) 
 72      return fractions 
 73       
 74   
 75   
 76 -def mergeOpacity(species,lowres='nom_res',highres='high_res'): 
 77       
 78      ''' 
 79      Merge high-res opacities into a grid of low-res opacities. 
 80       
 81      The wavelength range of the inserted high res opacities is taken from the  
 82      given high res grid.  
 83           
 84      @param species: The dust species for which this is done. This is also the  
 85                      name of the folder in ~/MCMax/DustOpacities/ that contains 
 86                      the data files. 
 87      @type species: string 
 88      @keyword lowres: The subfolder in ~/MCMax/DustOpacities/species containing 
 89                       the low resolution datafiles. 
 90                      
 91                       (default: low_res) 
 92      @type lowres: string 
 93      @keyword highres: The subfolder in ~/MCMax/DustOpacities/species containing 
 94                        the high resolution datafiles. 
 95                         
 96                        (default: high_res) 
 97      @type highres: string 
 98       
 99      ''' 
100       
101      path = os.path.join(cc.path.mopac,species) 
102      lowres_files = [f  
103                      for f in glob(os.path.join(path,lowres,'*'))  
104                      if f[-5:] == '.opac'] 
105      highres_files = [f  
106                       for f in glob(os.path.join(path,highres,'*'))  
107                       if f[-5:] == '.opac'] 
108      files = set([os.path.split(f)[1] for f in lowres_files] + \ 
109                  [os.path.split(f)[1] for f in highres_files]) 
110       
111      for f in files: 
112          hdfile = os.path.join(path,highres,f) 
113          ldfile = os.path.join(path,lowres,f) 
114          if os.path.isfile(ldfile) and os.path.isfile(hdfile): 
115              hd = DataIO.readCols(hdfile) 
116              ld = DataIO.readCols(ldfile) 
117              hdw = hd[0] 
118              ldw = ld[0] 
119              wmin = hdw[0] 
120              wmax = hdw[-1] 
121              ld_low = [list(col[ldw<wmin]) for col in ld] 
122              ld_high = [list(col[ldw>wmax]) for col in ld] 
123              hd = [list(col) for col in hd] 
124              merged = [ld_low[i] + hd[i] + ld_high[i]  
125                        for i in range(len(hd))] 
126              DataIO.writeCols(filename=os.path.join(path,f),cols=merged) 
127   
128   
129 -class CustomOpacity(): 
130      """ 
131      An interface for creating custom opacity files by taking the original and  
132      tinkering with it. 
133       
134      """ 
135       
136 -    def __init__(self,species): 
137           
138          """ 
139          Initializing an instance of the custom opacities interface. 
140       
141          @param species: the species short name 
142          @type species: string 
143           
144          """ 
145           
146          self.setSpecies(species) 
147          self.opacity_file = False 
148   
149       
150 -    def setSpecies(self,species): 
151           
152          """ 
153          Change the species of the current CustomOpacity instance. 
154           
155          @param species: the species short name 
156          @type species: string 
157           
158          """ 
159           
160          self.species = species 
161          self.index = DataIO.getInputData(keyword='SPECIES_SHORT',\ 
162                                           filename='Dust.dat')\ 
163                                          .index(self.species) 
164          self.filename =  DataIO.getInputData(keyword='PART_FILE',\ 
165                                               filename='Dust.dat',\ 
166                                               rindex=self.index) 
167          fn = os.path.join(cc.path.mopac,self.filename) 
168          self.input_data = DataIO.readFile(filename=fn,delimiter=' ')  
169           
170           
171           
172 -    def makeOpa(self,mode='ZERO',**args): 
173           
174          """ 
175          Making custom .particle files. 
176           
177          Every method called here will put the results in self.output_data. 
178           
179          @keyword mode: type of extrapolation (ZERO,FUNCTION,HONY,ZHANG) 
180                           
181                         (default: 'ZERO') 
182          @type mode: string 
183          @keyword args: Optional keywords required for the other methods of the  
184                         class 
185          @type args: dict 
186           
187          """ 
188           
189          self.output_data = [] 
190          mode = mode.upper() 
191          if hasattr(self,'do' + mode): 
192              getattr(self,'do' + mode)(**args) 
193              self.output_data = [' '.join(str(line))  
194                                  for line in self.output_data] 
195              output_filename = '_'.join(['customOpacity',mode] + \ 
196                                         sorted(args.values()) + \ 
197                                         [self.filename]) 
198              if self.opacity_file: 
199                  output_filename.replace('.particle','.opacity') 
200              DataIO.writeFile(filename=os.path.join(cc.path.mopac,\ 
201                                                     output_filename),\ 
202                               input_lines=self.output_data) 
203              new_short = self.species + mode 
204              #- filename is already present: the new file should have the same 
205              #- parameters and the short name can be kept,  
206              #- nothing is changed in the Dust.dat file 
207              try:     
208                  DataIO.getInputData(keyword='PART_FILE',filename='Dust.dat')\ 
209                                     .index(output_filename) 
210              #- filename is not present: do the normal procedure, ie check if  
211              #- short name is already present 
212              except ValueError:         
213                  i=0 
214                  while ' '.join(DataIO.getInputData(keyword='SPECIES_SHORT',\ 
215                                                     filename='Dust.dat'))\ 
216                                                    .find(new_short) != -1: 
217                      i+=1     
218                      new_short = new_short + str(i) 
219                  adding_line = [new_short] + \ 
220                                [str(DataIO.getInputData(keyword=key,\ 
221                                                         filename='Dust.dat',\ 
222                                                         rindex=self.index)) 
223                                 for key in ['SPEC_DENS','T_DES','T_DESA','T_DESB']] 
224                  adding_line.insert(2,output_filename) 
225                  adding_line = '\t\t'.join(adding_line) 
226                  DataIO.writeFile(os.path.join(cc.path.usr,'Dust.dat'),\ 
227                                   [adding_line+'\n'],mode='a') 
228          else: 
229              print 'Mode "' + mode + '" not available. Aborting.' 
230           
231           
232           
233 -    def doZERO(self,wl=10.0): 
234           
235          """ 
236          Replace all opacities by zero below a cut-off wavelength. 
237       
238          @keyword wl: the wavelength of the discontinuity to zero 
239           
240                       (default: 10.0)  
241          @type wl: float 
242           
243          """ 
244           
245          for line in self.input_data: 
246              if len(line) == 4 and float(line[0]) < wl: 
247                  self.output_data.append([line[0],str(2e-60),str(1e-60),str(1e-60)]) 
248              else: 
249                  self.output_data.append(line) 
250   
251   
252   
253 -    def doFUNCTION(self,wl_min=10.0,wl_max=20.0,function='linear'): 
254           
255          """ 
256          Replace all opacities by an extrapolation below a cut-off wavelength. 
257       
258          @keyword wl_min: lower boundary interpolation range and the cutoff  
259                           wavelength for the extrapolation 
260                            
261                           (default: 10.0) 
262          @type wl_min: float 
263          @keyword wl_max: upper boundary interpolation range 
264           
265                           (default: 20.0) 
266          @type wl_max: float 
267          @keyword function: type of function used for the inter/extrapolation 
268                             See Interpol.pEval for function types     
269                           
270                             (default: 'linear') 
271          @type function: string  
272           
273          """ 
274           
275          #raise TypeError("WARNING! The CustumOpa().doFUNCTION method is obsolete! New version NYI.") 
276           
277          #- Select relevant inputlines (not saving the scattering matrices) 
278          self.opacity_file = True 
279          inputsel = [line for line in self.input_data if len(line) == 4] 
280          inputsel = array([[float(f) for f in line] for line in inputsel]) 
281          wl = inputsel[:,0] 
282          function = function.lower() 
283           
284          #- Select the extrapolation and interpolation regions. 
285          wl_low = wl[wl < wl_min] 
286          wlinter = wl[(wl >= wl_min) * (wl <= wl_max)] 
287          abs_inter = inputsel[:,2][(wl >= wl_min) * (wl <= wl_max)] 
288          sca_inter = inputsel[:,3][(wl >= wl_min) * (wl <= wl_max)] 
289           
290          #- Fit the interpolation region and extrapolate this to lower wavelength 
291          abs_p0 = [abs_inter[0],wl[0],1,1] 
292          sca_p0 = [sca_inter[0],wl[0],1,1] 
293          abs_low = Interpol.fitFunction(x_in=wlinter,y_in=abs_inter,\ 
294                                         x_out=wl_low,func=function,\ 
295                                         initial=abs_p0,show=1)         
296          abs_low = array([a > 0 and a or 0 for a in abs_low]) 
297          sca_low = Interpol.fitFunction(x_in=wlinter,y_in=sca_inter,\ 
298                                         x_out=wl_low,func=function,\ 
299                                         initial=sca_p0,show=1) 
300          sca_low = array([s > 0 and s or 0 for s in sca_low]) 
301          ext_low = abs_low + sca_low 
302           
303          #- Set the output data 
304          self.output_data = [array([w,el,al,sl]) 
305                              for w,el,al,sl in zip(wl_low,ext_low,\ 
306                                                    abs_low,sca_low)] 
307          self.output_data.extend(inputsel[wl >= wl_min]) 
308                               
309                               
310                               
311 -    def doHONY(self,wl1 = 1.0,wl2=2.0,wl3=10.0,a_mod=0.05,q_cst=1.0): 
312           
313          """ 
314          Replace all opacities below a cut-off wavelength with values as assumed  
315          by Hony et al (2003, 2004) 
316       
317          @keyword wl1: first discontinuity in micron 
318           
319                        (default: 1.0) 
320          @type wl1: float 
321          @keyword wl2: second discontinuity in micron 
322                         
323                        (default: 2.0) 
324          @type wl2: float 
325          @keyword wl3: third discontinuity in micron 
326           
327                        (default: 10) 
328          @type wl3: float 
329          @keyword a_mod: model grain size, 0.01 according to Hony 2004 in micron 
330           
331                          (default: 0.05) 
332          @type a_mod: float 
333          @keyword q_cst: constant extinction efficiency at wavelengths < wl1 
334           
335                          (default: 1.0) 
336          @type q_cst: float 
337           
338          """ 
339   
340          spec_dens = DataIO.getInputData(keyword='SPEC_DENS',rindex=self.index,\ 
341                                          filename='Dust.dat') 
342          opa_cst = q_cst/4.0*3.0/spec_dens/(a_mod*10**(-4))       
343          for line in self.input_data: 
344              if len(line) == 4 and float(line[0]) < wl1: 
345                  self.output_data.append([line[0],str(opa_cst+1e-60),\ 
346                                           str(opa_cst),str(1e-60)]) 
347              elif len(line)==4 and float(line[0])>=wl1 and float(line[0]) < wl2: 
348                  opa = Interpol.linInterpol([wl1,wl2],[opa_cst,1e-60],\ 
349                                              float(line[0])) 
350                  self.output_data.append([line[0],str(opa+1e-60),str(opa),\ 
351                                           str(1e-60)]) 
352              elif len(line)==4 and float(line[0])>=wl2 and float(line[0]) < wl3: 
353                  self.output_data.append([line[0],str(2e-60),str(1e-60),\ 
354                                           str(1e-60)]) 
355              else: 
356                  self.output_data.append(line) 
357            
358            
359            
360 -    def doZHANG(self,wl=10.0,a_mod=0.05,q_cst=1.6,metallic=True): 
361           
362          """ 
363          Replace all opacities below a cut-off wavelength with values as assumed 
364          by Zhang et al (2009) 
365       
366          @keyword wl: wavelength of the discontinuity 
367           
368                       (default: 10.0) 
369          @type wl: float 
370          @keyword a_mod: model grain size, following Zhang et al 2009 
371           
372                          (default: 0.05) 
373          @type a_mod: float 
374          @keyword q_cst: the constant extinction efficieny at short wavelengths 
375           
376                          (default: 1.6) 
377          @type q_cst: float 
378          @keyword metallic: if metallic dust, cut off wavelength is calculated  
379                             differently than if dielectric dust (ie metallic=0) 
380                              
381                             (default: True) 
382          @type metallic: bool 
383           
384          """ 
385           
386          spec_dens = DataIO.getInputData(keyword='SPEC_DENS',rindex=self.index,\ 
387                                          filename='Dust.dat') 
388          wl1 = metallic and pi*a_mod or 2*pi*a_mod 
389          opa_cst = q_cst/4.0*3.0/spec_dens/(a_mod*10**(-4)) 
390          i = 0 
391          while float(self.input_data[i][0]) <= wl: 
392              i += 1 
393          opa_10 = float(self.input_data[i][1]) 
394          for line in self.input_data: 
395              if len(line) == 4 and float(line[0]) <= wl1: 
396                  self.output_data.append([line[0],str(opa_cst+1e-60),\ 
397                                           str(opa_cst),str(1e-60)]) 
398              elif len(line)==4 and float(line[0]) > wl1 and float(line[0]) < wl: 
399                  eps = min(1,float(line[0])/wl) 
400                  opa = (1-eps) * opa_cst*wl1/float(line[0]) + eps*opa_10 
401                   
402                  self.output_data.append([line[0],str(opa+1e-60),str(opa),\ 
403                                           str(1e-60)]) 
404              else: 
405                  self.output_data.append(line) 
406