Line data Source code
1 : /** \file stars.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Various utilities related to stars.
4 : * \ingroup astrofiles
5 : *
6 : */
7 :
8 : #include "../math/gslInterpolator.hpp"
9 :
10 : #include "units.hpp"
11 :
12 : #ifndef mx_astro_stars_hpp
13 : #define mx_astro_stars_hpp
14 :
15 : #include <cmath>
16 :
17 : #include "../ioutils/stringUtils.hpp"
18 : #include "../ioutils/readColumns.hpp"
19 :
20 : namespace mx
21 : {
22 : namespace astro
23 : {
24 :
25 : /// Parse a main sequence spectral type string into a numeric code.
26 : /** Expects standard spectral type strings such as "O5V" or "F2.5" or "G8.5V".
27 : * The numeric code starts at 0 for "O0V", 10 for "A0V", through 90 for "Y0". The subtypes are added to this value. For
28 : * instance, "G8.5V" becomes 48.5.
29 : *
30 : * Only works for main sequence types. Any other spectral type, such as MK class I-IV, will return -1.
31 : *
32 : * \ingroup stars
33 : */
34 : template <typename realT = float>
35 0 : realT numSpType( std::string spType /**< [in] The spectral type string to parse.*/ )
36 : {
37 :
38 0 : spType = ioutils::removeWhiteSpace( spType );
39 :
40 0 : spType = ioutils::toUpper( spType );
41 :
42 0 : realT num = 0;
43 :
44 0 : switch( spType[0] )
45 : {
46 0 : case 'O':
47 0 : num = 0;
48 0 : break;
49 0 : case 'B':
50 0 : num = 10;
51 0 : break;
52 0 : case 'A':
53 0 : num = 20;
54 0 : break;
55 0 : case 'F':
56 0 : num = 30;
57 0 : break;
58 0 : case 'G':
59 0 : num = 40;
60 0 : break;
61 0 : case 'K':
62 0 : num = 50;
63 0 : break;
64 0 : case 'M':
65 0 : num = 60;
66 0 : break;
67 0 : case 'L':
68 0 : num = 70;
69 0 : break;
70 0 : case 'T':
71 0 : num = 80;
72 0 : break;
73 0 : case 'Y':
74 0 : num = 90;
75 0 : break;
76 0 : default:
77 0 : return -1000;
78 : }
79 :
80 0 : if( spType.size() < 2 )
81 0 : return -1000;
82 0 : if( !isdigit( spType[1] ) )
83 0 : return -1000;
84 :
85 0 : int i = 1;
86 0 : while( i < spType.size() && ( isdigit( spType[i] ) || spType[i] == '.' ) )
87 0 : ++i;
88 :
89 0 : std::string subType = spType.substr( 1, i - 1 );
90 :
91 0 : num += ioutils::convertFromString<realT>( subType );
92 :
93 0 : if( spType.size() == i )
94 0 : return num;
95 0 : if( spType[i] == 'V' )
96 0 : return num;
97 :
98 0 : return -1 * num;
99 0 : }
100 :
101 : /// Provide various characteristics of main sequence stars according to their spectral type.
102 : /** Interpolates on "A Modern Mean Dwarf Stellar Color and Effective Temperature Sequence"
103 : * by Eric Mamajek, available at http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt.
104 : *
105 : * Loads values of the sequence at construction (hard coded), and then provides interpolation in between the points.
106 : *
107 : * Version of the sequence used: 2017.03.14
108 : *
109 : * \ingroup stars
110 : */
111 : template <typename realT>
112 : struct mainSequence
113 : {
114 : std::vector<double> m_numTypes; ///< The numeric spectral type codes from the sequence
115 : std::vector<double> m_numTypesR; ///< The numeric spectral type codes from the sequence, in reversed order.
116 :
117 : std::vector<double> m_Teffs; ///< The effective temperatures from the sequence
118 : std::vector<double> m_TeffsR; ///< The effective temperatures from the sequence, reveresed order.
119 :
120 : std::vector<double> m_logLs; ///< The log10 luminosities from the sequence
121 : std::vector<double> m_rads; ///< The radii from the sequence
122 : std::vector<double> m_masses; ///< The masses from the sequence
123 : std::vector<double> m_Mvs; ///< The absolute V magnitudes from the sequence
124 : std::vector<double> m_V_Rcs; ///< The V-Rc colors from the sequence
125 : std::vector<double> m_V_Ics; ///< The V-Ic colors from the sequence
126 : std::vector<double> m_V_Kss; ///< The V-Ks colors from the sequence
127 : std::vector<double> m_J_Hs; ///< The J-H colors from the sequence
128 : std::vector<double> m_H_Kss; ///< The H-Ks colors from the sequence
129 :
130 : math::gslInterpolator<math::gsl_interp_linear<double>> interpT; ///< The interpolator for effective Temperature
131 : double m_minT; ///< The minimum numeric spectral type for effective Temperature
132 : double m_maxT; ///< The maximum numeric spectral type for effective Temperature
133 :
134 : math::gslInterpolator<math::gsl_interp_linear<double>> interpRad; ///< The interpolator for effective Temperature
135 : double m_minRad; ///< The minimum numeric spectral type for effective Temperature
136 : double m_maxRad; ///< The maximum numeric spectral type for effective Temperature
137 :
138 : math::gslInterpolator<math::gsl_interp_linear<double>> interpL; ///< The interpolator for log luminosity
139 : double m_minL; ///< The minimum numeric spectral type for log luminosity
140 : double m_maxL; ///< The maximum numeric spectral type for log luminosity
141 :
142 : math::gslInterpolator<math::gsl_interp_linear<double>> interpMv; ///< The interpolator for absolute V magnitude
143 : double m_minMv; ///< The minimum numeric spectral type for absolute V magnitude
144 : double m_maxMv; ///< The maximum numeric spectral type for absolute V magnitude
145 :
146 : math::gslInterpolator<math::gsl_interp_linear<double>> interpVRc; ///< The interpolator for V-Rc color
147 : double m_minVRc; ///< The minimum numeric spectral type for V-Rc color
148 : double m_maxVRc; ///< The maximum numeric spectral type for V-Rc color
149 :
150 : math::gslInterpolator<math::gsl_interp_linear<double>> interpVIc; ///< The interpolator for V-Ic color
151 : double m_minVIc; ///< The minimum numeric spectral type for V-Ic color
152 : double m_maxVIc; ///< The maximum numeric spectral type for V-Ic color
153 :
154 : math::gslInterpolator<math::gsl_interp_linear<double>> interpVKs; ///< The interpolator for V-Ks color
155 : double m_minVKs; ///< The minimum numeric spectral type for V-Ks color
156 : double m_maxVKs; ///< The maximum numeric spectral type for V-Ks color
157 :
158 : math::gslInterpolator<math::gsl_interp_linear<double>> interpJH; ///< The interpolator for J-H color
159 : double m_minJH; ///< The minimum numeric spectral type for J-H color
160 : double m_maxJH; ///< The maximum numeric spectral type for J-H color
161 :
162 : math::gslInterpolator<math::gsl_interp_linear<double>> interpHKs; ///< The interpolator for H-Ks color
163 : double m_minHKs; ///< The minimum numeric spectral type for H-Ks color
164 : double m_maxHKs; ///< The maximum numeric spectral type for H-Ks color
165 :
166 : math::gslInterpolator<math::gsl_interp_linear<double>> interpSpTfmT; ///< The interpolator for effective Temperature
167 : double m_minSpTfmT; ///< The minimum numeric spectral type for effective Temperature
168 : double m_maxSpTfmT; ///< The maximum numeric spectral type for effective Temperature
169 :
170 0 : void findMinMax( std::vector<double> &seq, double &min, double &max, std::vector<double> &ref )
171 : {
172 0 : if( seq.size() != ref.size() )
173 : {
174 0 : std::cerr << "size mismatch\n";
175 0 : min = 0;
176 0 : max = 0;
177 : }
178 :
179 0 : min = ref[0];
180 : size_t n;
181 0 : for( n = 0; n < seq.size(); ++n )
182 : {
183 0 : if( seq[n] != -99 )
184 0 : break;
185 : }
186 :
187 0 : if( n >= seq.size() - 1 )
188 : {
189 0 : min = 0;
190 0 : max = 0;
191 0 : return;
192 : }
193 :
194 0 : min = ref[n];
195 :
196 0 : for( ; n < seq.size() - 1; ++n )
197 : {
198 0 : if( seq[n + 1] == -99 )
199 0 : break;
200 : }
201 :
202 0 : max = ref[n];
203 0 : ++n;
204 0 : for( ; n < seq.size(); ++n )
205 : {
206 0 : seq[n] = -99;
207 : }
208 : }
209 :
210 0 : mainSequence()
211 0 : {
212 : #include "mamajek.inc"
213 :
214 0 : m_numTypesR.assign( m_numTypes.rbegin(), m_numTypes.rend() );
215 :
216 0 : findMinMax( m_Teffs, m_minT, m_maxT, m_numTypes );
217 0 : m_TeffsR.assign( m_Teffs.rbegin(), m_Teffs.rend() ); // get reversed vector after normalizing -99s.
218 :
219 0 : findMinMax( m_rads, m_minRad, m_maxRad, m_numTypes );
220 :
221 0 : findMinMax( m_logLs, m_minL, m_maxL, m_numTypes );
222 0 : findMinMax( m_Mvs, m_minMv, m_maxMv, m_numTypes );
223 0 : findMinMax( m_V_Rcs, m_minVRc, m_maxVRc, m_numTypes );
224 0 : findMinMax( m_V_Ics, m_minVIc, m_maxVIc, m_numTypes );
225 0 : findMinMax( m_V_Kss, m_minVKs, m_maxVKs, m_numTypes );
226 0 : findMinMax( m_J_Hs, m_minJH, m_maxJH, m_numTypes );
227 0 : findMinMax( m_H_Kss, m_minHKs, m_maxHKs, m_numTypes );
228 :
229 0 : findMinMax( m_numTypesR, m_minSpTfmT, m_maxSpTfmT, m_TeffsR );
230 :
231 0 : interpT.setup( m_numTypes, m_Teffs );
232 0 : interpRad.setup( m_numTypes, m_rads );
233 0 : interpL.setup( m_numTypes, m_logLs );
234 0 : interpMv.setup( m_numTypes, m_Mvs );
235 0 : interpVRc.setup( m_numTypes, m_V_Rcs );
236 0 : interpVIc.setup( m_numTypes, m_V_Ics );
237 0 : interpVKs.setup( m_numTypes, m_V_Kss );
238 0 : interpJH.setup( m_numTypes, m_J_Hs );
239 0 : interpHKs.setup( m_numTypes, m_H_Kss );
240 :
241 0 : interpSpTfmT.setup( m_TeffsR, m_numTypesR );
242 0 : }
243 :
244 : /// Get the interpolated effective temperature
245 : /**
246 : * \returns the effective temperature.
247 : * \returns -999 if the input type code is out of range.
248 : */
249 0 : realT Teff( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
250 : {
251 0 : if( numType < m_minT || numType > m_maxT )
252 : {
253 0 : return -999;
254 : }
255 :
256 0 : return interpT( numType );
257 : }
258 :
259 : /// Get the interpolated effective temperature
260 : /**
261 : * \returns the effective temperature.
262 : * \returns -999 if the input type is out of range.
263 : */
264 0 : realT Teff( const std::string &spType /**< [in] the spectral type in standard format */ )
265 : {
266 0 : return Teff( numSpType( spType ) );
267 : }
268 :
269 : /// Get the interpolated radius
270 : /**
271 : * \returns the radius in Solar units.
272 : * \returns -999 if the input type code is out of range.
273 : */
274 0 : realT radius( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
275 : {
276 0 : if( numType < m_minRad || numType > m_maxRad )
277 : {
278 0 : return -999;
279 : }
280 :
281 0 : return interpRad( numType );
282 : }
283 :
284 : /// Get the interpolated radius
285 : /**
286 : * \returns the radius in Solar units.
287 : * \returns -999 if the input type is out of range.
288 : */
289 0 : realT radius( const std::string &spType /**< [in] the spectral type in standard format */ )
290 : {
291 0 : return radius( numSpType( spType ) );
292 : }
293 :
294 : /// Get the interpolated log of luminosity
295 : /**
296 : * \returns the log of luminosity.
297 : * \returns -999 if the input type code is out of range.
298 : */
299 0 : realT logL( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
300 : {
301 0 : if( numType < m_minL || numType > m_maxL )
302 : {
303 0 : return -999;
304 : }
305 :
306 0 : return interpL( numType );
307 : }
308 :
309 : /// Get the interpolated log of luminosity
310 : /**
311 : * \returns the log of luminosity.
312 : * \returns -999 if the input type is out of range.
313 : */
314 0 : realT logL( const std::string &spType /**< [in] the spectral type in standard format */ )
315 : {
316 0 : return logL( numSpType( spType ) );
317 : }
318 :
319 : /// Get the interpolated absolute V magnitude
320 : /**
321 : * \returns the aboslute V magnitude
322 : * \returns -999 if the input type code is out of range.
323 : */
324 0 : realT Mv( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
325 : {
326 0 : if( numType < m_minMv || numType > m_maxMv )
327 : {
328 0 : return -999;
329 : }
330 :
331 0 : return interpMv( numType );
332 : }
333 :
334 : /// Get the interpolated absolute V magnitude
335 : /**
336 : * \returns the aboslute V magnitude
337 : * \returns -999 if the input type is out of range.
338 : */
339 0 : realT Mv( const std::string &spType /**< [in] the spectral type in standard format */ )
340 : {
341 0 : return Mv( numSpType( spType ) );
342 : }
343 :
344 : /// Get the interpolated absolute V-Rc color
345 : /**
346 : * \returns the aboslute V-Rc color
347 : * \returns -999 if the input type code is out of range.
348 : */
349 0 : realT V_Rc( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
350 : {
351 0 : if( numType < m_minVRc || numType > m_maxVRc )
352 : {
353 0 : return -999;
354 : }
355 :
356 0 : return interpVRc( numType );
357 : }
358 :
359 : /// Get the interpolated absolute V-Rc color
360 : /**
361 : * \returns the aboslute V-Rc color
362 : * \returns -999 if the input type is out of range.
363 : */
364 0 : realT V_Rc( const std::string &spType /**< [in] the spectral type in standard format */ )
365 : {
366 0 : return V_Rc( numSpType( spType ) );
367 : }
368 :
369 : /// Get the interpolated absolute V-Ic color
370 : /**
371 : * \returns the aboslute V-Ic color
372 : * \returns -999 if the input type code is out of range.
373 : */
374 0 : realT V_Ic( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
375 : {
376 0 : if( numType < m_minVIc || numType > m_maxVIc )
377 : {
378 0 : return -999;
379 : }
380 :
381 0 : return interpVIc( numType );
382 : }
383 :
384 : /// Get the interpolated absolute V-Ic color
385 : /**
386 : * \returns the aboslute V-Ic color
387 : * \returns -999 if the input type is out of range.
388 : */
389 0 : realT V_Ic( const std::string &spType /**< [in] the spectral type in standard format */ )
390 : {
391 0 : return V_Ic( numSpType( spType ) );
392 : }
393 :
394 : /// Get the interpolated absolute V-Ks color
395 : /**
396 : * \returns the aboslute V-Ks color
397 : * \returns -999 if the input type code is out of range.
398 : */
399 0 : realT V_Ks( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
400 : {
401 0 : if( numType < m_minVKs || numType > m_maxVKs )
402 : {
403 0 : return -999;
404 : }
405 :
406 0 : return interpVKs( numType );
407 : }
408 :
409 : /// Get the interpolated absolute V-Ks color
410 : /**
411 : * \returns the aboslute V-Ks color
412 : * \returns -999 if the input type is out of range.
413 : */
414 0 : realT V_Ks( const std::string &spType /**< [in] the spectral type in standard format */ )
415 : {
416 0 : return V_Ks( numSpType( spType ) );
417 : }
418 :
419 : /// Get the interpolated absolute J-H color
420 : /**
421 : * \returns the aboslute J-H color
422 : * \returns -999 if the input type code is out of range.
423 : */
424 0 : realT J_H( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
425 : {
426 0 : if( numType < m_minJH || numType > m_maxJH )
427 : {
428 0 : return -999;
429 : }
430 :
431 0 : return interpJH( numType );
432 : }
433 :
434 : /// Get the interpolated absolute J-H color
435 : /**
436 : * \returns the aboslute J-H color
437 : * \returns -999 if the input type is out of range.
438 : */
439 0 : realT J_H( const std::string &spType /**< [in] the spectral type in standard format */ )
440 : {
441 0 : return J_H( numSpType( spType ) );
442 : }
443 :
444 : /// Get the interpolated absolute H-Ks color
445 : /**
446 : * \returns the aboslute H-Ks color
447 : * \returns -999 if the input type code is out of range.
448 : */
449 0 : realT H_Ks( realT numType /**< [in] the numeric spectral type code, see \ref mx::astro::numSpType() */ )
450 : {
451 0 : if( numType < m_minHKs || numType > m_maxHKs )
452 : {
453 0 : return -999;
454 : }
455 :
456 0 : return interpHKs( numType );
457 : }
458 :
459 : /// Get the interpolated absolute H-Ks color
460 : /**
461 : * \returns the aboslute H-Ks color
462 : * \returns -999 if the input type is out of range.
463 : */
464 0 : realT H_Ks( const std::string &spType /**< [in] the spectral type in standard format */ )
465 : {
466 0 : return H_Ks( numSpType( spType ) ); /// Get the interpolated absolute H-Ks color
467 : /**
468 : * \returns the aboslute H-Ks color
469 : * \returns -999 if the input type code is out of range.
470 : */
471 : }
472 :
473 0 : realT spTFromTeff( realT Teff )
474 : {
475 0 : if( Teff < m_minSpTfmT || Teff > m_maxSpTfmT )
476 : {
477 0 : return -999;
478 : }
479 :
480 0 : return std::floor( 2 * interpSpTfmT( Teff ) + 0.5 ) / 2; // round to nearest 0.5 types.
481 : }
482 : };
483 :
484 : namespace maintenance
485 : {
486 :
487 : /// Read in the main sequence table of Mamajek, and construct the vectors for input into the mainSequence class.
488 : /** This is used to create mamajek.inc.
489 : * The table should be copied to a text file, and all `...` replaced with -99, then any remaining . replaced with space.
490 : */
491 0 : void makeMSTable( const std::string &fn )
492 : {
493 0 : std::vector<std::string> SpT;
494 0 : std::vector<double> Teff;
495 0 : std::vector<double> logT;
496 0 : std::vector<double> logL;
497 0 : std::vector<double> Mbol;
498 0 : std::vector<double> BCv;
499 0 : std::vector<double> Mv;
500 0 : std::vector<double> B__V;
501 0 : std::vector<double> Bt__Vt;
502 0 : std::vector<double> G__V;
503 0 : std::vector<double> Bp__Rp;
504 0 : std::vector<double> G__Rp;
505 0 : std::vector<double> M_G;
506 0 : std::vector<double> b__y;
507 0 : std::vector<double> U__B;
508 0 : std::vector<double> V__Rc;
509 0 : std::vector<double> V__Ic;
510 0 : std::vector<double> V__Ks;
511 0 : std::vector<double> J__H;
512 0 : std::vector<double> H__Ks;
513 0 : std::vector<double> Ks__W1;
514 0 : std::vector<double> W1__W2;
515 0 : std::vector<double> W1__W3;
516 0 : std::vector<double> W1__W4;
517 0 : std::vector<double> M_J;
518 0 : std::vector<double> M_Ks;
519 0 : std::vector<double> i__z;
520 0 : std::vector<double> z__Y;
521 0 : std::vector<double> R_Rsun;
522 0 : std::vector<double> Msun;
523 :
524 0 : ioutils::readColumns( fn,
525 : SpT,
526 : Teff,
527 : logT,
528 : logL,
529 : Mbol,
530 : BCv,
531 : Mv,
532 : B__V,
533 : Bt__Vt,
534 : G__V,
535 : Bp__Rp,
536 : G__Rp,
537 : M_G,
538 : b__y,
539 : U__B,
540 : V__Rc,
541 : V__Ic,
542 : V__Ks,
543 : J__H,
544 : H__Ks,
545 : Ks__W1,
546 : W1__W2,
547 : W1__W3,
548 : W1__W4,
549 : M_J,
550 : M_Ks,
551 : i__z,
552 : z__Y,
553 : R_Rsun,
554 : Msun );
555 :
556 0 : std::cout << "//" << fn << "\n";
557 0 : std::cout << " m_numTypes = {" << numSpType( SpT[0] );
558 0 : for( int n = 1; n < SpT.size(); ++n )
559 0 : std::cout << ',' << numSpType( SpT[n] );
560 0 : std::cout << "};\n";
561 :
562 0 : std::cout << " m_Teffs = {" << Teff[0];
563 0 : for( int n = 1; n < SpT.size(); ++n )
564 0 : std::cout << ',' << Teff[n];
565 0 : std::cout << "};\n";
566 :
567 0 : std::cout << " m_logLs = {" << logL[0];
568 0 : for( int n = 1; n < SpT.size(); ++n )
569 0 : std::cout << ',' << logL[n];
570 0 : std::cout << "};\n";
571 :
572 0 : std::cout << " m_rads = {" << R_Rsun[0];
573 0 : for( int n = 1; n < SpT.size(); ++n )
574 0 : std::cout << ',' << R_Rsun[n];
575 0 : std::cout << "};\n";
576 :
577 0 : std::cout << " m_masses = {" << Msun[0];
578 0 : for( int n = 1; n < SpT.size(); ++n )
579 0 : std::cout << ',' << Msun[n];
580 0 : std::cout << "};\n";
581 :
582 0 : std::cout << " m_Mvs = {" << Mv[0];
583 0 : for( int n = 1; n < SpT.size(); ++n )
584 0 : std::cout << ',' << Mv[n];
585 0 : std::cout << "};\n";
586 :
587 0 : std::cout << " m_V_Rcs = {" << V__Rc[0];
588 0 : for( int n = 1; n < SpT.size(); ++n )
589 0 : std::cout << ',' << V__Rc[n];
590 0 : std::cout << "};\n";
591 :
592 0 : std::cout << " m_V_Ics = {" << V__Ic[0];
593 0 : for( int n = 1; n < SpT.size(); ++n )
594 0 : std::cout << ',' << V__Ic[n];
595 0 : std::cout << "};\n";
596 :
597 0 : std::cout << " m_V_Kss = {" << V__Ks[0];
598 0 : for( int n = 1; n < SpT.size(); ++n )
599 0 : std::cout << ',' << V__Ks[n];
600 0 : std::cout << "};\n";
601 :
602 0 : std::cout << " m_J_Hs = {" << J__H[0];
603 0 : for( int n = 1; n < SpT.size(); ++n )
604 0 : std::cout << ',' << J__H[n];
605 0 : std::cout << "};\n";
606 :
607 0 : std::cout << " m_H_Kss = {" << H__Ks[0];
608 0 : for( int n = 1; n < SpT.size(); ++n )
609 0 : std::cout << ',' << H__Ks[n];
610 0 : std::cout << "};\n";
611 0 : }
612 : } // namespace maintenance
613 : } // namespace astro
614 : } // namespace mx
615 :
616 : #endif
|
