1 #include "specclass.h"
2 #include "strings.h"
3 #include "mathdefs.h"
4 #include <map>
5
6 using std::string;
7
8 // These are only supposed to handle numbers in the range [0, 6].
int_to_roman(double MKtype)9 static string int_to_roman(double MKtype)
10 {
11 string result;
12 switch (FLOOR(MKtype)) {
13 case 0: result = "0"; case 1: result = "I"; case 2: result = "II";
14 case 3: result = "III"; case 4: result = "IV"; case 5: result = "V";
15 case 6: result = "VI";
16 default: result = "";
17 }
18 return result;
19 }
20
roman_to_int(const string & roman_num)21 static double roman_to_int(const string &roman_num)
22 {
23 if (starstrings::compare_n(roman_num, "VI", 2)) return 6.0;
24 else if (starstrings::compare_n(roman_num, "V", 1)) return 5.0;
25 else if (starstrings::compare_n(roman_num, "III", 3)) return 3.0;
26 else if (starstrings::compare_n(roman_num, "II", 2)) return 2.0;
27 else if (starstrings::compare_n(roman_num, "IV", 2)) return 4.0;
28 else if (starstrings::compare_n(roman_num, "Ia", 2)) return 0.5;
29 else if (starstrings::compare_n(roman_num, "Ib", 2)) return 1.5;
30 else if (starstrings::compare_n(roman_num, "I", 1)) return 1.0;
31 else return -1;
32 }
33
class_to_int(char major)34 int SpecClass::class_to_int(char major)
35 {
36 static bool initialized = false;
37 static std::map<char, int> table;
38 major = starstrings::toupper(major);
39 if (! initialized) {
40 table['O'] = table['W'] = 0;
41 table['B'] = 1;
42 table['A'] = 2;
43 table['F'] = 3;
44 table['G'] = 4;
45 table['K'] = table['R'] = 5;
46 table['M'] = table['N'] = table['S'] = table['C'] = 6;
47 initialized = true;
48 }
49 int temp = table[major];
50 if ((! temp) && major != 'O' && major != 'W') return -1;
51 else return temp;
52 }
53
54
55 // Minor and MK components of spectral type are not determined unless
56 // specifically requested by initialize() function, below.
SpecClass(const string & s)57 SpecClass::SpecClass(const string &s)
58 : sSpecstring(s), sMajor(0), sMinor(ERROR_VAL), sMKtype(ERROR_VAL),
59 sSpecial("")
60 {
61 if (! starstrings::isempty(sSpecstring))
62 sMajor = starstrings::toupper(s[s.find_first_not_of(WHITESPACE "ds")]);
63 }
64
SpecClass(const SpecClass & sc)65 SpecClass::SpecClass(const SpecClass &sc)
66 : sSpecstring(sc.sSpecstring), sMajor(sc.sMajor), sMinor(sc.sMinor),
67 sMKtype(sc.sMKtype), sSpecial(sc.sSpecial)
68 { }
69
SpecClass(char major,double minor,double MKtype,string special)70 SpecClass::SpecClass(char major, double minor, double MKtype, string special)
71 : sMajor(major), sMinor(minor), sMKtype(MKtype), sSpecial(special)
72 {
73 sSpecstring = string() + major + starstrings::ftoa(minor, 2) + " "
74 + int_to_roman(MKtype) + sSpecial;
75 }
76
77
operator =(const SpecClass & sc)78 SpecClass & SpecClass::operator = (const SpecClass &sc)
79 {
80 if (this != &sc) {
81 sSpecstring = sc.sSpecstring; sMajor = sc.sMajor;
82 sMinor = sc.sMinor; sMKtype = sc.sMKtype; sSpecial = sc.sSpecial;
83 }
84 return *this;
85 }
86
87
color() const88 color_t SpecClass::color() const
89 {
90 switch (sMajor) {
91 case 'W' : return WOLF_RAYET_COLOR;
92 case 'D' : return D_COLOR;
93 case '*' : return NON_STELLAR_COLOR;
94 default : break;
95 }
96 int classno = class_to_int(sMajor);
97 if (classno >= 0) return CLASS_COLORS[classno];
98 else return DEFAULT_COLOR;
99 }
100
initialize()101 void SpecClass::initialize()
102 {
103 size_t begin, end;
104 starstrings::stripspace(sSpecstring);
105 starstrings::utf8ize(sSpecstring);
106
107 // calculate minor part of spectral type
108 begin = sSpecstring.find_first_of(DIGITS);
109 end = sSpecstring.find_first_not_of(DIGITS ".", begin);
110 sMinor = (begin != end) ?
111 starmath::atof(sSpecstring.substr(begin, end - begin)) : 10;
112 if (end >= sSpecstring.size()) return;
113
114 // calculate luminosity class
115 begin = sSpecstring.find_first_of("IV", end);
116 end = sSpecstring.find_first_not_of("IVab", begin);
117 bool flagsplit = false;
118
119 if (begin >= sSpecstring.size() || begin == end) return;
120 sMKtype = roman_to_int(sSpecstring.substr(begin, end - begin));
121 if (sSpecstring[end] == '/' || sSpecstring[end] == '-') {
122 flagsplit = true;
123 end++;
124 }
125 if (flagsplit && sMKtype >= 0) sMKtype += 0.5;
126
127 if (end < sSpecstring.size()) {
128 sSpecial = sSpecstring.substr(end);
129 starstrings::stripspace(sSpecial);
130 }
131 }
132
133
134 // Lookup tables of data for use by absmag, diameter, temperature functions.
135 // These tables contain linearly interpolated data derived from Appendix E of
136 // Carroll and Ostlie, _An Introduction to Modern Astrophysics_ (1996).
137 // Data are in the form { value(minor == 0), d(value)/d(minor) } for each
138 // spectral class, at MK types V, III, and Iab respectively.
139
140 // Effective temperatures, in Kelvins (column 2 of tables in Appendix E)
141 static const double temp_table[3][7][2] = {
142 {
143 { 59000, -3000 }, { 30000, -2048 }, { 9520, -232 }, { 7200, -117 },
144 { 6030, -78 }, { 5250, -140 }, { 3850, -151 }
145 },
146 {
147 { 56000, -2700 }, { 29000, -1890 }, { 10100, -295 }, { 7150, -130 },
148 { 5850, -110 }, { 4750, -95 }, { 3800, -93 }
149 },
150 {
151 { 54600, -2860 }, { 26000, -1627 }, { 9730, -203 }, { 7700, -215 },
152 { 5550, -113 }, { 4420, -77 }, { 3650, -175 }
153 }
154 };
155
156 // Absolute visual magnitudes, M_v (column 8 of tables)
157 static const double mag_table[3][7][2] = {
158 {
159 { -7.4, 0.34 }, { -4.0, 0.46 }, { 0.6, 0.21 }, { 2.7, 0.17 },
160 { 4.4, 0.15 }, { 5.9, 0.29 }, { 8.8, 0.9 }
161 },
162 {
163 { -7.5, 0.24 }, { -5.1, 0.51 }, { 0.0, 0.15 }, { 1.5, -0.05 },
164 { 1.0, -0.03 }, { 0.7, -0.11 }, { -0.4, 0.03 }
165 },
166 {
167 { -6.8, 0.02 }, { -6.4, 0.01 }, { -6.3, -0.03 }, { -6.6, 0.02 },
168 { -6.4, 0.04 }, { -6.0, 0.04 }, { -5.6, 0.0 }
169 }
170 };
171
172 // Bolometric correction (column 7 of tables), BC == M_bolometric - M_v
173 static const double bolcorr_table[3][7][2] = {
174 {
175 { -5.64, 0.248 }, { -3.16, 0.286 }, { -0.30, 0.021 }, { -0.09, -0.009 },
176 { -0.18, -0.013 }, { -0.31, -0.107 }, { -1.38, -0.34 }
177 },
178 {
179 { -5.22, 0.234 }, { -2.88, 0.246 }, { -0.42, 0.031 }, { -0.11, -0.009 },
180 { -0.20, -0.03 }, { -0.50, -0.075 }, { -1.25, -0.247 }
181 },
182 {
183 { -5.25, 0.276 }, { -2.49, 0.208 }, { -0.41, 0.040 }, { -0.01, -0.014 },
184 { -0.15, -0.035 }, { -0.50, -0.079 }, { -1.29, -0.435 }
185 }
186 };
187
188
189 // lookup(): This function looks up values for the current spectral class
190 // in the given table. If MK type is between array elements, it interpolates.
191 // If no appropriate value is found, it returns a very large negative number
192 // to mark the error.
193
194 #define LOOKUP_RESULT(i) \
195 (table[(5 - (i)) / 2][class_to_int(sMajor)][0] + \
196 table[(5 - (i)) / 2][class_to_int(sMajor)][1] * sMinor)
197
lookup(const double table[3][7][2]) const198 double SpecClass::lookup(const double table[3][7][2]) const
199 {
200 int major_num = class_to_int(sMajor);
201 if (major_num < 0 || major_num > 6) return ERROR_VAL;
202
203 if (fabs(sMKtype - 5.0) < 0.1)
204 return LOOKUP_RESULT(5);
205 else if (fabs(sMKtype - 3.0) < 0.1)
206 return LOOKUP_RESULT(3);
207 else if (fabs(sMKtype - 1.0) < 0.1)
208 return LOOKUP_RESULT(1);
209
210 // if MKtype isn't exactly I, III or V, interpolate linearly between them
211 else if (sMKtype >= 0.0 && sMKtype < 3.0)
212 return ((3 - sMKtype) * LOOKUP_RESULT(1) +
213 (sMKtype - 1) * LOOKUP_RESULT(3)) / 2;
214 else if (sMKtype <= 6.0 && sMKtype > 3.0)
215 return ((5 - sMKtype) * LOOKUP_RESULT(3) +
216 (sMKtype - 3) * LOOKUP_RESULT(5)) / 2;
217 else
218 return ERROR_VAL;
219 }
220
221
222 // absmag(): returns a guess for the absolute magnitude of a star with
223 // the current spectral class. Returns a large negative number if an
224 // error occurs.
absmag() const225 double SpecClass::absmag() const
226 {
227 if (sMKtype < 0.0) {
228 // in this case, no MK type is available;
229 // try to save a few special cases from being thrown out
230 if (sMajor == 'D') return +14;
231 else return ERROR_VAL;
232 }
233 return lookup(mag_table);
234 }
235
236
237 // temperature(): returns the expected temperature in Kelvins
238 // of a star of this spectral class.
temperature() const239 double SpecClass::temperature() const
240 {
241 double temp = lookup(temp_table);
242 return (temp <= 0.0) ? ERROR_VAL : temp;
243 }
244
245
246 // diameter(): return the expected diameter, in light-years, of a star of
247 // this spectral class with the given absolute magnitude. We require the
248 // bolometric magnitude, M_bolo == BC + M_v, in order to get the correct
249 // diameter.
diameter(double absmag,bool bolometric) const250 double SpecClass::diameter(double absmag, bool bolometric) const
251 {
252 if (sMajor == 'D') return +1.5e-9; // white dwarfs ~ diameter of Earth
253
254 double temp = temperature();
255 double bolcorr = bolometric ? 0.0 : lookup(bolcorr_table);
256
257 if (temp <= 0.0 || bolcorr < -10.0 || absmag < -25.0)
258 return ERROR_VAL;
259
260 // Now use the Stefan-Boltzmann law:
261 // R* = RSun * sqrt(L* / LSun) * (TSun / T*)^2
262 // but L* / LSun = 10^(0.4 * (MSun - M*))
263 // so R* = RSun * 10^(0.2 * (MSun - M*)) * (TSun / T*)^2
264
265 using std::pow;
266 return SUN_DIAMETER * pow(SUN_TEMP / temp, 2) *
267 pow(10.0, (SUN_BOLO_MAG - absmag - bolcorr) / 5);
268 }
269
270