ComboCode.cc.modeling.tools.ColumnDensity

22 23 """ 24 Environment to calculate column densities and similar properties for a 25 Star() object. 26 27 """ 28

29 - def __init__(self,star):

30 31 """ 32 Initializing an instance of the ColumnDensity class. 33 34 To be calculated: 35 - r_des: The minimum radius, ie where the dust species starts to exist. 36 This is usually where the density of the species rises to 37 1% of the maximum density of the shell (in cm). 38 - t_des: The temperature at the destruction radius r_des of the species 39 - r_max: The maximum radius at which species exists, taken to be 40 outer radius, or where the density of species drops below 41 10**(-10) times the maximum density in the shell (in cm). 42 - t_min: The temperature at the maximum radius r_max of the species 43 - coldens: The column density in g/cm2 of the dust species between 44 r_des and r_max 45 46 @param star: The model parameter set 47 @type star: Star() 48 49 """ 50 51 self.star = star 52 self.Rsun = star.Rsun #in cm 53 self.Msun = star.Msun #in g 54 self.avogadro = cst.N_A.cgs.value 55 self.au = star.au 56 self.r_des = dict() 57 self.r_max = dict() 58 self.t_min = dict() 59 self.t_des = dict() 60 self.r_min_cd = dict() 61 self.r_max_cd = dict() 62 self.compd = dict() 63 self.rad = [] 64 self.coldens = dict() 65 self.fullcoldens = dict() 66 self.dustfractions = dict() 67 self.readDustInfo() 68 if int(self.star['MRN_DUST']): 69 raise IOError('No column densities can be calculated for now if MRN distribution is used in MCMax.')

70 71

72 - def readDustInfo(self):

73 74 """ 75 Read all column densities, min/max temperatures and min/max radii for 76 the species involved in the MCMax model. 77 78 Note that the self.coldens dictionary does not give real column 79 densities! This dict merely gives column densities in a prescribed 80 shell with given min and max radius, in order to compare with the H2 81 col density. 82 83 """ 84 85 dens = self.star.getDustDensity() 86 temp = self.star.getDustTemperature() 87 compf = os.path.join(cc.path.mcmax,self.star.path_mcmax,'models',\ 88 self.star['LAST_MCMAX_MODEL'],'composition.dat') 89 comp = DataIO.readCols(compf) 90 self.rad = comp.pop(0)*self.au 91 self.r_outer = self.rad[-1] 92 93 for species in self.star.getDustList(): 94 #- Save the actual density profile for this dust species, as well 95 #- as calculating the full column density of a dust species. 96 self.dustfractions[species] = comp.pop(0) 97 self.compd[species] = self.dustfractions[species]*dens 98 self.fullcoldens[species] = trapz(x=self.rad,y=self.compd[species]) 99 #- Determine the column density from 90% of the dust species formed 100 #- onward, based on the mass fractions! 101 #- Not before, because the comparison with H2 must be made, 102 #- and this will skew the result if not solely looking at where the 103 #- dust has (almost) all been formed. 104 #- We also save min amd max radii, for use with the H2 calculation 105 a_species = self.star['A_%s'%species] 106 maxdens = max(self.compd[species]) 107 mindens = maxdens*10**(-10) 108 radsel = self.rad[(self.dustfractions[species]>0.9*a_species)*\ 109 (self.compd[species]>mindens)] 110 denssel = self.compd[species]\ 111 [(self.dustfractions[species]>0.9*a_species)*\ 112 (self.compd[species]>mindens)] 113 self.coldens[species] = trapz(x=radsel,y=denssel) 114 if radsel.size: 115 self.r_min_cd[species] = radsel[0] 116 self.r_max_cd[species] = radsel[-1] 117 else: 118 print 'Threshold dust mass fraction not reached for %s.'%species 119 self.r_min_cd[species] = 0 120 self.r_max_cd[species] = 0 121 #- Determine the actual destruction radius and temperature. 122 #- Taken where the density reaches 1% of the maximum density 123 #- (not mass fraction). 124 self.r_des[species] = self.rad[self.compd[species]>(maxdens*0.01)][0] 125 self.t_des[species] = temp[self.compd[species]>(maxdens*0.01)][0] 126 127 #- e-10 as limit for minimum is ok, because if shell is 100000 R* 128 #- the mass conservation dictates ~ (10^5)^2 = 10^10 (r^2 law) 129 #- decrease in density. Shells this big dont occur anyway. 130 self.r_max[species] = self.rad[self.compd[species]>mindens][-1] 131 self.t_min[species] = temp[self.compd[species]>mindens][-1]

