ComboCode.cc.modeling.profilers.Temperature

1 # -*- coding: utf-8 -*- 2 3 """ 4 Module for calculating the temperature profile as a function of radius. 5 6 Author: R. Lombaert 7 8 """ 9 10 import numpy as np 11 import collections 12 from scipy.interpolate import InterpolatedUnivariateSpline as spline1d 13 14 from cc.data import Data 15 from cc.modeling.profilers import Profiler 16 from cc.tools.numerical import Operators as op 17 18 19

20 -def Teps(r,T0,r0,epsilon):

21 22 ''' 23 Define a temperature power law. 24 25 The functional form is: 26 Tg(r) = T0 * (R0/r)^epsilon 27 28 @param r: The radial grid in cm 29 @type r: array 30 @param T0: The initial temperature at Ri in K 31 @type T0: float 32 @param r0: The inner radius in cm 33 @type r0: float 34 @param epsilon: The exponent of the temp law 35 @type epsilon: float 36 37 @return: The temperature grid (K) 38 @rtype: array 39 40 ''' 41 42 return T0 * (r/r0)**-epsilon

43 44 45

46 -def Tdust(r,T0,r0,s=1):

47 48 ''' 49 Return a dust temperature power law of the form as suggested by 50 observational evidence. 51 52 See Thesis p32, where power is -2/(4+s) in 53 T(r) = T_eff*(2*r/R_STAR)**(-2/(4+s)) 54 and s is typically 1 in the Rayleigh limit and optically thin case, with 55 T_eff = T0, R_STAR = r0. 56 57 @param r: The radial grid for the power law in Rstar 58 @type r: array 59 @param T0: The stellar effective temperature 60 @type T0: float 61 @param r0: The initial value for the radius, typically the stellar radius. 62 @type r0 63 64 @keyword s: The s parameter in the power law T(r) given above. 65 66 (default: 1) 67 @type s: int 68 69 @return: The temperature grid (K) 70 @rtype: array 71 72 ''' 73 74 return T0*(2.*r/r0)**(-2./(4+s))

75 76 77

78 -class Temperature(Profiler.Profiler):

79 80 ''' 81 An interface for a temperature profile and its derivative 82 83 ''' 84

85 - def __init__(self,r,func,dfunc=None,order=3,inner=0,inner_eps=0.5,\ 86 *args,**kwargs):

