1 // star.cpp
2 //
3 // Copyright (C) 2001, Chris Laurel <claurel@shatters.net>
4 //
5 // This program is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU General Public License
7 // as published by the Free Software Foundation; either version 2
8 // of the License, or (at your option) any later version.
9 
10 #include <celmath/mathlib.h>
11 #include <cstring>
12 #include <cassert>
13 #include <cstdio>
14 #include "celestia.h"
15 #include "astro.h"
16 #include "orbit.h"
17 #include "star.h"
18 #include "texmanager.h"
19 
20 using namespace std;
21 
22 
23 // The value of the temperature of the sun is actually 5780, but the
24 // stellar class tables list the temperature of a G2V star as 5860.  We
25 // use the latter value so that the radius of the sun is computed correctly
26 // as one times SOLAR_RADIUS . . .  the high metallicity of the Sun is
27 // probably what accounts for the discrepancy in temperature.
28 // #define SOLAR_TEMPERATURE    5780.0f
29 #define SOLAR_TEMPERATURE    5860.0f
30 #define SOLAR_RADIUS         696000
31 
32 
33 struct SpectralTypeInfo
34 {
35     char* name;
36     float temperature;
37     float rotationPeriod;
38 };
39 
40 
41 static StarDetails** normalStarDetails = NULL;
42 static StarDetails** whiteDwarfDetails = NULL;
43 static StarDetails*  neutronStarDetails = NULL;
44 static StarDetails*  blackHoleDetails = NULL;
45 static StarDetails*  barycenterDetails = NULL;
46 
47 static string DEFAULT_INFO_URL("");
48 
49 StarDetails::StarTextureSet StarDetails::starTextures;
50 
51 // Star temperature data from Lang's _Astrophysical Data: Planets and Stars_
52 // Temperatures from missing (and typically not used) types in those
53 // tables were just interpolated.
54 static float tempO[3][10] =
55 {
56     { 52500, 52500, 52500, 52500, 48000, 44500, 41000, 38000, 35800, 33000 },
57     { 50000, 50000, 50000, 50000, 45500, 42500, 39500, 37000, 34700, 32000 },
58     { 47300, 47300, 47300, 47300, 44100, 42500, 39500, 37000, 34700, 32000 },
59 };
60 
61 static float tempB[3][10] =
62 {
63     { 30000, 25400, 22000, 18700, 17000, 15400, 14000, 13000, 11900, 10500 },
64     { 29000, 24000, 20300, 17100, 16000, 15000, 14100, 13200, 12400, 11000 },
65     { 26000, 20800, 18500, 16200, 15100, 13600, 13000, 12200, 11200, 10300 },
66 };
67 
68 static float tempA[3][10] =
69 {
70     {  9520, 9230, 8970, 8720, 8460, 8200, 8020, 7850, 7580, 7390 },
71     { 10100, 9480, 9000, 8600, 8300, 8100, 7850, 7650, 7450, 7250 },
72     {  9730, 9230, 9080, 8770, 8610, 8510, 8310, 8150, 7950, 7800 },
73 };
74 
75 static float tempF[3][10] =
76 {
77     { 7200, 7050, 6890, 6740, 6590, 6440, 6360, 6280, 6200, 6110 },
78     { 7150, 7000, 6870, 6720, 6570, 6470, 6350, 6250, 6150, 6080 },
79     { 7700, 7500, 7350, 7150, 7000, 6900, 6500, 6300, 6100, 5800 },
80 };
81 
82 static float tempG[3][10] =
83 {
84     { 6030, 5940, 5860, 5830, 5800, 5770, 5700, 5630, 5570, 5410 },
85     { 5850, 5650, 5450, 5350, 5250, 5150, 5050, 5070, 4900, 4820 },
86     { 5550, 5350, 5200, 5050, 4950, 4850, 4750, 4660, 4600, 4500 },
87 };
88 
89 static float tempK[3][10] =
90 {
91     { 5250, 5080, 4900, 4730, 4590, 4350, 4200, 4060, 3990, 3920 },
92     { 4750, 4600, 4420, 4200, 4000, 3950, 3900, 3850, 3830, 3810 },
93     { 4420, 4330, 4250, 4080, 3950, 3850, 3760, 3700, 3680, 3660 },
94 };
95 
96 static float tempM[3][10] =
97 {
98     { 3850, 3720, 3580, 3470, 3370, 3240, 3050, 2940, 2640, 2000 },
99     { 3800, 3720, 3620, 3530, 3430, 3330, 3240, 3240, 3240, 3240 },
100     { 3650, 3550, 3450, 3200, 2980, 2800, 2600, 2600, 2600, 2600 },
101 };
102 
103 // Wolf-Rayet temperatures.  From Lang's Astrophysical Data: Planets and
104 // Stars.
105 static float tempWN[10] =
106 {
107     50000, 50000, 50000, 50000, 47000, 43000, 39000, 32000, 29000, 29000
108 };
109 
110 static float tempWC[10] =
111 {
112     60000, 60000, 60000, 60000, 60000, 60000, 60000, 54000, 46000, 38000
113 };
114 
115 // Brown dwarf temperatures
116 static float tempL[10] =
117 {
118     1960, 1930, 1900, 1850, 1800, 1740, 1680, 1620, 1560, 1500
119 };
120 
121 static float tempT[10] =
122 {
123     1425, 1350, 1275, 1200, 1140, 1080, 1020, 900, 800, 750
124 };
125 
126 // White dwarf temperaturs
127 static float tempWD[10] =
128 {
129     100000.0f, 50400.0f, 25200.0f, 16800.0f, 12600.0f,
130     10080.0f, 8400.0f, 7200.0f, 6300.0f, 5600.0f,
131 };
132 
133 
134 // Tables with adjustments for estimating absolute bolometric magnitude from
135 // visual magnitude, from Lang's "Astrophysical Data: Planets and Stars".
136 // Gaps in the tables from unused spectral classes were filled in with linear
137 // interpolation--not accurate, but these shouldn't appear in real catalog
138 // data anyway.
139 static float bmag_correctionO[3][10] =
140 {
141     // Lum class V (main sequence)
142     {
143         -4.75f, -4.75f, -4.75f, -4.75f, -4.45f,
144         -4.40f, -3.93f, -3.68f, -3.54f, -3.33f,
145     },
146     // Lum class III
147     {
148         -4.58f, -4.58f, -4.58f, -4.58f, -4.28f,
149         -4.05f, -3.80f, -3.58f, -3.39f, -3.13f,
150     },
151     // Lum class I
152     {
153         -4.41f, -4.41f, -4.41f, -4.41f, -4.17f,
154         -3.87f, -3.74f, -3.48f, -3.35f, -3.18f,
155     }
156 };
157 
158 static float bmag_correctionB[3][10] =
159 {
160     // Lum class V (main sequence)
161     {
162         -3.16f, -2.70f, -2.35f, -1.94f, -1.70f,
163         -1.46f, -1.21f, -1.02f, -0.80f, -0.51f,
164     },
165     // Lum class III
166     {
167         -2.88f, -2.43f, -2.02f, -1.60f, -1.45f,
168         -1.30f, -1.13f, -0.97f, -0.82f, -0.71f,
169     },
170     // Lum class I
171     {
172         -2.49f, -1.87f, -1.58f, -1.26f, -1.11f,
173         -0.95f, -0.88f, -0.78f, -0.66f, -0.52f,
174     }
175 };
176 
177 static float bmag_correctionA[3][10] =
178 {
179     // Lum class V (main sequence)
180     {
181         -0.30f, -0.23f, -0.20f, -0.17f, -0.16f,
182         -0.15f, -0.13f, -0.12f, -0.10f, -0.09f,
183     },
184     // Lum class III
185     {
186         -0.42f, -0.29f, -0.20f, -0.17f, -0.15f,
187         -0.14f, -0.12f, -0.10f, -0.10f, -0.10f,
188     },
189     // Lum class I
190     {
191         -0.41f, -0.32f, -0.28f, -0.21f, -0.17f,
192         -0.13f, -0.09f, -0.06f, -0.03f, -0.02f,
193     }
194 };
195 
196 static float bmag_correctionF[3][10] =
197 {
198     // Lum class V (main sequence)
199     {
200         -0.09f, -0.10f, -0.11f, -0.12f, -0.13f,
201         -0.14f, -0.14f, -0.15f, -0.16f, -0.17f,
202     },
203     // Lum class III
204     {
205         -0.11f, -0.11f, -0.11f, -0.12f, -0.13f,
206         -0.13f, -0.15f, -0.15f, -0.16f, -0.18f,
207     },
208     // Lum class I
209     {
210         -0.01f,  0.00f,  0.00f, -0.01f, -0.02f,
211         -0.03f, -0.05f, -0.07f, -0.09f, -0.12f,
212     }
213 };
214 
215 static float bmag_correctionG[3][10] =
216 {
217     // Lum class V (main sequence)
218     {
219         -0.18f, -0.19f, -0.20f, -0.20f, -0.21f,
220         -0.21f, -0.27f, -0.33f, -0.40f, -0.36f,
221     },
222     // Lum class III
223     {
224         -0.20f, -0.24f, -0.27f, -0.29f, -0.32f,
225         -0.34f, -0.37f, -0.40f, -0.42f, -0.46f,
226     },
227     // Lum class I
228     {
229         -0.15f, -0.18f, -0.21f, -0.25f, -0.29f,
230         -0.33f, -0.36f, -0.39f, -0.42f, -0.46f,
231     }
232 };
233 
234 static float bmag_correctionK[3][10] =
235 {
236     // Lum class V (main sequence)
237     {
238         -0.31f, -0.37f, -0.42f, -0.50f, -0.55f,
239         -0.72f, -0.89f, -1.01f, -1.13f, -1.26f,
240     },
241     // Lum class III
242     {
243         -0.50f, -0.55f, -0.61f, -0.76f, -0.94f,
244         -1.02f, -1.09f, -1.17f, -1.20f, -1.22f,
245     },
246     // Lum class I
247     {
248         -0.50f, -0.56f, -0.61f, -0.75f, -0.90f,
249         -1.01f, -1.10f, -1.20f, -1.23f, -1.26f,
250     }
251 };
252 
253 static float bmag_correctionM[3][10] =
254 {
255     // Lum class V (main sequence)
256     {
257         -1.38f, -1.62f, -1.89f, -2.15f, -2.38f,
258         -2.73f, -3.21f, -3.46f, -4.10f, -4.40f,
259     },
260     // Lum class III
261     {
262         -1.25f, -1.44f, -1.62f, -1.87f, -2.22f,
263         -2.48f, -2.73f, -2.73f, -2.73f, -2.73f,
264     },
265     // Lum class I
266     {
267         -1.29f, -1.38f, -1.62f, -2.13f, -2.75f,
268         -3.47f, -3.90f, -3.90f, -3.90f, -3.90f,
269     }
270 };
271 
272 // Brown dwarf data from Grant Hutchison
273 static float bmag_correctionL[10] =
274 {
275     -4.6f, -4.9f, -5.0f, -5.2f, -5.4f, -5.9f, -6.1f, -6.7f, -7.4f, -8.2f,
276 };
277 
278 static float bmag_correctionT[10] =
279 {
280     -8.9f, -9.6f, -10.8f, -11.9f, -13.1f, -14.4f, -16.1f, -17.9f, -19.6f, -19.6f,
281 };
282 
283 // White dwarf data from Grant Hutchison; value for hypothetical
284 // 0 subclass is just duplicated from subclass 1.
285 static float bmag_correctionWD[10] =
286 {
287     -4.15f, -4.15f, -2.22f, -1.24f, -0.67f,
288     -0.32f, -0.13f, -0.04f, -0.03f, -0.09f,
289 };
290 
291 
292 // Stellar rotation by spectral and luminosity class.
293 // Tables from Grant Hutchison:
294 // "Most data are from Lang's _Astrophysical Data: Planets and Stars_ (I
295 // calculated from theoretical radii and observed rotation velocities), but
296 // with some additional information gleaned from elsewhere.
297 // A big scatter in rotation periods, of course, particularly in the K and
298 // early M dwarfs. I'm not hugely happy with the supergiant and giant rotation
299 // periods for K and M, either - they may be considerably slower yet, but it's
300 // obviously difficult to come by the data when the rotation velocity is too
301 // slow to obviously affect the spectra."
302 //
303 // I add missing values by interpolating linearly--certainly not the best
304 // technique, but adequate for our purposes.  The rotation rate of the Sun
305 // was used for spectral class G2.
306 
307 static float rotperiod_O[3][10] =
308 {
309     { 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f, 2.0f },
310     { 6.3f, 6.3f, 6.3f, 6.3f, 6.3f, 6.3f, 6.3f, 6.3f, 6.3f, 6.3f },
311     { 15.0f, 15.0f, 15.0f, 15.0f, 15.0f, 15.0f, 15.0f, 15.0f, 15.0f, 15.0f },
312 };
313 
314 static float rotperiod_B[3][10] =
315 {
316     { 2.0f, 1.8f, 1.6f, 1.4f, 1.1f, 0.8f, 0.8f, 0.8f, 0.8f, 0.7f },
317     { 6.3f, 5.6f, 5.0f, 4.3f, 3.7f, 3.1f, 2.9f, 2.8f, 2.7f, 2.6f },
318     { 15.0f, 24.0f, 33.0f, 42.0f, 52.0f, 63.0f, 65.0f, 67.0f, 70.0f, 72.0f },
319 };
320 
321 static float rotperiod_A[3][10] =
322 {
323     { 0.7f, 0.7f, 0.6f, 0.6f, 0.5f, 0.5f, 0.5f, 0.6f, 0.6f, 0.7f },
324     { 2.5f, 2.3f, 2.1f, 1.9f, 1.7f, 1.6f, 1.6f, 1.7f, 1.7f, 1.8f },
325     { 75.0f, 77.0f, 80.0f, 82.0f, 85.0f, 87.0f, 95.0f, 104.0f, 115.0f, 125.0f },
326 };
327 
328 static float rotperiod_F[3][10] =
329 {
330     { 0.7f, 0.7f, 0.6f, 0.6f, 0.5f, 0.5f, 0.5f, 0.6f, 0.6f, 0.7f },
331     { 1.9f, 2.5f, 3.0f, 3.5f, 4.0f, 4.6f, 5.6f, 6.7f, 7.8f, 8.9f },
332     { 135.0f, 141.0f, 148.0f, 155.0f, 162.0f, 169.0f, 175.0f, 182.0f, 188.0f, 195.0f },
333 };
334 
335 static float rotperiod_G[3][10] =
336 {
337     { 11.1f, 18.2f, 25.4f, 24.7f, 24.0f, 23.3f, 23.0f, 22.7f, 22.3f, 21.9f },
338     { 10.0f, 13.0f, 16.0f, 19.0f, 22.0f, 25.0f, 28.0f, 31.0f, 33.0f, 35.0f },
339     { 202.0f, 222.0f, 242.0f, 262.0f, 282.0f,
340       303.0f, 323.0f, 343.0f, 364.0f, 384.0f },
341 };
342 
343 static float rotperiod_K[3][10] =
344 {
345     { 21.5f, 20.8f, 20.2f, 19.4f, 18.8f, 18.2f, 17.6f, 17.0f, 16.4f, 15.8f },
346     { 38.0f, 43.0f, 48.0f, 53.0f, 58.0f, 63.0f, 71.0f, 78.0f, 86.0f, 93.0f },
347     { 405.0f, 526.0f, 648.0f, 769.0f, 891.0f,
348       1012.0f, 1063.0f, 1103.0f, 1154.0f, 1204.0f },
349 };
350 
351 static float rotperiod_M[3][10] =
352 {
353     { 15.2f, 12.4f, 9.6f, 6.8f, 4.0f, 1.3f, 1.0f, 0.7f, 0.4f, 0.2f },
354     { 101.0f, 101.0f, 101.0f, 101.0f, 101.0f, 101.0f, 101.0f, 101.0f, 101.0f, 101.0f },
355     { 1265.0f, 1265.0f, 1265.0f, 1265.0f, 1265.0f,
356       1265.0f, 1265.0f, 1265.0f, 1265.0f, 1265.0f },
357 };
358 
359 
360 const char* LumClassNames[StellarClass::Lum_Count] = {
361     "I-a0", "I-a", "I-b", "II", "III", "IV", "V", "VI", ""
362 };
363 
364 const char* SubclassNames[11] = {
365     "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", ""
366 };
367 
368 const char* SpectralClassNames[StellarClass::NormalClassCount] = {
369     "O", "B", "A", "F", "G", "K", "M", "R",
370     "S", "N", "WC", "WN", "?", "L", "T", "C",
371 };
372 
373 const char* WDSpectralClassNames[StellarClass::WDClassCount] = {
374     "DA", "DB", "DC", "DO", "DQ", "DZ", "D", "DX",
375 };
376 
377 
378 StarDetails*
GetStarDetails(const StellarClass & sc)379 StarDetails::GetStarDetails(const StellarClass& sc)
380 {
381     switch (sc.getStarType())
382     {
383     case StellarClass::NormalStar:
384         return GetNormalStarDetails(sc.getSpectralClass(),
385                                     sc.getSubclass(),
386                                     sc.getLuminosityClass());
387 
388     case StellarClass::WhiteDwarf:
389         return GetWhiteDwarfDetails(sc.getSpectralClass(),
390                                     sc.getSubclass());
391     case StellarClass::NeutronStar:
392         return GetNeutronStarDetails();
393     case StellarClass::BlackHole:
394         return GetBlackHoleDetails();
395     default:
396         return NULL;
397     }
398 }
399 
400 
401 StarDetails*
CreateStandardStarType(const std::string & specTypeName,float _temperature,float _rotationPeriod)402 StarDetails::CreateStandardStarType(const std::string& specTypeName,
403                                     float _temperature,
404                                     float _rotationPeriod)
405 
406 {
407     StarDetails* details = new StarDetails();
408 
409     details->setTemperature(_temperature);
410     details->setSpectralType(specTypeName);
411 
412     details->setRotationModel(new UniformRotationModel(_rotationPeriod,
413                                                        0.0f,
414                                                        astro::J2000,
415                                                        0.0f,
416                                                        0.0f));
417 
418     return details;
419 }
420 
421 
422 StarDetails*
GetNormalStarDetails(StellarClass::SpectralClass specClass,unsigned int subclass,StellarClass::LuminosityClass lumClass)423 StarDetails::GetNormalStarDetails(StellarClass::SpectralClass specClass,
424                                   unsigned int subclass,
425                                   StellarClass::LuminosityClass lumClass)
426 {
427     if (normalStarDetails == NULL)
428     {
429         unsigned int nTypes = StellarClass::Spectral_Count * 11 *
430             StellarClass::Lum_Count;
431         normalStarDetails = new StarDetails*[nTypes];
432         for (unsigned int i = 0; i < nTypes; i++)
433             normalStarDetails[i] = NULL;
434     }
435 
436     if (subclass > StellarClass::Subclass_Unknown)
437         subclass = StellarClass::Subclass_Unknown;
438 
439     uint index = subclass + (specClass + lumClass * StellarClass::Spectral_Count) * 11;
440     if (normalStarDetails[index] == NULL)
441     {
442         char name[16];
443         if ((lumClass == StellarClass::Lum_VI) &&
444             (specClass >= StellarClass::Spectral_O) && (specClass <= StellarClass::Spectral_A))
445         {
446             // Hot subdwarfs are prefixed with "sd", while cool subdwarfs use
447             // luminosity class VI, per recommendations in arXiv:0805.2567v1
448             sprintf(name, "sd%s%s",
449                     SpectralClassNames[specClass],
450                     SubclassNames[subclass]);
451         }
452         else
453         {
454             sprintf(name, "%s%s%s",
455                     SpectralClassNames[specClass],
456                     SubclassNames[subclass],
457                     LumClassNames[lumClass]);
458         }
459 
460         // Use the same properties for an unknown subclass as for subclass 5
461         if (subclass == StellarClass::Subclass_Unknown)
462         {
463             // Since early O and Wolf-Rayet stars are exceedingly rare,
464             // use temperature of the more common late types when the subclass
465             // is unspecified in the spectral type.  For other stars, default
466             // to subclass 5.
467             switch (specClass)
468             {
469             case StellarClass::Spectral_O:
470             case StellarClass::Spectral_WN:
471             case StellarClass::Spectral_WC:
472                 subclass = 9;
473                 break;
474             default:
475                 subclass = 5;
476                 break;
477             }
478         }
479 
480         unsigned int lumIndex = 0;
481         switch (lumClass)
482         {
483         case StellarClass::Lum_Ia0:
484         case StellarClass::Lum_Ia:
485         case StellarClass::Lum_Ib:
486         case StellarClass::Lum_II:
487             lumIndex = 2;
488             break;
489         case StellarClass::Lum_III:
490         case StellarClass::Lum_IV:
491             lumIndex = 1;
492             break;
493         case StellarClass::Lum_V:
494         case StellarClass::Lum_VI:
495         case StellarClass::Lum_Unknown:
496             lumIndex = 0;
497             break;
498 
499         default: break;  // Do nothing, but prevent GCC4 warnings (Beware: potentially dangerous)
500         }
501 
502         float temp = 0.0f;
503         switch (specClass)
504         {
505         case StellarClass::Spectral_O:
506             temp = tempO[lumIndex][subclass];
507             break;
508         case StellarClass::Spectral_B:
509             temp = tempB[lumIndex][subclass];
510             break;
511         case StellarClass::Spectral_Unknown:
512         case StellarClass::Spectral_A:
513             temp = tempA[lumIndex][subclass];
514             break;
515         case StellarClass::Spectral_F:
516             temp = tempF[lumIndex][subclass];
517             break;
518         case StellarClass::Spectral_G:
519             temp = tempG[lumIndex][subclass];
520             break;
521         case StellarClass::Spectral_K:
522             temp = tempK[lumIndex][subclass];
523             break;
524         case StellarClass::Spectral_M:
525             temp = tempM[lumIndex][subclass];
526             break;
527         case StellarClass::Spectral_R:
528             temp = tempK[lumIndex][subclass];
529             break;
530         case StellarClass::Spectral_S:
531             temp = tempM[lumIndex][subclass];
532             break;
533         case StellarClass::Spectral_N:
534             temp = tempM[lumIndex][subclass];
535             break;
536         case StellarClass::Spectral_C:
537             temp = tempM[lumIndex][subclass];
538             break;
539         case StellarClass::Spectral_WN:
540             temp = tempWN[subclass];
541             break;
542         case StellarClass::Spectral_WC:
543             temp = tempWC[subclass];
544             break;
545         case StellarClass::Spectral_L:
546             temp = tempL[subclass];
547             break;
548         case StellarClass::Spectral_T:
549             temp = tempT[subclass];
550             break;
551 
552         default: break;  // Do nothing, but prevent GCC4 warnings (Beware: potentially dangerous)
553         }
554 
555         float bmagCorrection = 0.0f;
556         float period = 1.0f;
557         switch (specClass)
558         {
559         case StellarClass::Spectral_O:
560             period = rotperiod_O[lumIndex][subclass];
561             bmagCorrection = bmag_correctionO[lumIndex][subclass];
562             break;
563         case StellarClass::Spectral_B:
564             period = rotperiod_B[lumIndex][subclass];
565             bmagCorrection = bmag_correctionB[lumIndex][subclass];
566             break;
567         case StellarClass::Spectral_Unknown:
568         case StellarClass::Spectral_A:
569             period = rotperiod_A[lumIndex][subclass];
570             bmagCorrection = bmag_correctionA[lumIndex][subclass];
571             break;
572         case StellarClass::Spectral_F:
573             period = rotperiod_F[lumIndex][subclass];
574             bmagCorrection = bmag_correctionF[lumIndex][subclass];
575             break;
576         case StellarClass::Spectral_G:
577             period = rotperiod_G[lumIndex][subclass];
578             bmagCorrection = bmag_correctionG[lumIndex][subclass];
579             break;
580         case StellarClass::Spectral_K:
581             period = rotperiod_K[lumIndex][subclass];
582             bmagCorrection = bmag_correctionK[lumIndex][subclass];
583             break;
584         case StellarClass::Spectral_M:
585             period = rotperiod_M[lumIndex][subclass];
586             bmagCorrection = bmag_correctionM[lumIndex][subclass];
587             break;
588 
589         case StellarClass::Spectral_R:
590         case StellarClass::Spectral_S:
591         case StellarClass::Spectral_N:
592         case StellarClass::Spectral_C:
593             period = rotperiod_M[lumIndex][subclass];
594             bmagCorrection = bmag_correctionM[lumIndex][subclass];
595             break;
596 
597         case StellarClass::Spectral_WC:
598         case StellarClass::Spectral_WN:
599             period = rotperiod_O[lumIndex][subclass];
600             bmagCorrection = bmag_correctionO[lumIndex][subclass];
601             break;
602 
603         case StellarClass::Spectral_L:
604             // Assume that brown dwarfs are fast rotators like late M dwarfs
605             period = 0.2f;
606             bmagCorrection = bmag_correctionL[subclass];
607             break;
608 
609         case StellarClass::Spectral_T:
610             // Assume that brown dwarfs are fast rotators like late M dwarfs
611             period = 0.2f;
612             bmagCorrection = bmag_correctionT[subclass];
613             break;
614 
615         default: break;  // Do nothing, but prevent GCC4 warnings (Beware: potentially dangerous)
616         }
617 
618         normalStarDetails[index] = CreateStandardStarType(name, temp, period);
619         normalStarDetails[index]->setBolometricCorrection(bmagCorrection);
620 
621         MultiResTexture starTex = starTextures.starTex[specClass];
622         if (!starTex.isValid())
623             starTex = starTextures.defaultTex;
624         normalStarDetails[index]->setTexture(starTex);
625     }
626 
627     return normalStarDetails[index];
628 }
629 
630 
631 StarDetails*
GetWhiteDwarfDetails(StellarClass::SpectralClass specClass,unsigned int subclass)632 StarDetails::GetWhiteDwarfDetails(StellarClass::SpectralClass specClass,
633                                   unsigned int subclass)
634 {
635     // Hack assumes all WD types are consecutive
636     unsigned int scIndex = static_cast<unsigned int>(specClass) -
637         StellarClass::FirstWDClass;
638 
639     if (whiteDwarfDetails == NULL)
640     {
641         unsigned int nTypes =
642             StellarClass::WDClassCount * StellarClass::SubclassCount;
643         whiteDwarfDetails = new StarDetails*[nTypes];
644         for (unsigned int i = 0; i < nTypes; i++)
645             whiteDwarfDetails[i] = NULL;
646     }
647 
648     if (subclass > StellarClass::Subclass_Unknown)
649         subclass = StellarClass::Subclass_Unknown;
650 
651     uint index = subclass + (scIndex * StellarClass::SubclassCount);
652     if (whiteDwarfDetails[index] == NULL)
653     {
654         char name[16];
655         sprintf(name, "%s%s",
656                 WDSpectralClassNames[scIndex],
657                 SubclassNames[subclass]);
658 
659         float temp;
660         float bmagCorrection;
661         // subclass is always >= 0:
662         if (subclass <= 9)
663         {
664             temp = tempWD[subclass];
665             bmagCorrection = bmag_correctionWD[subclass];
666         }
667         else
668         {
669             // Treat unknown as subclass 5
670             temp = tempWD[5];
671             bmagCorrection = bmag_correctionWD[5];
672         }
673 
674         // Assign white dwarfs a rotation period of half an hour; very
675         // rough, as white rotation rates vary a lot.
676         float period = 1.0f / 48.0f;
677 
678         whiteDwarfDetails[index] = CreateStandardStarType(name, temp, period);
679         MultiResTexture starTex = starTextures.starTex[StellarClass::Spectral_D];
680         if (!starTex.isValid())
681             starTex = starTextures.defaultTex;
682         whiteDwarfDetails[index]->setTexture(starTex);
683         whiteDwarfDetails[index]->setBolometricCorrection(bmagCorrection);
684     }
685 
686     return whiteDwarfDetails[index];
687 }
688 
689 
690 StarDetails*
GetNeutronStarDetails()691 StarDetails::GetNeutronStarDetails()
692 {
693     if (neutronStarDetails == NULL)
694     {
695         // The default neutron star has a rotation period of one second,
696         // surface temperature of five million K.
697         neutronStarDetails = CreateStandardStarType("Q", 5000000.0f,
698                                                     1.0f / 86400.0f);
699         neutronStarDetails->setRadius(10.0f);
700         neutronStarDetails->addKnowledge(KnowRadius);
701         MultiResTexture starTex = starTextures.neutronStarTex;
702         if (!starTex.isValid())
703             starTex = starTextures.defaultTex;
704         neutronStarDetails->setTexture(starTex);
705     }
706 
707     return neutronStarDetails;
708 }
709 
710 
711 StarDetails*
GetBlackHoleDetails()712 StarDetails::GetBlackHoleDetails()
713 {
714     if (blackHoleDetails == NULL)
715     {
716         // Default black hole parameters are based on a one solar mass
717         // black hole.
718         // The temperature is computed from the equation:
719         //      T=h_bar c^3/(8 pi G k m)
720         blackHoleDetails = CreateStandardStarType("X", 6.15e-8f,
721                                                   1.0f / 86400.0f);
722         blackHoleDetails->setRadius(2.9f);
723         blackHoleDetails->addKnowledge(KnowRadius);
724     }
725 
726     return blackHoleDetails;
727 }
728 
729 
730 StarDetails*
GetBarycenterDetails()731 StarDetails::GetBarycenterDetails()
732 {
733 
734     if (barycenterDetails == NULL)
735     {
736         barycenterDetails = CreateStandardStarType("Bary", 1.0f, 1.0f);
737         barycenterDetails->setRadius(0.001f);
738         barycenterDetails->addKnowledge(KnowRadius);
739         barycenterDetails->setVisibility(false);
740     }
741 
742     return barycenterDetails;
743 }
744 
745 
746 void
SetStarTextures(const StarTextureSet & _starTextures)747 StarDetails::SetStarTextures(const StarTextureSet& _starTextures)
748 {
749     starTextures = _starTextures;
750 }
751 
752 
StarDetails()753 StarDetails::StarDetails() :
754     radius(0.0f),
755     temperature(0.0f),
756     bolometricCorrection(0.0f),
757     knowledge(0u),
758     visible(true),
759     texture(texture),
760     geometry(InvalidResource),
761     orbit(NULL),
762     orbitalRadius(0.0f),
763     barycenter(NULL),
764     rotationModel(NULL),
765     semiAxes(1.0f, 1.0f, 1.0f),
766     infoURL(NULL),
767     orbitingStars(NULL),
768     isShared(true)
769 {
770     spectralType[0] = '\0';
771 }
772 
773 
StarDetails(const StarDetails & sd)774 StarDetails::StarDetails(const StarDetails& sd) :
775     radius(sd.radius),
776     temperature(sd.temperature),
777     bolometricCorrection(sd.bolometricCorrection),
778     knowledge(sd.knowledge),
779     visible(sd.visible),
780     texture(sd.texture),
781     geometry(sd.geometry),
782     orbit(sd.orbit),
783     orbitalRadius(sd.orbitalRadius),
784     barycenter(sd.barycenter),
785     rotationModel(sd.rotationModel),
786     semiAxes(sd.semiAxes),
787     infoURL(NULL),
788     orbitingStars(NULL),
789     isShared(false)
790 {
791     assert(sd.isShared);
792     memcpy(spectralType, sd.spectralType, sizeof(spectralType));
793     if (sd.infoURL != NULL)
794         infoURL = new string(*sd.infoURL);
795 }
796 
797 
~StarDetails()798 StarDetails::~StarDetails()
799 {
800     delete orbitingStars;
801     delete infoURL;
802 }
803 
804 
805 /*! Return the InfoURL. If the InfoURL has not been set, this method
806  *  returns an empty string.
807  */
808 const std::string&
getInfoURL() const809 StarDetails::getInfoURL() const
810 {
811     if (infoURL != NULL)
812         return *infoURL;
813     else
814         return DEFAULT_INFO_URL;
815 }
816 
817 
818 void
setRadius(float _radius)819 StarDetails::setRadius(float _radius)
820 {
821     radius = _radius;
822 }
823 
824 
825 void
setTemperature(float _temperature)826 StarDetails::setTemperature(float _temperature)
827 {
828     temperature = _temperature;
829 }
830 
831 
832 void
setSpectralType(const std::string & s)833 StarDetails::setSpectralType(const std::string& s)
834 {
835     strncpy(spectralType, s.c_str(), sizeof(spectralType));
836     spectralType[sizeof(spectralType) - 1] = '\0';
837 }
838 
839 
840 void
setKnowledge(uint32 _knowledge)841 StarDetails::setKnowledge(uint32 _knowledge)
842 {
843     knowledge = _knowledge;
844 }
845 
846 
847 void
addKnowledge(uint32 _knowledge)848 StarDetails::addKnowledge(uint32 _knowledge)
849 {
850     knowledge |= _knowledge;
851 }
852 
853 
854 void
setBolometricCorrection(float correction)855 StarDetails::setBolometricCorrection(float correction)
856 {
857     bolometricCorrection = correction;
858 }
859 
860 
861 void
setTexture(const MultiResTexture & tex)862 StarDetails::setTexture(const MultiResTexture& tex)
863 {
864     texture = tex;
865 }
866 
867 
868 void
setGeometry(ResourceHandle rh)869 StarDetails::setGeometry(ResourceHandle rh)
870 {
871     geometry = rh;
872 }
873 
874 
875 void
setOrbit(Orbit * o)876 StarDetails::setOrbit(Orbit* o)
877 {
878     orbit = o;
879     computeOrbitalRadius();
880 }
881 
882 
883 void
setOrbitBarycenter(Star * bc)884 StarDetails::setOrbitBarycenter(Star* bc)
885 {
886     barycenter = bc;
887     computeOrbitalRadius();
888 }
889 
890 
891 void
setOrbitalRadius(float r)892 StarDetails::setOrbitalRadius(float r)
893 {
894     if (orbit != NULL)
895         orbitalRadius = r;
896 }
897 
898 
899 void
computeOrbitalRadius()900 StarDetails::computeOrbitalRadius()
901 {
902     if (orbit == NULL)
903     {
904         orbitalRadius = 0.0f;
905     }
906     else
907     {
908         orbitalRadius = (float) astro::kilometersToLightYears(orbit->getBoundingRadius());
909         if (barycenter != NULL)
910             orbitalRadius += barycenter->getOrbitalRadius();
911     }
912 }
913 
914 
915 void
setVisibility(bool b)916 StarDetails::setVisibility(bool b)
917 {
918     visible = b;
919 }
920 
921 
922 void
setRotationModel(const RotationModel * rm)923 StarDetails::setRotationModel(const RotationModel* rm)
924 {
925     rotationModel = rm;
926 }
927 
928 
929 /*! Set the InfoURL for this star.
930 */
931 void
setInfoURL(const string & _infoURL)932 StarDetails::setInfoURL(const string& _infoURL)
933 {
934     if (_infoURL.empty())
935     {
936         // Save space in the common case--no InfoURL--by not
937         // allocating a string.
938         delete infoURL;
939         infoURL = NULL;
940     }
941     else
942     {
943         // Allocate the new string before freeing the old one, so we don't crash
944         // in the event the caller does something like:
945         // star->setInfoURL(star->getInfoURL());
946         string* oldURL = infoURL;
947         infoURL = new string(_infoURL);
948         delete oldURL;
949     }
950 }
951 
952 
~Star()953 Star::~Star()
954 {
955     // TODO: Implement reference counting for StarDetails objects so that
956     // we can enable this.
957 #if 0
958     if (!details->shared())
959         delete details;
960 #endif
961 }
962 
963 
964 // Return the radius of the star in kilometers
getRadius() const965 float Star::getRadius() const
966 {
967     if (details->getKnowledge(StarDetails::KnowRadius))
968         return details->getRadius();
969 
970 #ifdef NO_BOLOMETRIC_MAGNITUDE_CORRECTION
971     // Use the Stefan-Boltzmann law to estimate the radius of a
972     // star from surface temperature and luminosity
973     return SOLAR_RADIUS * (float) sqrt(getLuminosity()) *
974         square(SOLAR_TEMPERATURE / getTemperature());
975 #else
976     // Calculate the luminosity of the star from the bolometric, not the
977     // visual magnitude of the star.
978     float solarBMag = SOLAR_ABSMAG + bmag_correctionG[0][2];
979     float bmag = getBolometricMagnitude();
980     float boloLum = (float) exp((solarBMag - bmag) / LN_MAG);
981 
982     // Use the Stefan-Boltzmann law to estimate the radius of a
983     // star from surface temperature and luminosity
984     return SOLAR_RADIUS * (float) sqrt(boloLum) *
985         square(SOLAR_TEMPERATURE / getTemperature());
986 #endif
987 }
988 
989 
990 void
setEllipsoidSemiAxes(const Vec3f & v)991 StarDetails::setEllipsoidSemiAxes(const Vec3f& v)
992 {
993     semiAxes = v;
994 }
995 
996 
997 bool
shared() const998 StarDetails::shared() const
999 {
1000     return isShared;
1001 }
1002 
1003 
1004 void
addOrbitingStar(Star * star)1005 StarDetails::addOrbitingStar(Star* star)
1006 {
1007     assert(!shared());
1008     if (orbitingStars == NULL)
1009         orbitingStars = new vector<Star*>();
1010     orbitingStars->push_back(star);
1011 }
1012 
1013 
1014 /*! Get the position of the star in the universal coordinate system.
1015  */
1016 UniversalCoord
getPosition(double t) const1017 Star::getPosition(double t) const
1018 {
1019     const Orbit* orbit = getOrbit();
1020     if (!orbit)
1021     {
1022         return UniversalCoord(position.x * 1.0e6,
1023                               position.y * 1.0e6,
1024                               position.z * 1.0e6);
1025     }
1026     else
1027     {
1028         const Star* barycenter = getOrbitBarycenter();
1029 
1030         if (barycenter == NULL)
1031         {
1032             Point3d barycenterPos(position.x * 1.0e6,
1033                                   position.y * 1.0e6,
1034                                   position.z * 1.0e6);
1035 
1036             return UniversalCoord(barycenterPos) +
1037                 ((orbit->positionAtTime(t) - Point3d(0.0, 0.0, 0.0)) *
1038                  astro::kilometersToMicroLightYears(1.0));
1039         }
1040         else
1041         {
1042             return barycenter->getPosition(t) +
1043                 ((orbit->positionAtTime(t) - Point3d(0.0, 0.0, 0.0)) *
1044                  astro::kilometersToMicroLightYears(1.0));
1045         }
1046     }
1047 }
1048 
1049 
1050 UniversalCoord
getOrbitBarycenterPosition(double t) const1051 Star::getOrbitBarycenterPosition(double t) const
1052 {
1053     const Star* barycenter = getOrbitBarycenter();
1054 
1055     if (barycenter == NULL)
1056     {
1057         Point3d barycenterPos(position.x * 1.0e6,
1058                               position.y * 1.0e6,
1059                               position.z * 1.0e6);
1060 
1061         return UniversalCoord(barycenterPos);
1062     }
1063     else
1064     {
1065         return barycenter->getPosition(t);
1066     }
1067 }
1068 
1069 
1070 /*! Get the velocity of the star in the universal coordinate system.
1071  */
1072 Vec3d
getVelocity(double t) const1073 Star::getVelocity(double t) const
1074 {
1075     const Orbit* orbit = getOrbit();
1076     if (!orbit)
1077     {
1078 		// The star doesn't have a defined orbit, so the velocity is just
1079 		// zero. (This will change when stellar proper motion is implemented.)
1080 		return Vec3d(0.0, 0.0, 0.0);
1081     }
1082     else
1083     {
1084         const Star* barycenter = getOrbitBarycenter();
1085 
1086         if (barycenter == NULL)
1087         {
1088 			// Star orbit is defined around a fixed point, so the total velocity
1089 			// is just the star's orbit velocity.
1090 			return orbit->velocityAtTime(t);
1091         }
1092         else
1093         {
1094 			// Sum the star's orbital velocity and the velocity of the barycenter.
1095             return barycenter->getVelocity(t) + orbit->velocityAtTime(t);
1096         }
1097     }
1098 }
1099 
1100 
1101 MultiResTexture
getTexture() const1102 Star::getTexture() const
1103 {
1104     return details->getTexture();
1105 }
1106 
1107 
1108 ResourceHandle
getGeometry() const1109 Star::getGeometry() const
1110 {
1111     return details->getGeometry();
1112 }
1113 
1114 
1115 /*! Return the InfoURL. If the InfoURL has not been set, this method
1116 *  returns an empty string.
1117 */
1118 const string&
getInfoURL() const1119 Star::getInfoURL() const
1120 {
1121     return details->getInfoURL();
1122 }
1123 
1124 
setCatalogNumber(uint32 n)1125 void Star::setCatalogNumber(uint32 n)
1126 {
1127     catalogNumber = n;
1128 }
1129 
setPosition(float x,float y,float z)1130 void Star::setPosition(float x, float y, float z)
1131 {
1132     position = Point3f(x, y, z);
1133 }
1134 
setPosition(Point3f p)1135 void Star::setPosition(Point3f p)
1136 {
1137     position = p;
1138 }
1139 
setAbsoluteMagnitude(float mag)1140 void Star::setAbsoluteMagnitude(float mag)
1141 {
1142     absMag = mag;
1143 }
1144 
1145 
getApparentMagnitude(float ly) const1146 float Star::getApparentMagnitude(float ly) const
1147 {
1148     return astro::absToAppMag(absMag, ly);
1149 }
1150 
1151 
getLuminosity() const1152 float Star::getLuminosity() const
1153 {
1154     return astro::absMagToLum(absMag);
1155 }
1156 
setLuminosity(float lum)1157 void Star::setLuminosity(float lum)
1158 {
1159     absMag = astro::lumToAbsMag(lum);
1160 }
1161 
getDetails() const1162 StarDetails* Star::getDetails() const
1163 {
1164     return details;
1165 }
1166 
setDetails(StarDetails * sd)1167 void Star::setDetails(StarDetails* sd)
1168 {
1169     // TODO: delete existing details if they aren't shared
1170     details = sd;
1171 }
1172 
setOrbitBarycenter(Star * s)1173 void Star::setOrbitBarycenter(Star* s)
1174 {
1175     if (details->shared())
1176         details = new StarDetails(*details);
1177     details->setOrbitBarycenter(s);
1178 }
1179 
computeOrbitalRadius()1180 void Star::computeOrbitalRadius()
1181 {
1182     details->computeOrbitalRadius();
1183 }
1184 
1185 void
setRotationModel(const RotationModel * rm)1186 Star::setRotationModel(const RotationModel* rm)
1187 {
1188     details->setRotationModel(rm);
1189 }
1190 
1191 void
addOrbitingStar(Star * star)1192 Star::addOrbitingStar(Star* star)
1193 {
1194     if (details->shared())
1195         details = new StarDetails(*details);
1196     details->addOrbitingStar(star);
1197 }
1198