132 133 134

135 - def dustFullColDens(self,species):

136 137 """ 138 Calculate the full column density of a dust species in the shell. 139 140 This is NOT the value used for determining the abundance of the species 141 142 @param species: The dust species 143 @type species: string 144 145 @return: The column density of the dust species in the full envelope 146 in g/cm2 147 @rtype: float 148 149 """ 150 151 if not species in self.star.getDustList(): 152 return 0 153 return self.fullcoldens[species]

154 155

156 - def dustFullNumberColDens(self,species):

157 158 """ 159 Calculate the full NUMBER column density of a dust species in the shell. 160 161 This is NOT the number used in determining the equivalent molecular 162 abundance! 163 164 @param species: The dust species 165 @type species: string 166 167 @return: The number column density of the dust species in the full 168 envelope in cm-2 169 @rtype: float 170 171 """ 172 173 if not species in self.star.getDustList(): 174 return 0 175 if not self.star.dust[species]['molar']: 176 print 'No molar weight given for dust species %s in Dust.dat.'\ 177 %species 178 return 0 179 cndsp = self.fullcoldens[species]*self.avogadro\ 180 /self.star.dust[species]['molar'] 181 return cndsp

182 183 184

185 - def dustMolecAbun(self,species):

186 187 """ 188 Calculate the molecular abundance of a dust species with respect to H2. 189 190 Note that the self.coldens dictionary does not give real column 191 densities! This dict merely gives column densities in a prescribed 192 shell with given min and max radius, in order to compare with the H2 193 col density. 194 195 @param species: the dust species (from Dust.dat) 196 @type species: string 197 198 @return: The molecular abundance with respect to H2 of the dust species 199 @rtype: float 200 201 """ 202 203 if not species in self.star.getDustList(): 204 return 0 205 if not self.star.dust[species]['molar']: 206 print 'No molar weight given for dust species %s in Dust.dat.'\ 207 %species 208 return 0 209 if self.r_min_cd[species] == 0: 210 print 'No significant amount of dust species %s found.'%species 211 return 0 212 cndspecies = self.coldens[species]*self.avogadro\ 213 /self.star.dust[species]['molar'] 214 cndh2 = self.hydrogenColDens(species) 215 return cndspecies/cndh2

216 217 218

219 - def hydrogenColDens(self,species):

220 221 """ 222 Calculate the column number density of molecular hydrogen between two 223 radial distances. 224 225 If a GASTRoNOoM model is calculated, the h2 density profile is taken 226 from there. Otherwise, the value is derived from the total mass-loss 227 rate, assuming everything is H2. 228 229 The use of a proper velocity profile can make a big difference for low 230 mass-loss rate sources! 231 232 @param species: the dust species for which the H2 col dens is 233 calculated. Required for the radial information. 234 @type species: string 235 236 @return: The molecular hydrogen column number density between the given 237 radii (cm-2) 238 @rtype: float 239 240 """ 241 242 modelid = self.star['LAST_GASTRONOOM_MODEL'] 243 rin = self.r_min_cd[species] 244 rout = self.r_max_cd[species] 245 if modelid: 246 rad = self.star.getGasRad(ftype='fgr') 247 nh2 = self.star.getGasNumberDensity(ftype='fgr') 248 cndh2 = trapz(x=rad[(rad<rout)*(rad>rin)],\ 249 y=nh2[(rad<rout)*(rad>rin)]) 250 else: 251 mdot = (float(self.star['MDOT_GAS'])*u.M_sun/u.yr).to(u.g/u.s).value 252 vexp = float(self.star['VEL_INFINITY_GAS']) * 100000 253 h2_molar = 2. 254 sigma = (1./rin-1./rout)*mdot/vexp/4./math.pi 255 cndh2 = sigma * self.avogadro / h2_molar 256 return cndh2

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