87 88 ''' 89 Create an instance of the Temperature() class. Requires a radial grid 90 and a temperature function. 91 92 The function can also be given as an interp1d object. 93 94 The optional args and kwargs give the additional arguments for the 95 temperature function, which are ignored in case func is an interp1d 96 object. 97 98 An additional option concerns the extrapolation to smaller radii, e.g. 99 in the inner wind between the stellar surface and the dust condensation 100 radius. In this case, the eval/diff methods will differ between inner 101 wind (r<r0) and the rest of the wind (r>r0) in terms of evaluation. 102 103 At r>r0: The given func (or interpolator) is used. 104 At r<r0: A Teps power law is used, for which 1 out of Tstar or 105 epsilon can be defined. The default sets epsilon == 0.5 (the optically 106 thin case), but inner_epsilon can be passed to the initialisation if 107 needed. r0, T0 must be defined in kwargs either for the 108 main wind's Teps law, or as an additional argument if a file is read. 109 110 r0 and T0 are usually passed as kwargs for the func. THey are used here 111 explicitly for the inner wind power law. If they are not given, the 112 defaults are used: r[0] and -- after T has been evaluated -- T[0]. In 113 case the func is interp_file, r0 and T0 can still be passed as kwargs 114 but they are then removed from the kwargs dict after extraction. 115 116 @param r: The radial points (cm) 117 @type r: array 118 @param func: The function that describes the temperature profile. Can be 119 given as an interp1d object. 120 @type func: function 121 122 @keyword dfunc: The function that describes the derivative of the 123 profile with respect to r. Can be given as an interp1d 124 object. If None, a generic central difference is taken 125 and interpolated. 126 127 (default: None) 128 @type dfunc: function/interp1d object 129 @keyword order: Order of the spline interpolation. Default is cubic. 130 131 (default: 3) 132 @type order: int 133 @keyword inner: Applies a power law to the inner region, as a means of 134 extrapolation. Off by default, but can be turned on. 135 In this case, r0 is taken from kwargs, ie the r0 from 136 the power law of the wind (if Teps is chosen), or as an 137 additional keyword if a file is read. 138 139 (default: 0) 140 @type inner: bool 141 @keyword inner_eps: The exponent for the power law in the inner wind. 142 143 (default: 0.5) 144 @type inner_eps: float 145 146 @keyword args: Additional parameters passed to the functions when eval 147 or diff are called. 148 149 (default: []) 150 @type args: tuple 151 @keyword kwargs: Additional keywords passed to the functions when eval 152 or diff are called. 153 154 (default: {}) 155 @type kwargs: dict 156 157 ''' 158 159 #-- Check if a file is requested to be interpolated: use r0 and T0 from 160 # the file 161 if func == 'interp_file' or func == Profiler.interp_file: 162 self.r0 = None 163 self.T0 = None 164 else: 165 self.r0 = kwargs.get('r0') 166 self.T0 = kwargs.get('T0') 167 168 #-- Do not name func, dfunc, etc in function call, or *args won't work 169 super(Temperature, self).__init__(r,func,dfunc,order,*args,**kwargs) 170 171 self.inner = inner 172 173 #-- No extrapolation to the inner wind requested 174 if not self.inner: 175 self.r = self.x 176 self.T = self.y 177 self.dTdr = self.dydx 178 return 179 180 #-- Check if r0 and T0 are set, if not: interpolation of a file was 181 # requested, so take first values from file 182 if self.r0 is None: self.r0 = self.xin[0] 183 if self.T0 is None: self.T0 = self.yin[0] 184 185 #-- Alternatively, a power law extrapolation is requested 186 self.setInnerEps(inner_eps=inner_eps)

187 188 189

190 - def setInnerEps(self,inner_eps):

191 192 ''' 193 Set the inner wind epsilon for the T law. Only relevant when self.inner 194 is 1. 195 196 Changes the standard T law as well, in case inner eps is different from 197 when the instance of the class was created. 198 199 @param inner_eps: the epsilon of the inner wind power law. 200 @type inner_eps: float 201 202 ''' 203 204 self.inner_eps = inner_eps 205 206 #-- re-initialise some parameters 207 self.r = None 208 self.T = None 209 self.y = None 210 self.dTdr = None 211 212 #-- Run the class' own eval() method, with r=self.x, which will not be 213 # equal to the class' self.r attribute (which is now None), and thus 214 # it will evaluated properly. Reevaluated if epsilon is different upon 215 # eval call. 216 self.T = self.eval(self.x,inner_eps=self.inner_eps,warn=0) 217 self.y = self.T 218 219 #-- Now also redo the derivative method. Then call diff to set the 220 # standard derivative of the instance 221 if self.interp_dfunc: 222 self.dfunc = spline1d(x=self.x,y=op.diff_central(self.y,self.x),\ 223 k=self.order) 224 225 #-- Evaluate the derivative with the default grid 226 self.dydx = self.diff(self.x,inner_eps=self.inner_eps,warn=0) 227 228 #-- Now set r and dTdr. 229 self.r = self.x 230 self.dTdr = self.dydx

231 232 233

234 - def __call__(self,r=None,warn=1,*args,**kwargs):

235 236 ''' 237 Evaluate the profile function at a coordinate point. 238 239 r can be any value or array. 240 241 The default y-grid is returned if r is None. 242 243 Passes the call to eval. Additional keywords for the power law 244 extrapolation can be passed this way as well. 245 246 @keyword r: The primary coordinate point(s). If None, the default 247 coordinate grid is used. 248 249 (default: None) 250 @type r: array/float 251 @keyword warn: Warn when extrapolation occurs. 252 253 (default: 1) 254 @type warn: bool 255 256 @return: The profile evaluated at r 257 @rtype: array/float 258 259 ''' 260 261 return self.eval(r,warn=warn,*args,**kwargs)

