Line data Source code
1 : /** \file astroSpectrum.hpp
2 : * \author Jared R. Males
3 : * \brief A class for working with astronomical spectra.
4 : * \ingroup astrophot
5 : *
6 : */
7 :
8 : #ifndef mx_astro_astroSpectrum_hpp
9 : #define mx_astro_astroSpectrum_hpp
10 :
11 : #include <cmath>
12 :
13 : #include "../sys/environment.hpp"
14 : #include "../ioutils/readColumns.hpp"
15 : #include "../math/gslInterpolation.hpp"
16 : #include "constants.hpp"
17 : #include "units.hpp"
18 :
19 : namespace mx
20 : {
21 : namespace astro
22 : {
23 :
24 : /// Base spectrum class which provides manipulation and characterization functionality.
25 : /**
26 : * \ingroup astrophot
27 : */
28 : template <typename realT>
29 : struct baseSpectrum
30 : {
31 : std::vector<realT> _spectrum; ///< Contains the spectrum after it is set.
32 :
33 : /// Get the current size of the spectrum
34 : /**
35 : * \returns the number of points in the spectrum
36 : */
37 0 : size_t size()
38 : {
39 0 : return _spectrum.size();
40 : }
41 :
42 : /// Access a single point in the spectrum, specified by its vector index.
43 : /**
44 : * \returns a reference to the value at i.
45 : */
46 0 : realT &operator[]( size_t i /**< [in] the index of the spectral point*/ )
47 : {
48 0 : return _spectrum[i];
49 : }
50 :
51 : /// Access a single point in the spectrum, specified by its vector index.
52 : /**
53 : * \returns the value of point i.
54 : */
55 0 : const realT operator[]( size_t i /**< [in] the index of the spectral point*/ ) const
56 : {
57 0 : return _spectrum[i];
58 : }
59 :
60 : /// Multiply two spectra together.
61 : /**
62 : * \returns a new baseSpectrum
63 : */
64 : template <typename compSpectrumT>
65 : baseSpectrum<realT> operator*( const compSpectrumT &spec /**< [in] the spectrum to multiply by */ )
66 : {
67 : baseSpectrum outVec;
68 : outVec._spectrum.resize( _spectrum.size() );
69 :
70 : for( int i = 0; i < _spectrum.size(); ++i )
71 : {
72 : outVec[i] = _spectrum[i] * spec[i];
73 : }
74 :
75 : return outVec;
76 : }
77 :
78 : /// Calculate the mean value of the spectrum.
79 : /**
80 : * \returns the mean value of the spectrum
81 : */
82 0 : realT mean()
83 : {
84 0 : realT sum = 0;
85 :
86 0 : for( int i = 0; i < _spectrum.size(); ++i )
87 0 : sum += _spectrum[i];
88 :
89 0 : return sum / _spectrum.size();
90 : }
91 :
92 : /// Calculate the mean value of the spectrum when mutiplied by another.
93 : /** The result is normalized by the mean value of the input spectrum, equivalent to:
94 : * \f[
95 : \mu = \frac{\int T(\lambda)S(\lambda) d\lambda}{\int T(\lambda) d\lambda}
96 : \f]
97 : * For instance, use this to get the mean value of a spectrum in a filter.
98 : *
99 : * \returns the mean value of the multiplied spectrum.
100 : *
101 : * \tparam compSpectrumT the vector-like type of the comparison spectrum.
102 : */
103 : template <typename compSpectrumT>
104 : realT mean( const compSpectrumT &T /**< [in] the spectrum to multiply by*/ )
105 : {
106 : realT sum1 = 0, sum2 = 0;
107 :
108 : for( int i = 0; i < _spectrum.size(); ++i )
109 : {
110 : sum1 += _spectrum[i] * T[i];
111 : sum2 += T[i];
112 : }
113 :
114 : return sum1 / sum2;
115 : }
116 :
117 : /// Characterize the spectrum as a filter transmission curve.
118 : /** For a photonic transmission curve given by \f$ S(\lambda ) \f$ The mean photon wavelength is defined as
119 : * \f[
120 : \lambda_0 = \frac{1}{\Delta\lambda_{0}}\int \frac{S(\lambda )}{S_{max}} \lambda d\lambda
121 : \f]
122 : * which is Equation A14 of Bessel 2012.
123 : *
124 : * where the effective width is defined by
125 : \f[
126 : \Delta\lambda_{o} = \int \frac{S(\lambda )} { S_{max}} d\lambda
127 : \f]
128 : *
129 : * The full-width at half-maximum, FWHM, is the distance between the points at 50% of maximum \f$ S(\lambda) \f$.
130 : *
131 : */
132 0 : void charTrans( realT &lambda0, ///< [out] the central wavelength of the filter
133 : realT &weff, ///< [out] the effective width of the filter
134 : realT &max, ///< [out] the maximum value of the transmission curve
135 : realT &fwhm, ///< [out] the full-width at half-maximum of the filter profile
136 : std::vector<realT> &lambda ///< [in] the wavelength scale, should correspond to the spectrum.
137 : )
138 : {
139 0 : size_t N = lambda.size();
140 0 : realT dl = lambda[1] - lambda[0];
141 0 : realT half = static_cast<realT>( 0.5 );
142 :
143 0 : max = 0;
144 0 : for( int i = 0; i < N; ++i )
145 : {
146 0 : if( _spectrum[i] > max )
147 0 : max = _spectrum[i];
148 : }
149 :
150 0 : weff = half * _spectrum[0];
151 0 : for( int i = 1; i < N - 1; ++i )
152 0 : weff += _spectrum[i];
153 0 : weff += half * _spectrum[N - 1];
154 0 : weff *= dl / max;
155 :
156 0 : lambda0 = half * lambda[0] * _spectrum[0];
157 0 : for( int i = 1; i < N - 1; ++i )
158 : {
159 0 : lambda0 += lambda[i] * _spectrum[i];
160 : }
161 0 : lambda0 += half * lambda[N - 1] * _spectrum[N - 1];
162 0 : lambda0 *= dl / max / weff;
163 :
164 : realT left, right;
165 :
166 0 : int i = 1;
167 0 : while( i < lambda.size() )
168 : {
169 0 : if( _spectrum[i] >= 0.5 * max )
170 0 : break;
171 0 : ++i;
172 : }
173 :
174 : // Interpolate
175 0 : left = lambda[i - 1] + ( 0.5 * max - _spectrum[i - 1] ) * ( dl / ( _spectrum[i] - _spectrum[i - 1] ) );
176 :
177 0 : i = N - 2;
178 0 : while( _spectrum[i] < 0.5 * max )
179 : {
180 0 : --i;
181 0 : if( i < 0 )
182 0 : break;
183 : }
184 :
185 0 : right = lambda[i] + ( _spectrum[i] - 0.5 * max ) * ( dl / ( _spectrum[i] - _spectrum[i + 1] ) );
186 :
187 0 : fwhm = right - left;
188 0 : }
189 :
190 : /// Characterize the flux densities of the spectrum w.r.t. a filter transmission curve
191 : /** To obtain the flux (e.g. W/m^2) multiply these quantities by the effective width calculated using
192 : * \ref charTrans.
193 : *
194 : * This implements Equations A11, A12, and A13 of Bessel 2012.
195 : *
196 : * \warning this only produces correct fphot0 for a spectrum in W/m^3. DO NOT USE FOR ANYTHING ELSE.
197 : *
198 : * \todo use unit conversions to make it work for everything.
199 : * \todo check on integration method, should it be trap?
200 : *
201 : */
202 : template <class filterT>
203 : void charFlux( realT &flambda0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in W/m^3
204 : realT &fnu0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in W/m^2/Hz
205 : realT &fphot0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in photons/sec/m^3
206 : const realT &lambda_0, ///< [in] the mean photon wavelength lambda_0 (from charTrans).
207 : const std::vector<realT> &lambda, ///< [in] the wavelength scale of this spectrum.
208 : const filterT &trans ///< [in] the filter transmission curve over which to characterize, on the same
209 : ///< wavelength grid.
210 : )
211 : {
212 : charFlux( flambda0, fnu0, fphot0, lambda_0, lambda, trans._spectrum );
213 : }
214 :
215 : /// Characterize the flux densities of the spectrum w.r.t. a filter transmission curve
216 : /** To obtain the flux (e.g. W/m^2) multiply these quantities by the effective width calculated using
217 : * \ref charTrans.
218 : *
219 : * This implements Equations A11, A12, and A13 of Bessel 2012.
220 : *
221 : * \warning this only produces correct fphot0 for a spectrum in W/m^3. DO NOT USE FOR ANYTHING ELSE.
222 : *
223 : * \todo use unit conversions to make it work for everything.
224 : * \todo check on integration method, should it be trap?
225 : *
226 : */
227 0 : void charFlux( realT &flambda0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in W/m^3
228 : realT &fnu0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in W/m^2/Hz
229 : realT &fphot0, ///< [out] the flux of the star at \f$ \lambda_0 \f$ in photons/sec/m^3
230 : const realT &lambda_0, ///< [in] the mean photon wavelength lambda_0 (from charTrans).
231 : const std::vector<realT> &lambda, ///< [in] the wavelength scale of this spectrum.
232 : const std::vector<realT> &trans ///< [in] the filter transmission curve over which to characterize,
233 : ///< on the same wavelength grid.
234 : )
235 : {
236 0 : constexpr realT h = constants::h<units::si<realT>>();
237 0 : constexpr realT c = constants::c<units::si<realT>>();
238 :
239 0 : size_t N = lambda.size();
240 0 : realT dl = lambda[1] - lambda[0];
241 0 : realT half = static_cast<realT>( 0.5 );
242 :
243 : // flambda0 -- Eqn A11
244 :
245 0 : realT denom = half * trans[0] * lambda[0];
246 0 : for( int i = 1; i < N - 1; ++i )
247 : {
248 0 : denom += trans[i] * lambda[i];
249 : }
250 0 : denom += half * trans[N - 1] * lambda[N - 1];
251 0 : denom *= dl;
252 :
253 0 : realT num = half * _spectrum[0] * trans[0] * lambda[0];
254 0 : for( int i = 1; i < N - 1; ++i )
255 : {
256 0 : num += _spectrum[i] * trans[i] * lambda[i];
257 : }
258 0 : num += half * _spectrum[N - 1] * trans[N - 1] * lambda[N - 1];
259 0 : num *= dl;
260 :
261 0 : flambda0 = num / denom;
262 :
263 : // fnu0 -- Eqn A12
264 :
265 0 : denom = half * trans[0] / lambda[0];
266 0 : for( int i = 1; i < N - 1; ++i )
267 : {
268 0 : denom += trans[i] / lambda[i];
269 : }
270 0 : denom += half * trans[N - 1] / lambda[N - 1];
271 0 : denom *= dl * c;
272 :
273 0 : fnu0 = num / denom;
274 :
275 : // fphot0 -- Eqn A13
276 :
277 0 : fphot0 = flambda0 * lambda_0 / ( h * c );
278 0 : }
279 : };
280 :
281 : /// Class to manage an astronomical spectrum.
282 : /** Details of the spectra, including units, file-location, and how they are read from disk, are specified by template
283 : * parameter _spectrumT: see \ref astrophot_spectra "the available spectrum definitions." Spectra are loaded and
284 : * immediately interpolated onto a wavelength grid. This is to facilitate manipulation of spectra, i.e. for multiplying
285 : * by filters, etc.
286 : *
287 : * Inherits functions and operators to manipulate and characterize spectra from \ref mx::astro::baseSpectrum
288 : * "baseSpectrum".
289 : *
290 : * \tparam _spectrumT is the underlying spectrum type, which provides (through static interfaces) the specifications of
291 : * how to read or calculate the spectrum.
292 : * \tparam freq specify whether this spectrum uses frequency instead of
293 : * wavelength. Default is false.
294 : *
295 : * \ingroup astrophot
296 : */
297 : template <typename _spectrumT, bool freq = false, typename verboseT = verbose::d>
298 : struct astroSpectrum : public baseSpectrum<typename _spectrumT::units::realT>
299 : {
300 : typedef _spectrumT spectrumT;
301 : typedef typename spectrumT::units units;
302 : typedef typename units::realT realT;
303 : typedef typename spectrumT::paramsT paramsT;
304 :
305 : std::string _dataDir; ///< The directory containing the spectrum
306 :
307 : paramsT _params; ///< The parameters of the spectrum, e.g. its name or the numerical parameters needed to generate
308 : ///< the name.
309 :
310 : /// Default c'tor
311 0 : astroSpectrum()
312 0 : {
313 0 : }
314 :
315 : /// Constructor specifying name, the enivronment will be queried for data directory.
316 0 : explicit astroSpectrum( const paramsT ¶ms /**< [in] The name of the spectrum */ )
317 0 : {
318 0 : setParameters( params );
319 0 : }
320 :
321 : /// Constructor specifying name and data directory.
322 0 : astroSpectrum( const paramsT ¶ms, ///< [in] The name of the spectrum
323 : const std::string &dataDir ///< [in] The directory containing the spectrum
324 : )
325 0 : {
326 0 : setParameters( params, dataDir );
327 0 : }
328 :
329 : /// Set the parameters of the spectrum, using the underlying spectrums parameter type.
330 0 : int setParameters( const paramsT ¶ms /**< [in] The name of the spectrum */ )
331 : {
332 0 : _params = params;
333 :
334 : if( spectrumT::dataDirEnvVar )
335 : {
336 0 : _dataDir = sys::getEnv( spectrumT::dataDirEnvVar );
337 : }
338 :
339 0 : return 0;
340 : }
341 :
342 : /// Set the parameters of the spectrum, using the underlying spectrum's parameter type.
343 : /** This version also sets the data directory, instead of using the enivronment variable.
344 : *
345 : * \overload
346 : */
347 0 : int setParameters( const paramsT ¶ms, ///< [in] The name of the spectrum
348 : const std::string &dataDir ///< [in] The directory containing the spectrum
349 : )
350 : {
351 0 : _params = params;
352 0 : _dataDir = dataDir;
353 :
354 0 : return 0;
355 : }
356 :
357 : /// Load the spectrum and interpolate it onto a wavelength scale.
358 : template <typename gridT>
359 0 : error_t setSpectrum( gridT &lambda /**< [in] the wavelength scale on which to interpolate.
360 : Same units as in source file. */
361 : )
362 : {
363 0 : std::vector<realT> rawLambda;
364 0 : std::vector<realT> rawSpectrum;
365 :
366 0 : std::string fileName = spectrumT::fileName( _params );
367 :
368 0 : std::string path;
369 :
370 0 : if( fileName.size() < 1 )
371 : {
372 0 : return internal::mxlib_error_report<verboseT>( error_t::paramnotset, "filename is empty" );
373 : }
374 :
375 0 : if( _dataDir == "" && fileName[0] == '/' )
376 : {
377 0 : path = fileName;
378 : }
379 : else
380 : {
381 0 : if( _dataDir == "" )
382 0 : _dataDir = ".";
383 :
384 0 : path = _dataDir + "/" + fileName;
385 : }
386 :
387 0 : error_t rv = spectrumT::readSpectrum( rawLambda, rawSpectrum, path, _params );
388 0 : if( rv != error_t::noerror )
389 : {
390 0 : return internal::mxlib_error_report<verboseT>( rv );
391 : }
392 :
393 : // Unit conversions
394 0 : for( int i = 0; i < rawLambda.size(); ++i )
395 : {
396 0 : rawLambda[i] /= spectrumT::wavelengthUnits / ( units::length );
397 0 : rawSpectrum[i] /= spectrumT::fluxUnits /
398 : ( units::energy / ( units::time * units::length * units::length * units::length ) );
399 : }
400 :
401 0 : math::gsl_interpolate( ::gsl_interp_linear, rawLambda, rawSpectrum, lambda, this->_spectrum );
402 :
403 0 : for( int i = 0; i < lambda.size(); ++i )
404 : {
405 0 : if( !std::isnormal( this->_spectrum[i] ) )
406 0 : this->_spectrum[i] = 0;
407 : }
408 :
409 0 : return error_t::noerror;
410 0 : }
411 : };
412 :
413 : } // namespace astro
414 :
415 : } // namespace mx
416 :
417 : #endif // mx_astro_astroSpectrum_hpp
|
