1
2
3 """
4 Module for calculating the spectral radiance as a function of wavelength.
5
6 Provides the alternative expressions such as flux density, luminosity, etc as
7 well.
8
9 Note that I avoid the term specific or total intensity, as there is confusion
10 between the precise definition of the term. Spectral intensity refers to
11 radiance integrated over a surface, ie in units of erg/s/Hz/sr. Hence, I
12 consistently refer to radiance when specific intensity is concerned.
13
14 Also note that radiance is also known as the brightness. That term however
15 refers to the subjective brightness of an object, and is not very appropriate in
16 this context.
17
18 For an overview of terminology:
19 https://en.wikipedia.org/wiki/Radiance
20
21 For an overview on relation between brightness (or specific intensity), flux
22 density, and luminosity:
23 http://www.cv.nrao.edu/course/astr534/Brightness.html
24
25 Author: R. Lombaert
26
27 """
28
29 import numpy as np
30 from astropy import constants as cst
31
32 from cc.modeling.profilers import Profiler
33
34
35
37
38 '''
39 Calculate a blackbody spectrum (ie brightness) as a function of
40 frequency and for a given temperature.
41
42 Gives the energy per unit time per unit area of emitting surface in the
43 normal direction per unit solid angle per unit frequency.
44
45 Note that GASTRoNOoM works in brightness! Not flux density!
46
47 @param f: The frequency grid (cm)
48 @type f: array/float
49 @param T: The temperature value (not array, K)
50 @type T: float
51
52 @return: The blackbody spectrum in erg/s/cm2/Hz
53 @rtype: array/float
54
55 '''
56
57
58 h = cst.h.cgs.value
59 c = cst.c.cgs.value
60 k = cst.k_B.cgs.value
61
62 return 2*h*f**3/c**2/(np.exp((h*f)/(k*T))-1.)
63
64
65
67
68 '''
69 An interface for a spectral radiance and its derivative.
70
71 '''
72
73 - def __init__(self,func,f=None,l=None,dfunc=None,order=3,*args,**kwargs):
74
75 '''
76 Create an instance of the Radiance() class. Requires a frequency grid
77 and a spectral radiance function.
78
79 The function can also be given as an interpolation object.
80
81 The optional args and kwargs give the additional arguments for the
82 temperature function, which are ignored in case func is an interp1d
83 object.
84
85 @param func: The function that describes the intensity profile. Can be
86 given as an interp1d object.
87 @type func: function
88
89 @keyword f: The frequency points (Hz). This or wavelength array must
90 be given!
91
92 (default: None)
93 @type f: array
94 @keyword l: The wavelength points (cm). This or frequency array must
95 be given!
96
97 (default: None)
98 @type l: array
99 @keyword dfunc: The function that describes the derivative of the
100 profile with respect to r. Can be given as an interp
101 object. If None, a generic central difference is taken
102 and interpolated.
103
104 (default: None)
105 @type dfunc: function/interpolation object
106 @keyword order: Order of the spline interpolation. Default is cubic.
107
108 (default: 3)
109 @type order: int
110
111 @keyword args: Additional parameters passed to the functions when eval
112 or diff are called.
113
114 (default: [])
115 @type args: tuple
116 @keyword kwargs: Additional keywords passed to the functions when eval
117 or diff are called.
118
119 (default: {})
120 @type kwargs: dict
121
122 '''
123
124 if not (l is None or f is None):
125 raise ValueError('Both wavelength and frequency given. ' + \
126 'Define only one.')
127
128
129 self.c = cst.c.cgs.value
130 if f is None:
131 f = self.c/l[::-1]
132
133
134 super(Radiance, self).__init__(f,func,dfunc,order,*args,**kwargs)
135
136
137
138 self.f = self.x
139 self.R = self.y
140 self.dRdl = self.dydx
141
142
143 self.Lsun = cst.L_sun.cgs.value
144 self.Rsun = cst.R_sun.cgs.value
145
146
147
149
150 '''
151 Define the surface area of the star.
152
153 @param rstar: The stellar radius (cm)
154 @type rstar: float
155
156 '''
157
158 self.surface = 4*np.pi*rstar**2
159
160
161
162 - def eval(self,f=None,l=None,ftype='Rnu',warn=1):
163
164 '''
165 Evaluate the radiance function at a coordinate point.
166
167 l/f can be any value or array. If func is an interpolation object, it is
168 in principle limited by the l/f-range of the interpolator. However,
169 extrapolation is enabled, but it is advised not to extend much
170 beyond the given l/f-range.
171
172 Can be evaluated both in frequency and wavelength (l or f), and with
173 respect to both frequency and wavelength (ftype).
174
175 @keyword f: The frequency grid (Hz). Passed to self.eval(). Can be None
176 in which case the default grid is used. Can also be given as
177 wavelength
178
179 (default: None)
180 @type f: array/float
181 @keyword l: The wavelength grid (cm). Converted to frequency before
182 evaluation.
183
184 (default: None)
185 @type l: array/float
186 @keyword ftype: Return the spectral radiance with respect to wavelength
187 or frequency. Given as '*nu' or '*lambda', where * is
188 R, I, F, ... Default is always '*nu', also
189 when ftype is not recognized. Rlambda is calculated from
190 Rlambda = Rnu*nu/lambda.
191
192 (default: 'Rnu')
193 @type ftype: str
194 @keyword warn: Warn when extrapolation occurs.
195
196 (default: 1)
197 @type warn: bool
198
199 @return: The profile evaluated at l/f (erg/s/cm2/Hz/sr or erg/s/cm3/sr)
200 @rtype: array/float
201
202 '''
203
204
205 if not (l is None or f is None):
206 msg = 'Both wavelength and frequency given. Define only one.'
207 raise ValueError(msg)
208
209
210
211
212 if not l is None:
213 f = self.c/l
214
215 else:
216 l = self.c/self.f if f is None else self.c/f
217
218
219
220 radiance = super(Radiance,self).eval(x=f,warn=warn)
221
222
223 if ftype.lower()[1:] == 'lambda':
224 if f is None: f = self.f
225 radiance = radiance*f/l
226
227 return radiance
228
229
230
231 - def diff(self,f=None,l=None,warn=1):
232
233 '''
234 Evaluate the derivative of the radiance function at a coordinate point.
235
236 For now only possible with respect to frequency.
237
238 x can be any value or array. If func is an interpolation object, it is
239 in principle limited by the x-range of the interpolator. However,
240 linear extrapolation is enabled, but it is advised not to extend much
241 beyond the given x-range.
242
243 @keyword f: The frequency grid (Hz). Passed to self.eval(). Can be None
244 in which case the default grid is used. Can also be given as
245 wavelength
246
247 (default: None)
248 @type f: array/float
249 @keyword l: The wavelength grid (cm). Converted to frequency before
250 evaluation.
251
252 (default: None)
253 @type l: array/float
254 @keyword warn: Warn when extrapolation occurs.
255
256 (default: 1)
257 @type warn: bool
258
259 @return: The derivative evaluated at l or f
260 @rtype: array/float
261
262 '''
263
264
265 if not (l is None or f is None):
266 msg = 'Both wavelength and frequency given. Define only one.'
267 raise ValueError(msg)
268
269
270
271
272 if f is None:
273 f = self.c/l
274
275
276
277 return super(Radiance,self).diff(x=f,warn=warn)
278
279
280
281 - def getLuminosity(self,surface=None,f=None,l=None,ftype='Lnu',warn=1):
282
283 '''
284 Return the spectral luminosity, i.e. integrated over the surface
285 and solid angle, as a function of frequency or wavelength.
286
287 Integrated over frequency, this leads to the (total) luminosity, or
288 also the "flux". Do not integrate over wavelength! Units not right for
289 that.
290
291 Note that the integration over solid angle leads to the factor of pi.
292 Fnu = Int(cos(theta) dOmega) = Int Int (cos(theta)sin(theta)dtheta dpsi)
293 which has primitive function sin^2(t)/, psi in [0,2pi], theta in
294 [0,pi/2]. Note that this is the outgoing flux, along the direction of
295 the ray (opposite the direction would be [pi/2,pi], since theta is wrt
296 the normal of the surface).
297
298 Can be evaluated both in frequency and wavelength (l or f), and with
299 respect to both frequency and wavelength (ftype).
300
301 @keyword surface: The surface over which the energy is being emitted
302 (cm^2). If None, self.setSurface must have been called
303 before. self.surface gives the stellar surface area.
304
305 (default: None)
306 @type surface: float
307 @keyword f: The frequency grid (Hz). Passed to self.eval(). Can be None
308 in which case the default grid is used. Can also be given as
309 wavelength
310
311 (default: None)
312 @type f: array/float
313 @keyword l: The wavelength grid (cm). Converted to frequency before
314 evaluation.
315
316 (default: None)
317 @type l: array/float
318 @keyword ftype: Return spectral luminosity with respect to wavelength
319 or frequency. Given as 'Lnu' or 'Llambda'. Default is
320 always 'Lnu', also when ftype is not recognized. Llambda
321 is calculated from Llambda = Lnu*nu/lambda.
322
323 (default: 'Lnu')
324 @type ftype: str
325 @keyword warn: Warn when extrapolation occurs.
326
327 (default: 1)
328 @type warn: bool
329
330 @return: The spectral luminosity (erg/s/Hz or erg/s/cm)
331 @rtype: array
332
333 '''
334
335
336 radiance = self.eval(l=l,f=f,ftype=ftype,warn=warn)
337 if surface is None: surface = self.surface
338 return np.pi*radiance*surface
339
340
341
343
344 '''
345 Return the spectral flux density, i.e. integrated over solid
346 angle, as a function of wavelength.
347
348 Integrated over wavelength, this leads to the (total) flux density.
349 For integration over wavelength, convert to Flambda = Fnu*nu/lambda
350 first.
351
352 Note that the integration over solid angle leads to the factor of pi.
353 Fnu = Int(cos(theta) dOmega) = Int Int (cos(theta)sin(theta)dtheta dpsi)
354 which has primitive function sin^2(t)/2, psi in [0,2pi], theta in
355 [0,pi/2]. Note that this is the outgoing flux, along the direction of
356 the ray (opposite the direction would be [pi/2,pi], since theta is wrt
357 the normal of the surface).
358
359 Can be evaluated both in frequency and wavelength (l or f), and with
360 respect to both frequency and wavelength (ftype).
361
362 @keyword f: The frequency grid (Hz). Passed to self.eval(). Can be None
363 in which case the default grid is used. Can also be given as
364 wavelength
365
366 (default: None)
367 @type f: array/float
368 @keyword l: The wavelength grid (cm). Converted to frequency before
369 evaluation.
370
371 (default: None)
372 @type l: array/float
373 @keyword ftype: Return the spectral flux dens with respect to wavelength
374 or frequency. Given as 'Fnu' or 'Flambda'. Default is
375 always 'Fnu', also when ftype is not recognized. Flambda
376 is calculated from Flambda = Fnu*nu/lambda.
377
378 (default: 'Fnu')
379 @type ftype: str
380 @keyword warn: Warn when extrapolation occurs.
381
382 (default: 1)
383 @type warn: bool
384
385 @return: The spectral flux density (erg/s/cm2/Hz or erg/s/cm3)
386 @rtype: array
387
388 '''
389
390
391 radiance = self.eval(l=l,f=f,ftype=ftype,warn=warn)
392 return np.pi*radiance
393
394
395
396 - def getIntensity(self,surface=None,f=None,l=None,ftype='Inu',warn=1):
397
398 '''
399 Return the spectral intensity, i.e. integrated over the surface, as a
400 function of wavelength or frequency.
401
402 Integrated over frequency, this leads to the (total) intensity. Do not
403 integrate over wavelength! Units not right for that.
404
405 Not to be confused with the specific intensity, which is the radiance or
406 brightness of the source!
407
408 Can be evaluated both in frequency and wavelength (l or f), and with
409 respect to both frequency and wavelength (ftype).
410
411 @keyword surface: The surface over which the energy is being emitted
412 (cm^2). If None, self.setSurface must have been called
413 before. self.surface gives the stellar surface area.
414
415 (default: None)
416 @type surface: float
417 @keyword f: The frequency grid (Hz). Passed to self.eval(). Can be None
418 in which case the default grid is used. Can also be given as
419 wavelength
420
421 (default: None)
422 @type f: array/float
423 @keyword l: The wavelength grid (cm). Converted to frequency before
424 evaluation.
425
426 (default: None)
427 @type l: array/float
428 @keyword ftype: Return the spectral intensity with respect to wavelength
429 or frequency. Given as 'Inu' or 'Ilambda'. Default is
430 always 'Inu', also when ftype is not recognized. Ilambda
431 is calculated from Ilambda = Inu*nu/lambda.
432
433 (default: 'Inu')
434 @type ftype: str
435 @keyword warn: Warn when extrapolation occurs.
436
437 (default: 1)
438 @type warn: bool
439
440 @return: The spectral intensity (erg/s/Hz/sr or erg/s/cm/sr)
441 @rtype: array
442
443 '''
444
445 radiance = self.eval(l=l,f=f,ftype=ftype,warn=warn)
446 if surface is None: surface = self.surface
447 return radiance*surface
448