ComboCode.cc.tools.readers.MlineReader

1 # -*- coding: utf-8 -*- 2 3 """ 4 A class for reading and managing Mline output. 5 6 Author: R. Lombaert and M. van de Sande 7 8 """ 9 10 import os 11 import numpy as np 12 13 from astropy import constants as cst 14 15 from cc.tools.readers.PopReader import PopReader 16 from cc.tools.readers.MolReader import MolReader 17 from cc.tools.io import DataIO 18 19

20 -class MlineReader(MolReader,PopReader):

21 22 ''' 23 A Reader for Mline output files. 24 25 This class inherits from both PopReader and MolReader, and thus provides 26 spectroscopic information of the molecule as given in the mline output. In 27 principle this information should be identical to what is read by 28 RadiatReader from the GASTRoNOoM spectroscopic input files. This information 29 is read from ml1. 30 31 In addition, 32 mline output gives other information calculated by the radiative transfer, 33 and that is available here as well. The level populations specifically are 34 managed through methods in PopReader. This class then also provides 35 additional data such as source function, optical depth, etc. This 36 information is read from ml3. 37 38 MlineReader only depends on GASTRoNOoM output and does not require 39 information from the databases or input files. 40 41 ''' 42

43 - def __init__(self,fn,*args,**kwargs):

44 45 ''' 46 Creating an Mline object ready for reading Mline output. 47 48 Reading mline output through full filename, including filepath. 49 50 The mlinefile number is given as *, which is then replaced as needed. 51 52 The wildcard character does not have to be present, and is inserted 53 automatically. It is handled down the line by the parent classes. 54 55 Additional args and kwargs are passed to the dict creation (parent of 56 Reader) 57 58 Note that properties based on the filename are saved as instance 59 properties (such as id, molecule name, filename). 60 61 @param fn: The mline filename, including filepath. The mline 62 file number is given by *. 63 @type fn: string 64 65 ''' 66 67 #-- Insert wildcard character to encompass all mline output files 68 fn = fn.replace('ml1','ml*').replace('ml2','ml*').replace('ml3','ml*') 69 super(MlineReader, self).__init__(fn=fn,*args,**kwargs) 70 71 #-- Save molecule and id info. 72 self.molecule = os.path.splitext(self.fn)[0].split('_')[-1] 73 self.id = os.path.splitext(self.fn)[0].split('_')[-2].replace('ml*','') 74 75 #-- Read the files. Note that parent classes do this also, but only if 76 # they are not used as parent class through inheritance, and thus work 77 # standalone 78 self.read()

79 80 81

82 - def read(self):

83 84 ''' 85 Read the mline output files ml1 and ml3. 86 87 The information is stored in the instance dictionary. 88 89 ''' 90 91 #-- Read the main output files. 92 self.__readML1() 93 self.__readML3()

94 95 96

97 - def __readML1(self):

98 99 ''' 100 Read molecule spectroscopic properties and basic circumstellar radial 101 profiles from the ml1 file. 102 103 This is the information that is returned by MolReader methods. 104 105 Note that the level indexing is J + 1 for simple molecules like CO. For 106 any molecule, the 0-energy level has index 1! 107 108 Note also that the vel in the ml1 file is the velocity profile divided 109 by the stochastic velocity. (see source_common/vel.f) The velocity 110 profile saved here is the real velocity profile, corrected for this. 111 112 ''' 113 114 #-- Molecule settings are in ml1 115 fn = self.fn.replace('ml*','ml1') 116 117 #-- First block contains keywords, so use readDict 118 d1 = DataIO.readDict(filename=fn,start_row=8,end_row=38,\ 119 convert_floats=1,convert_ints=1,\ 120 key_modifier='lower') 121 self['pars'] = d1 122 123 #-- Number of header lines: number of pars, 9 lines in first header, 1 124 # line with column names. 125 nhdr = len(d1) + 9 + 1 126 127 #-- Second block contains transitions 128 d2 = np.genfromtxt(fn,usecols=[0,1,2,3,4],skip_header=nhdr,\ 129 dtype='i8,i8,i8,f8,f8',max_rows=d1['nline'],\ 130 names=['index','lup','llow','frequency','einsteinA']) 131 132 #-- Add 1 to the lup and llow columns so they match up with the indices 133 # of the levels (for some reason j = 0 level is given as index 0 for 134 # llow in the transition definition, while the level index is actually 135 # 1). While for J-numbers this works, this approach breaks down for eg 136 # H2O. 137 d2['lup'] += 1 138 d2['llow'] += 1 139 self['trans'] = d2 140 141 #-- Third block contains levels (note hdr contains line with col names) 142 nhdr += d1['nline'] + 1 143 d3 = np.genfromtxt(fn,skip_header=nhdr,max_rows=d1['ny'],\ 144 dtype='i8,f8,f8',names=['index','weight','energy']) 145 146 self['level'] = d3 147 148 #-- Fourth block contains circumstellar properties as a function of 149 # impact parameter (note that P1 = R_outer, Pmax = R_star, so reverse) 150 nhdr += d1['ny'] 151 colnames = ['p_rstar','vel','nmol','nh2','amol','Tg','Td'] 152 d4 = np.genfromtxt(fn,names=colnames,\ 153 skip_header=nhdr,max_rows=d1['n_impact'])[::-1] 154 d4['vel'] *= d1['vsto'] 155 self['props'] = d4 156 157 #-- Add a convenient reference to the instance dict so PopReader and 158 # MolReader parent methods can work with the impact parameter grid. 159 # R_STAR is given in cm. 160 self['p'] = d4['p_rstar'] * self['pars']['r_star']