262 263 264

265 - def eval(self,r=None,warn=1,inner_eps=None):

266 267 ''' 268 Evaluate the profile function at a coordinate point. 269 270 r can be any value or array. 271 272 The default y-grid is returned if r is None and inner_eps is 0.5 273 (the default). 274 275 @keyword r: The coordinate point(s). If None, the default 276 coordinate grid is used. 277 278 (default: None) 279 @type r: array/float 280 @keyword warn: Warn when extrapolation occurs. 281 282 (default: 1) 283 @type warn: bool 284 @keyword inner_eps: The exponent of the Teps power law inner wind 285 extrapolation. If default, the inner_eps defined 286 upon initialisation is used. 287 288 (default: None) 289 @type inner_eps: float 290 291 @return: The profile evaluated at r 292 @rtype: array/float 293 294 ''' 295 296 #-- First retrieve the original profile. If this is the first call from 297 # the __init__ method of Temperature(), this will be the standard grid 298 y = super(Temperature,self).eval(r,warn=warn) 299 300 #-- No inner power law requested, just pass on the original 301 if not self.inner: 302 return y 303 304 #-- Return self.y since r was given as None, if the eps is correct 305 if inner_eps is None: inner_eps = self.inner_eps 306 if r is None and inner_eps == self.inner_eps: 307 return self.y 308 309 #-- r can still be None. Apparently different inner_eps requested. 310 # So calc the profile anew with the inner wind law. Need r defined. 311 if r is None: 312 r = self.r 313 314 #-- Replace the extrapolated values in the inner wind with the new power 315 # law. Make sure r is an array for this. 316 rarr = Data.arrayify(r) 317 y[rarr<self.r0] = Teps(rarr[rarr<self.r0],T0=self.T0,r0=self.r0,\ 318 epsilon=inner_eps) 319 320 return y if isinstance(r,collections.Iterable) else y[0]

321 322 323

324 - def diff(self,r=None,warn=1,inner_eps=None):

325 326 ''' 327 Evaluate the derivative of the profile function at a coordinate point. 328 329 r can be any value or array. 330 331 The default y-grid is returned if r is None and inner_eps is 0.5 332 (the default). 333 334 @keyword r: The coordinate point(s). If None, the default 335 coordinate grid is used. 336 337 (default: None) 338 @type r: array/float 339 @keyword warn: Warn when extrapolation occurs. 340 341 (default: 1) 342 @type warn: bool 343 @keyword inner_eps: The exponent of the Teps power law inner wind 344 extrapolation. If default, the inner_eps defined 345 upon initialisation is used. 346 347 (default: None) 348 @type inner_eps: float 349 350 @return: The derivative evaluated at r 351 @rtype: array/float 352 353 ''' 354 355 #-- First retrieve the original profile. If this is the first call from 356 # the __init__ method of Opacity(), this will be the standard grid. 357 dydx = super(Temperature,self).diff(r,warn=warn) 358 359 #-- No inner power law requested, just pass on the original 360 if not self.inner: 361 return dydx 362 363 #-- Return self.dydr since r was given as None, if the eps is correct 364 if inner_eps is None: inner_eps = self.inner_eps 365 if r is None and inner_eps == self.inner_eps: 366 return self.y 367 368 #-- r can still be None. Apparently different inner_eps requested. 369 # So calc the profile anew with the inner wind law. Need r defined. 370 if r is None: 371 r = self.r 372 373 #-- Replace the extrapolated values in the inner wind with the new power 374 # law. Make sure r is an array for this. 375 rarr = Data.arrayify(r) 376 dy_fac = Teps(rarr[rarr<self.r0],T0=self.T0,r0=self.r0,\ 377 epsilon=inner_eps-1) 378 dydx[rarr<self.r0] = -dy_fac*inner_eps/self.r0 379 380 return dydx if isinstance(r,collections.Iterable) else dydx[0]

381

Source Code for Module ComboCode.cc.modeling.profilers.Temperature