Source Code for Module ComboCode.cc.statistics.ChemStats

  1  # -*- coding: utf-8 -*- 
  2   
  3  """ 
  4  Examination of the Chemistry analysis routine output. 
  5   
  6  Author: M. Van de Sande 
  7   
  8  """ 
  9   
 10  import os 
 11  import scipy 
 12  from scipy import argmin,ones 
 13  from scipy import array 
 14  from scipy import sum 
 15  import operator 
 16  import types 
 17  import numpy as np 
 18   
 19  import cc.path 
 20  from cc.tools.io import DataIO 
 21  #from cc.modeling.tools import CodeIO 
 22  from cc.statistics.Statistics import Statistics 
 23  import matplotlib.pyplot as plt 
 24  from cc.modeling.objects import Transition 
 25  from cc.statistics import BasicStats     
 26  from time import gmtime 
 27   
 28 -class ChemStats(object): 
 29       
 30      """ 
 31      Environment with several tools to perform statistics on the peak and  
 32      integrated intensities of resolved molecular lines. 
 33       
 34      """ 
 35           
 36 -    def __init__(self,star_name,path_code='OutputClumpy',\ 
 37                   star_grid=[],molecules = []):         
 38           
 39          """  
 40           
 41          @param star_name: Star name from Star.dat 
 42          @type star_name: string 
 43           
 44          @keyword code: the code used for producing your output  
 45           
 46                         (default: 'GASTRoNOoM') 
 47          @type code: string 
 48          @keyword path_code: Output folder in the code's home folder 
 49                          
 50                              (default: 'codeSep2010') 
 51          @type path_code: string 
 52          @keyword molecules: Molecules for which the analysis output needs to be 
 53                              retrieved 
 54                              (default: []) 
 55          @type molecules: list 
 56                   
 57          """ 
 58           
 59          #super(ChemStats,self).__init__(star_name=star_name,\ 
 60                                         #code=code,path_code=path_code) 
 61   
 62          self.path_code = path_code 
 63          self.molecules = molecules 
 64           
 65          #-- Keep track of each Star() model that is valid. 
 66          self.star_grid = dict() 
 67          self.star_grid = [star for star in star_grid 
 68                            if star['LAST_CHEMISTRY_MODEL'] 
 69                            and star['PERFORM_ROUTINE']==1] 
 70           
 71          self.molecules = molecules 
 72          #-- Convenience path 
 73          cc.path.cout = os.path.join(cc.path.chemistry,self.path_code) 
 74           
 75           
 76 -    def analyseMolecule(self, star, molec, print_main=1, print_sec=0): 
 77           
 78          analyse = self.readAnalyse(star) 
 79           
 80          if print_main == 1: 
 81              self.printMainPathways(analyse,molec) 
 82           
 83          if print_sec == 1: 
 84              sec = self.getSecondaryPathways(analyse,molec,print_sec=1) 
 85           
 86           
 87           
 88 -    def readAnalyse(self, star): 
 89           
 90          ''' 
 91          Reads the output of the analyse routine performed. 
 92           
 93          @keyword star: Star object containing the chemistry model 
 94          @type star: Star() 
 95           
 96          @return: Dictionary containing the analyse output for every 
 97                           molecule, eg analyse['SiO']. The analyse radius is  
 98                           also included (analyse['radius']). 
 99          @rtype: dict() 
100           
101          ''' 
102           
103          #- Read in analyse.out 
104          filename = os.path.join(cc.path.cout,'models',\ 
105                                  star['LAST_CHEMISTRY_MODEL'],'analyse.out') 
106          data = DataIO.readFile(filename) 
107          data = [x.split() for x in data] 
108           
109          #- Initialise output dictionary 
110          analyse = dict() 
111           
112          #- Determine and save routine radius 
113          radius = float(data[0][-2]) 
114          index = np.where([len(i)==2 for i in data])[0] 
115          analyse['radius'] = radius 
116           
117          #- Save analyse routine output in dictionary per molecule 
118          for ii in range(len(index[:-1])): 
119              name = data[index[ii]][1][1:] 
120              analyse[name] = dict() 
121              analyse[name]['MAIN'] = data[index[ii]+1:index[ii+1]-1] 
122              analyse[name]['DRATE'] = data[index[ii+1]-1][1] 
123              analyse[name]['PRATE'] = data[index[ii+1]-1][3] 
124           
125          return analyse 
126   
127       
128 -    def getSecondaryPathways(self,analyse,molec,print_sec=1,\ 
129                               destruction=1,formation=1): 
130           
131          others = analyse.keys() 
132          others.remove('radius') 
133          others.remove(molec) 
134           
135          sec = dict() 
136           
137          for x in others: 
138              sec[x] = dict() 
139              if destruction == 1 and formation == 1: 
140                  dummy = [st for st in analyse[x]['MAIN'] if molec in st] 
141              elif destruction == 1 and formation == 0: 
142                  dummy = [st for st in analyse[x]['MAIN'] if molec in st[:3]] 
143              elif destruction == 0 and formation == 1: 
144                  dummy = [st for st in analyse[x]['MAIN'] if molec in st[3:]] 
145              if dummy != []: 
146                  sec[x]['MAIN'] = dummy 
147                  sec[x]['PRATE'] = analyse[x]['PRATE'] 
148                  sec[x]['DRATE'] = analyse[x]['DRATE'] 
149              else: 
150                  sec.pop(x) 
151           
152          if print_sec == 1: 
153              print '***   Minor pathways at radius '+str(analyse['radius'])+', '+molec+'  ***'          
154              for x in sec.keys(): 
155                  print '\t ***      '+x+'       ***' 
156                  for st in sec[x]['MAIN']: 
157                      print '\t'.join(st) 
158                  print '\t '+'DRATE = '+sec[x]['DRATE']+\ 
159                      '\t'+'PRATE = '+sec[x]['PRATE']+'\n' 
160   
161          return sec 
162           
163   
164 -    def printMainPathways(self,analyse,molec): 
165           
166          print '***   Main pathways at radius '+str(analyse['radius'])+', '+molec+'  ***'          
167          for m in analyse[molec]['MAIN']: 
168              print '\t'.join(m) 
169          print '\t '+'DRATE = '+analyse[molec]['DRATE']+\ 
170                  '\t'+'PRATE = '+analyse[molec]['PRATE']+'\n' 
171   
172   
173 -    def printSecondaryPathways(self,analyse,molec,\ 
174                               destruction=1,formation=1): 
175           
176          others = analyse.keys() 
177          others.remove('radius') 
178          others.remove(molec) 
179           
180          sec = dict() 
181           
182          for x in others: 
183              sec[x] = dict() 
184              if destruction == 1 and formation == 1: 
185                  dummy = [st for st in analyse[x]['MAIN'] if molec in st] 
186              elif destruction == 1 and formation == 0: 
187                  dummy = [st for st in analyse[x]['MAIN'] if molec in st[:3]] 
188              elif destruction == 0 and formation == 1: 
189                  dummy = [st for st in analyse[x]['MAIN'] if molec in st[3:]] 
190              if dummy != []: 
191                  sec[x]['MAIN'] = dummy 
192                  sec[x]['PRATE'] = analyse[x]['PRATE'] 
193                  sec[x]['DRATE'] = analyse[x]['DRATE'] 
194              else: 
195                  sec.pop(x) 
196           
197          print '***   Minor pathways at radius '+str(analyse['radius'])+', '+molec+'  ***'          
198          for x in sec.keys(): 
199              print '\t ***      '+x+'       ***' 
200              for st in sec[x]['MAIN']: 
201                  print '\t'.join(st) 
202              print '\t '+'DRATE = '+sec[x]['DRATE']+\ 
203                  '\t'+'PRATE = '+sec[x]['PRATE']+'\n' 
204   
205           
206      #def printMainPathways(self,star): 
207           
208          ##for star in self.star_grid: 
209          #analyse = self.readAnalyse(star) 
210       
211          #print '\n' 
212          #print '***************************************************' 
213          #print '***   Main formation and destruction pathways   ***' 
214          #print '*** \t \t at radius '+str(analyse['radius'])+' \t  \t***' 
215          #print '***************************************************' 
216           
217          #for molec in self.molecules: 
218              #print '\t ***      '+molec+'       ***' 
219              #for m in analyse[molec]['MAIN']: 
220                  #print '\t'.join(m) 
221              #print '\t '+'DRATE = '+analyse[molec]['DRATE']+\ 
222                      #'\t'+'PRATE = '+analyse[molec]['PRATE']+'\n' 
223