161 162 163

164 - def __readML3(self):

165 166 ''' 167 Read the information from the ML3 file. 168 169 This includes in six blocks, with subblocks for each impact parameter: 170 - Scattering integral (si) 171 - Source function (sf) 172 - Line opacities (lo): if negative (or 1e-60 when 173 use_no_maser_option=1) the line is masing. 174 - Number densities (pop) 175 - Derivative of scat int wrt line opacity times line opacity at line 176 center (DsiDloXlo) 177 (see GASTRoNOoM src code source_mline/lamit.f for details: D1) 178 - Derivative of scat int wrt src function times src function at line 179 center (DsiDsfXsf) 180 (see GASTRoNOoM src code source_mline/lamit.f for details: D2) 181 182 ''' 183 184 #-- 1) First three blocks list values for each transition and each 185 # impact parameter. No white lines. If a white line occurs, it is 186 # because nline%8 == 0. Hence number of text lines is 187 # int(nline)/8+1 per impact parameter. Number of sub-blocks is then 188 # the number of impact parameters. Total length of the block is 189 # nline_l*n_impact. Data starts at previous n_block + 3*N_block, 190 # with N_block in [0->5] for [scat int, sf, line opac, pop, DsiDlo, 191 # DsiDsf]. 3 stands for 2 whitespace lines and one line with dashes 192 # 193 # 2) Data is fixed format, with 14 characters per entry (cannot use 194 # spaces as delimiters, because they might be filled with - sign) 195 # 196 # 3) Data is stored based on level index (for level pops), or trans 197 # index (for everything else), and each index gives an array as a 198 # function of impact parameter, which is in units of R_STAR, and is 199 # given in self['props']['P']. Converted to cgs it is 200 # given in self['props']['P_cm']. 201 # 202 #-- Define the function that reads the fixed-format file per block 203 def readBlock(N_block): 204 205 ''' 206 Read a block of data from ml3. 207 208 @param N_block: The index of the block (0 => 5) 209 @type N_block: int 210 211 @return: The data 212 @rtype: array 213 214 ''' 215 216 #-- Number of lines for this block 217 nblock = nblock_ny if N_block == 3 else nblock_nline 218 219 #-- Set starting index, based on number of block and block length 220 ind0 = (N_block-1)*nblock_nline + 3*N_block 221 222 #-- If N_block > 3, one block has only ny=ny_up+ny_low values, 223 # otherwise add another block with nline values per subblock 224 if N_block > 3: 225 ind0 += nblock_ny 226 else: 227 ind0 += nblock_nline 228 229 #-- Read data and return. fn defined in mother function 230 return np.genfromtxt(fn,skip_header=ind0,max_rows=nblock,\ 231 delimiter=[14]*8)

232 233 #-- Grab filename and set the contents values for the six blocks. 234 fn = self.fn.replace('ml*','ml3') 235 props = ['si','sf','lo','pop','DsiDloXlo','DsiDsfXsf'] 236 for pr in props: 237 self[pr] = dict() 238 239 #-- Gather some relevant parameters. Length of block based on n_impact 240 # and number of values (nline or ny). Assuming 8 columns. 241 n_impact = self['pars']['n_impact'] 242 nline = self['pars']['nline'] 243 nline_l = int(nline)/8+1 244 nblock_nline = n_impact * nline_l 245 ny = self['pars']['ny'] 246 ny_l = int(ny)/8+1 247 nblock_ny = n_impact * ny_l 248 249 #-- Read data, looping over the 6 blocks. 250 dds = [readBlock(N) for N in xrange(6)] 251 252 #-- Set all the transition-specific information. Note indexing (0-based 253 # in python, 1-based in fortran). Note that the arrays are reversed to 254 # match the increasing impact parameter grid. 255 for dd,pr in zip(dds,props): 256 if pr == 'pop': continue 257 for i in xrange(nline): 258 self[pr][i+1] = dd[i/8::nline_l,i%8][::-1] 259 260 #-- Set the level populations. Note indexing (0-based in python, 1-based 261 # in fortran). Note that the arrays are reversed to match the 262 # increasing impact parameter grid. 263 for i in xrange(ny): 264 self['pop'][i+1] = dds[3][i/8::ny_l,i%8][::-1] 265 266 267

268 - def getProp(self,prop):

269 270 ''' 271 Return a radial circumstellar profile associated with the impact 272 parameter grid for the molecule of this mline model. 273 274 The impact parameter grid in cm is available through getP(). 275 276 @param prop: The requested property. One of p_rstar, vel (km/s), nmol 277 (cm^-3), nh2 (cm^-3), amol, Tg (K), Td (K). 278 @type prop: str 279 280 @return: The abundance profile 281 @rtype: array 282 283 ''' 284 285 return self['props'][prop]

286

Source Code for Module ComboCode.cc.tools.readers.MlineReader