开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > OpenGL > 查看源码
buildstardb.cpp
资源名称：celestia-1.6.0.tar.gz　[点击查看]
上传用户：center1979
上传日期：2022-07-26
资源大小：50633k
文件大小：38k
源码类别：

OpenGL

开发平台：

Visual C++

buildstardb.cpp：源码内容
  1. // buildstardb.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.
  10. #include <string>
  11. #include <vector>
  12. #include <map>
  13. #include <algorithm>
  14. #include <iostream>
  15. #include <fstream>
  16. #include <cstdio>
  17. #include <assert.h>
  18. #include <unistd.h>
  19. #include "celengine/stardb.h"
  21. using namespace std;
  24. #ifdef _WIN32
  25. static string MainDatabaseFile("hip_main.dat");
  26. static string TychoDatabaseFile("tyc_main.dat");
  27. static string ComponentDatabaseFile("h_dm_com.dat");
  28. static string OrbitalDatabase("hip_dm_o.dat");
  29. #else
  30. #ifndef MACOSX_PB
  31. #include <config.h>
  32. #endif /*MACOSX_PB*/
  33. static string MainDatabaseFile(HIP_DATA_DIR "/hip_main.dat");
  34. static string TychoDatabaseFile(HIP_DATA_DIR "/tyc_main.dat");
  35. static string ComponentDatabaseFile(HIP_DATA_DIR "/h_dm_com.dat");
  36. static string OrbitalDatabase(HIP_DATA_DIR "/hip_dm_o.dat");
  37. #endif
  39. static const int HipStarRecordLength = 451;
  40. static const int HipComponentRecordLength = 239;
  41. static const int TycStarRecordLength = 351;
  43. static uint32 NullCCDMIdentifier = 0xffffffff;
  44. static uint32 NullCatalogNumber = 0xffffffff;
  45. static int verbose, dropstars;
  48. class HipparcosStar
  49. {
  50. public:
  51.     HipparcosStar();
  53.     void write(ostream&);
  54.     void analyze();
  56.     uint32 HIPCatalogNumber;
  57.     uint32 HDCatalogNumber;
  58.     float ascension;
  59.     float declination;
  60.     float parallax;
  61.     float appMag;
  62.     StellarClass stellarClass;
  64.     uint32 CCDMIdentifier;
  65.     uint8 starsWithCCDM;
  66.     uint8 nComponents;
  67.     uint8 parallaxError;
  68.     uint tycline;
  69.     uint status;
  70.     float e_RA;  // Errors in Right Ascension,
  71.     float e_DE;  // Declination,
  72.     float e_Mag; // and Magnitude
  73. };
  75. HipparcosStar::HipparcosStar() :
  76.     HIPCatalogNumber(NullCatalogNumber),
  77.     HDCatalogNumber(NullCatalogNumber),
  78.     ascension(0.0f),
  79.     declination(0.0f),
  80.     parallax(0.0f),
  81.     appMag(0.0f),
  82.     CCDMIdentifier(NullCCDMIdentifier),
  83.     starsWithCCDM(0),
  84.     nComponents(1),
  85.     parallaxError(0),
  86.     tycline(0),
  87.     e_RA(0.0f),
  88.     e_DE(0.0f),
  89.     e_Mag(0.0f)
  90. {
  91. }
  93. template<class T> void binwrite(ostream& out, T x)
  94. {
  95.     out.write(reinterpret_cast<char*>(&x), sizeof(T));
  96. }
  98. void HipparcosStar::write(ostream& out)
  99. {
  100.     if (status>=2)
  101.         return;
  102. #if 0
  103.     if (HDCatalogNumber != NullCatalogNumber)
  104.         binwrite(out, HDCatalogNumber);
  105.     else
  106.         binwrite(out, HIPCatalogNumber | 0x10000000);
  107. #endif
  108.     binwrite(out, HIPCatalogNumber);
  109.     binwrite(out, HDCatalogNumber);
  110.     binwrite(out, ascension);
  111.     binwrite(out, declination);
  112.     binwrite(out, parallax);
  113.     binwrite(out, (short) (appMag * 256));
  114.     binwrite(out, stellarClass);
  115.     binwrite(out, parallaxError);
  116. }
  118.                     // Statistic Values:
  119. float s_er;         // Sum of Error in RA
  120. float s_erq;        // Sum of Squares of Errors in RA
  121. unsigned int n_er;  // Number of Error Values
  122. float s_ed, s_edq;  // Ditto for Declination
  123. unsigned int n_ed;
  124. unsigned int n_drop,n_dub; // number of stars to drop, number dubious
  126. void HipparcosStar::analyze()
  127. {
  128.     int dubious=0;
  129.     status=0;
  130.     if ((parallax) <= 0.4 && ((dropstars==0) || (dropstars==1 && appMag<6.0)))
  131.         parallax=0.4;  /* fix strange paralaxes so the stars aren't *way*
  132.                           out there. A parallax of 0.4 will put them at
  133.                           just a touch above 8154 LY away. */
  134.     if (parallax <= 0.0)
  135.         dubious+=400;
  136.     else if (parallax<0.2)
  137.         dubious+=4;
  138.     else if (parallax<0.4)
  139.         dubious+=2;
  140.     if (parallax<=parallaxError)
  141.         dubious+=2;
  142.     else if (parallax<=(2*parallaxError))
  143.         dubious++;
  144.     if (e_RA<1000.0)
  145.     {
  146.         s_er+=e_RA;
  147.         s_erq+=square(e_RA);
  148.         n_er++;
  149.         if (e_RA>20.0)
  150.             dubious+=100;
  151.         else if (e_RA>15.0)
  152.             dubious+=2;
  153.         else if (e_RA>10.0)
  154.             dubious++;
  155.     }
  156.     else
  157.         dubious+=4;  /* No error given, assume it's rather dubious */
  159.     if (e_DE<1000.0)
  160.     {
  161.         s_ed+=e_DE;
  162.         s_edq+=square(e_DE);
  163.         n_ed++;
  164.         if (e_DE>20.0)
  165.             dubious+=100;
  166.         else if (e_DE>15.0)
  167.             dubious+=2;
  168.         else if (e_DE>10.0)
  169.             dubious++;
  170.     }
  171.     else
  172.         dubious+=4;  /* No error given, assume it's rather dubious */
  174.     if (dubious>3)
  175.     {
  176.         status=1;
  177.         n_dub++;
  178.     }
  179.     if ((dubious>5) && (dropstars) && (!(dropstars==1 && appMag<6.0)))
  180.     {
  181.         status=2;
  182.         n_drop++;
  183.     }
  184. }
  187. bool operator<(const HipparcosStar& a, const HipparcosStar& b)
  188. {
  189.     return a.HIPCatalogNumber < b.HIPCatalogNumber;
  190. }
  192. struct HIPCatalogComparePredicate
  193. {
  194.     HIPCatalogComparePredicate() : dummy(0)
  195.     {
  196.     }
  198.     bool operator()(const HipparcosStar* star0, const HipparcosStar* star1) const
  199.     {
  200.         return star0->HIPCatalogNumber < star1->HIPCatalogNumber;
  201.     }
  203.     bool operator()(const HipparcosStar* star0, uint32 hip)
  204.     {
  205.         return star0->HIPCatalogNumber < hip;
  206.     }
  208.     int dummy;
  209. };
  212. class MultistarSystem
  213. {
  214. public:
  215.     int nStars; // Never greater than four in the HIPPARCOS catalog
  216.     HipparcosStar* stars[4];
  217. };
  220. class HipparcosComponent
  221. {
  222. public:
  223.     HipparcosComponent();
  225.     HipparcosStar* star;
  226.     char componentID;
  227.     char refComponentID;
  228.     float ascension;
  229.     float declination;
  230.     float appMag;
  231.     float bMag;
  232.     float vMag;
  233.     bool hasBV;
  234.     float positionAngle;
  235.     float separation;
  236. };
  238. HipparcosComponent::HipparcosComponent() :
  239.     star(NULL),
  240.     componentID('A'),
  241.     refComponentID('A'),
  242.     appMag(0.0f),
  243.     bMag(0.0f),
  244.     vMag(0.0f),
  245.     hasBV(false),
  246.     positionAngle(0.0f),
  247.     separation(0.0f)
  248. {
  249. }
  252. vector<HipparcosStar> stars;
  253. vector<HipparcosStar> companions;
  254. vector<HipparcosComponent> components;
  255. vector<HipparcosStar*> starIndex;
  257. typedef map<uint32, MultistarSystem*> MultistarSystemCatalog;
  258. MultistarSystemCatalog starSystems;
  261. HipparcosStar* findStar(uint32 hip)
  262. {
  263.     HIPCatalogComparePredicate pred;
  265.     vector<HipparcosStar*>::iterator iter = lower_bound(starIndex.begin(),
  266.                                                         starIndex.end(),
  267.                                                         hip, pred);
  268.     if (iter == starIndex.end())
  269.         return NULL;
  270.     HipparcosStar* star = *iter;
  271.     if (star->HIPCatalogNumber == hip)
  272.         return star;
  273.     else
  274.         return NULL;
  275. }
  278. StellarClass ParseStellarClass(char *starType)
  279. {
  280.     StellarClass::StarType type = StellarClass::NormalStar;
  281.     StellarClass::SpectralClass specClass = StellarClass::Spectral_A;
  282.     StellarClass::LuminosityClass lum = StellarClass::Lum_V;
  283.     unsigned short number = 5;
  284.     int i = 0;
  286.     // Subdwarves (luminosity class VI) are prefixed with sd
  287.     if (starType[i] == 's' && starType[i + 1] == 'd')
  288.     {
  289.         lum = StellarClass::Lum_VI;
  290.         i += 2;
  291.     }
  293.     switch (starType[i])
  294.     {
  295.     case 'O':
  296.         specClass = StellarClass::Spectral_O;
  297.         break;
  298.     case 'B':
  299.         specClass = StellarClass::Spectral_B;
  300.         break;
  301.     case 'A':
  302.         specClass = StellarClass::Spectral_A;
  303.         break;
  304.     case 'F':
  305.         specClass = StellarClass::Spectral_F;
  306.         break;
  307.     case 'G':
  308.         specClass = StellarClass::Spectral_G;
  309.         break;
  310.     case 'K':
  311.         specClass = StellarClass::Spectral_K;
  312.         break;
  313.     case 'M':
  314.         specClass = StellarClass::Spectral_M;
  315.         break;
  316.     case 'R':
  317.         specClass = StellarClass::Spectral_R;
  318.         break;
  319.     case 'N':
  320.         specClass = StellarClass::Spectral_S;
  321.         break;
  322.     case 'S':
  323.         specClass = StellarClass::Spectral_N;
  324.         break;
  325.     case 'W':
  326.         i++;
  327.         if (starType[i] == 'C')
  328.              specClass = StellarClass::Spectral_WC;
  329.         else if (starType[i] == 'N')
  330.             specClass = StellarClass::Spectral_WN;
  331.         else
  332.             i--;
  333.         break;
  335.     case 'D':
  336.         type = StellarClass::WhiteDwarf;
  337.         return StellarClass(type, specClass, 0, lum);
  339.     default:
  340.         specClass = StellarClass::Spectral_Unknown;
  341.         break;
  342.     }
  344.     i++;
  345.     if (starType[i] >= '0' && starType[i] <= '9')
  346.     {
  347.         number = starType[i] - '0';
  348.     }
  349.     else
  350.     {
  351.         // No number given for spectral class; assume it's a 5 unless
  352.         // the star is type O, as O5 stars are exceedingly rare.
  353.         if (specClass == StellarClass::Spectral_O)
  354.             number = 9;
  355.         else
  356.             number = 5;
  357.     }
  359.     if (lum != StellarClass::Lum_VI)
  360.     {
  361.         i++;
  362.         lum = StellarClass::Lum_V;
  363.         while (i < 13 && starType[i] != '') {
  364.             if (starType[i] == 'I') {
  365.                 if (starType[i + 1] == 'I') {
  366.                     if (starType[i + 2] == 'I') {
  367.                         lum = StellarClass::Lum_III;
  368.                     } else {
  369.                         lum = StellarClass::Lum_II;
  370.                     }
  371.                 } else if (starType[i + 1] == 'V') {
  372.                     lum = StellarClass::Lum_IV;
  373.                 } else if (starType[i + 1] == 'a') {
  374.                     if (starType[i + 2] == '0')
  375.                         lum = StellarClass::Lum_Ia0;
  376.                     else
  377.                         lum = StellarClass::Lum_Ia;
  378.                 } else if (starType[i + 1] == 'b') {
  379.                     lum = StellarClass::Lum_Ib;
  380.                 } else {
  381.                     lum = StellarClass::Lum_Ib;
  382.                 }
  383.                 break;
  384.             } else if (starType[i] == 'V') {
  385.                 if (starType[i + 1] == 'I')
  386.                     lum = StellarClass::Lum_VI;
  387.                 else
  388.                     lum = StellarClass::Lum_V;
  389.                 break;
  390.             }
  391.             i++;
  392.         }
  393.     }
  395.     return StellarClass(type, specClass, number, lum);
  396. }
  399. HipparcosStar TheSun()
  400. {
  401.     HipparcosStar star;
  403.     star.HDCatalogNumber = 0;
  404.     star.HIPCatalogNumber = 0;
  405.     star.ascension = 0.0f;
  406.     star.declination = 0.0f;
  407.     star.parallax = 1000000.0f;
  408.     star.appMag = -15.17f;
  409.     star.stellarClass = StellarClass(StellarClass::NormalStar,
  410.                                      StellarClass::Spectral_G, 2,
  411.                                      StellarClass::Lum_V);
  413.     return star;
  414. }
  417. StellarClass guessSpectralType(float colorIndex, float /*absMag*/)
  418. {
  419.     StellarClass::SpectralClass specClass = StellarClass::Spectral_Unknown;
  420.     float subclass = 0.0f;
  422.     if (colorIndex < -0.25f)
  423.     {
  424.         specClass = StellarClass::Spectral_O;
  425.         subclass = (colorIndex + 0.5f) / 0.25f;
  426.     }
  427.     else if (colorIndex < 0.0f)
  428.     {
  429.         specClass = StellarClass::Spectral_B;
  430.         subclass = (colorIndex + 0.25f) / 0.25f;
  431.     }
  432.     else if (colorIndex < 0.25f)
  433.     {
  434.         specClass = StellarClass::Spectral_A;
  435.         subclass = (colorIndex - 0.0f) / 0.25f;
  436.     }
  437.     else if (colorIndex < 0.6f)
  438.     {
  439.         specClass = StellarClass::Spectral_F;
  440.         subclass = (colorIndex - 0.25f) / 0.35f;
  441.     }
  442.     else if (colorIndex < 0.85f)
  443.     {
  444.         specClass = StellarClass::Spectral_G;
  445.         subclass = (colorIndex - 0.6f) / 0.25f;
  446.     }
  447.     else if (colorIndex < 1.4f)
  448.     {
  449.         specClass = StellarClass::Spectral_K;
  450.         subclass = (colorIndex - 0.85f) / 0.55f;
  451.     }
  452.     else
  453.     {
  454.         specClass = StellarClass::Spectral_M;
  455.         subclass = (colorIndex - 1.4f) / 1.0f;
  456.     }
  458.     if (subclass < 0.0f)
  459.         subclass = 0.0f;
  460.     else if (subclass > 1.0f)
  461.         subclass = 1.0f;
  463.     return StellarClass(StellarClass::NormalStar,
  464.                         specClass,
  465.                         (unsigned int) (subclass * 9.99f),
  466.                         StellarClass::Lum_V);
  467. }
  471. static unsigned int okStars, lineno, changes, tested;
  473. bool CheckStarRecord(istream& in)
  474. {
  475.     HipparcosStar star,*hipstar;
  476.     char buf[HipStarRecordLength];
  477.     bool ok=true;
  479.     in.read(buf, TycStarRecordLength);
  480.     lineno++;
  482.     if (sscanf(buf + 210, "%d", &star.HIPCatalogNumber) != 1)
  483.     {
  484.         // Not in Hipparcos, skip it.
  485.         if (verbose>1)
  486.             cout << "Error reading HIPPARCOS catalog number.n";
  487.         return false;
  488.     }
  490.     if (verbose>1)
  491.         cout << "Testing HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  492.     if (!(hipstar=findStar(star.HIPCatalogNumber)))
  493.     {
  494.         cout << "Error finding HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  495.         return false;
  496.     }
  498.     if (hipstar->tycline)
  499.     {
  500.         if (verbose>0)
  501.             cout << "Duplicate Tycho Line for HIP " << star.HIPCatalogNumber << " from Line " << lineno << ", earlier Line at " << hipstar->tycline << " ." << endl;
  502.     }
  503.     else
  504.         hipstar->tycline=lineno;
  506.     tested++;
  507.     sscanf(buf + 309, "%d", &star.HDCatalogNumber);
  509.     if (sscanf(buf + 224, "%f", &star.e_Mag) != 1)
  510.         /* Tycho Database gives no error in for VMag, so we use error on BTmag
  511.            instead.*/
  512.     {
  513.             /* no standard Error given even in BTmag, give it a large value so
  514.                CheckStarRecord() will use the Tycho value if possible */
  515.             star.e_Mag = 1000.0f;
  516.     }
  517.     else if (star.e_Mag >1000.0)
  518.     {
  519.         cout << "Huge BTmag error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  520.     }
  522.     if (sscanf(buf + 41, "%f", &star.appMag) != 1)
  523.     {
  524.         if (verbose>0)
  525.             cout << "Error reading magnitude for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  526.     }
  527.     else if (star.e_Mag < hipstar->e_Mag)
  528.     {
  529.         hipstar->appMag=star.appMag;
  530.         hipstar->e_Mag=star.e_Mag;
  531.         changes++;
  532.         if (verbose > 2)
  533.             cout << "  Change Mag.n";
  534.     }
  536.     float parallaxError = 0.0f;
  537.     if (sscanf(buf + 119, "%f", &parallaxError) != 0)
  538.     {
  539.         if (star.parallax < 0.0f || parallaxError / star.parallax > 1.0f)
  540.             star.parallaxError = 255u;
  541.         else
  542.             star.parallaxError = (uint8) (parallaxError / star.parallax * 200);
  543.     }
  545.     if (sscanf(buf + 79, "%f", &star.parallax) != 1)
  546.     {
  547.         if (verbose>0)
  548.             cout << "Error reading parallax for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  549.         ok=false;
  550.     }
  551.     else if (star.parallax< 0.0)
  552.     {
  553.         if (hipstar->parallax< 0.0)
  554.         {
  555.             if (verbose>0)
  556.                 cout << "Negative parallax for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  557.         ok=false;
  558.         }
  559.     }
  560.     else if (star.parallaxError < hipstar->parallaxError)
  561.     {
  562.         hipstar->parallax=star.parallax;
  563.         hipstar->parallaxError=star.parallaxError;
  564.         changes++;
  565.         if (verbose > 2)
  566.             cout << "  Change Parallax.n";
  567.     }
  570.     if (sscanf(buf + 105, "%f", &star.e_RA) != 1)
  571.     {
  572.             /* no standard Error givenfor Right Ascension , give it a large
  573.                value so original HIPPARCOS value will be used */
  574.             if (verbose>0)
  575.                 cout << "No RA error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  576.             star.e_RA = 2000.0f;
  577.     }
  578.     else if (star.e_RA >1000.0)
  579.     {
  580.         cout << "Huge RA error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  581.     }
  582.     if (sscanf(buf + 112, "%f", &star.e_DE) != 1)
  583.     {
  584.             /* no standard Error given for Declination, give it a large value
  585.                so original HIPPARCOS value will be used. */
  586.             if (verbose>0)
  587.                 cout << "No DE error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  588.             star.e_DE = 2000.0f;
  589.     }
  590.     else if (star.e_DE >1000.0)
  591.     {
  592.         cout << "Huge DE error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  593.     }
  595.     bool coordReadError = false;
  596.     if (sscanf(buf + 51, "%f", &star.ascension) != 1)
  597.         coordReadError = true;
  598.     if (sscanf(buf + 64, "%f", &star.declination) != 1)
  599.         coordReadError = true;
  600.     star.ascension = (float) (star.ascension * 24.0 / 360.0);
  602.     // Read the lower resolution coordinates in hhmmss form if we failed
  603.     // to read the coordinates in degrees.  Not sure why the high resolution
  604.     // coordinates are occasionally missing . . .
  605.     if (coordReadError)
  606.     {
  607.         int hh = 0;
  608.         int mm = 0;
  609.         float seconds;
  610.         coordReadError=false;
  611.         if (sscanf(buf + 17, "%d %d %f", &hh, &mm, &seconds) != 3)
  612.         {
  613.             cout << "Error reading ascension for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  614.             coordReadError=true;;
  615.         }
  616.         else
  617.         {
  618.             star.ascension = hh + (float) mm / 60.0f + (float) seconds / 3600.0f;
  620.             char decSign;
  621.             int deg;
  622.             if (sscanf(buf + 29, "%c%d %d %f", &decSign, &deg, &mm, &seconds) != 4)
  623.             {
  624.                 cout << "Error reading declination for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  625.                 coordReadError=true;;
  626.             }
  627.             else
  628.             {
  629.                 star.declination = deg + (float) mm / 60.0f + (float) seconds / 3600.0f;
  630.                 if (decSign == '-')
  631.                     star.declination = -star.declination;
  632.             }
  633.         }
  634.     }
  636.     if (!((coordReadError) || ((star.ascension==hipstar->ascension) && (star.declination==hipstar->declination))))
  637.     {
  638.         float ast=star.e_RA * star.e_DE;
  639.         float ahi=hipstar->e_RA * hipstar->e_DE;
  640.         if ((ast<ahi) || ((ast==ahi) && ((star.e_RA + star.e_DE) < (hipstar->e_RA + hipstar->e_DE))))
  641.         //if ((star.e_RA * star.e_DE) < (hipstar->e_RA * hipstar->e_DE))
  642.         {
  643.             // Error on the Tycho value is smaller, use it.
  644.             hipstar->ascension=star.ascension;
  645.             hipstar->declination=star.declination;
  646.             hipstar->e_RA=star.e_RA;
  647.             hipstar->e_DE=star.e_DE;
  648.             changes++;
  649.             if (verbose > 2)
  650.                 cout << "  Change Position.n";
  651.         }
  652.     }
  654.     int asc = 0;
  655.     int dec = 0;
  656.     char decSign = ' ';
  657.     if (sscanf(buf + 299, "%d%c%d", &asc, &decSign, &dec) == 3)
  658.     {
  659.         if (decSign == '-')
  660.             dec = -dec;
  661.         star.CCDMIdentifier = asc << 16 | (dec & 0xffff);
  663.         if (star.CCDMIdentifier != hipstar->CCDMIdentifier)
  664.         {
  665.             cout << "Diffrence in CCDM Identifier for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  666.         ok=false;
  667.         }
  668.     }
  670.     if (ok)
  671.         okStars++;
  673.     return true;
  674. }
  677. bool ReadStarRecord(istream& in)
  678. {
  679.     HipparcosStar star;
  680.     char buf[HipStarRecordLength];
  682.     in.read(buf, HipStarRecordLength);
  684.     if (sscanf(buf + 8, "%d", &star.HIPCatalogNumber) != 1)
  685.     {
  686.         cout << "Error reading catalog number.n";
  687.         return false;
  688.     }
  690.     sscanf(buf + 390, "%d", &star.HDCatalogNumber);
  691.     star.tycline=0;
  693.     if (sscanf(buf + 41, "%f", &star.appMag) != 1)
  694.     {
  695.         if (verbose>0)
  696.             cout << "Error reading magnitude for HIP " << star.HIPCatalogNumber << " ." << endl;
  697.         return false;
  698.     }
  700.     if (sscanf(buf + 79, "%f", &star.parallax) != 1)
  701.     {
  702.         if (verbose>0)
  703.             cout << "Error reading parallax for HIP " << star.HIPCatalogNumber << " ." << endl;
  704.         return false;
  705.     }
  706.     if (star.parallax< 0.0)
  707.     {
  708.         if (verbose>0)
  709.             cout << "Negative parallax for HIP " << star.HIPCatalogNumber << " ." << endl;
  710.     }
  712.     bool coordReadError = false;
  713.     if (sscanf(buf + 51, "%f", &star.ascension) != 1)
  714.         coordReadError = true;
  715.     if (sscanf(buf + 64, "%f", &star.declination) != 1)
  716.         coordReadError = true;
  717.     star.ascension = (float) (star.ascension * 24.0 / 360.0);
  719.     // Read the lower resolution coordinates in hhmmss form if we failed
  720.     // to read the coordinates in degrees.  Not sure why the high resolution
  721.     // coordinates are occasionally missing . . .
  722.     if (coordReadError)
  723.     {
  724.         int hh = 0;
  725.         int mm = 0;
  726.         float seconds;
  727.         if (sscanf(buf + 17, "%d %d %f", &hh, &mm, &seconds) != 3)
  728.         {
  729.             cout << "Error reading ascension for HIP " << star.HIPCatalogNumber << " ." << endl;
  730.             return false;
  731.         }
  732.         star.ascension = hh + (float) mm / 60.0f + (float) seconds / 3600.0f;
  734.         char decSign;
  735.         int deg;
  736.         if (sscanf(buf + 29, "%c%d %d %f", &decSign, &deg, &mm, &seconds) != 4)
  737.         {
  738.             cout << "Error reading declination for HIP " << star.HIPCatalogNumber << " ." << endl;
  739.             return false;
  740.         }
  741.         star.declination = deg + (float) mm / 60.0f + (float) seconds / 3600.0f;
  742.         if (decSign == '-')
  743.             star.declination = -star.declination;
  744.     }
  746.     int asc = 0;
  747.     int dec = 0;
  748.     char decSign = ' ';
  749.     int n = 1;
  750.     if (sscanf(buf + 327, "%d%c%d", &asc, &decSign, &dec) == 3)
  751.     {
  752.         if (decSign == '-')
  753.             dec = -dec;
  754.         star.CCDMIdentifier = asc << 16 | (dec & 0xffff);
  756.         sscanf(buf + 340, "%d", &n);
  757.         star.starsWithCCDM = (uint8) n;
  758.         sscanf(buf + 343, "%d", &n);
  759.         star.nComponents = (uint8) n;
  760.     }
  762.     char* spectralType = buf + 435;
  763.     spectralType[12] = '';
  764.     star.stellarClass = ParseStellarClass(spectralType);
  765.     if (star.stellarClass.getSpectralClass() == StellarClass::Spectral_Unknown)
  766.     {
  767.         float bmag,vmag;
  768.         if ((sscanf(buf + 217, "%f", &bmag) == 1) &&
  769.             (sscanf(buf + 230, "%f", &vmag) == 1))
  770.             {
  771.             if (verbose>0)
  772.                 cout << "Guessing Type " << spectralType << "for HIP " << star.HIPCatalogNumber << " ." << endl;
  773.             star.stellarClass = guessSpectralType(bmag - vmag, 0.0f);
  774.             }
  775.         else if (verbose>0)
  776.             cout << "Unparsable stellar class " << spectralType << "for HIP " << star.HIPCatalogNumber << " ." << endl;
  777.     }
  779.     float parallaxError = 0.0f;
  780.     if (sscanf(buf + 119, "%f", &parallaxError) != 0)
  781.     {
  782.         if (star.parallax < 0.0f || parallaxError / star.parallax > 1.0f)
  783.             star.parallaxError = 255u;
  784.         else
  785.             star.parallaxError = (uint8) (parallaxError / star.parallax * 200);
  786.     }
  788.     if (sscanf(buf + 105, "%f", &star.e_RA) != 1)
  789.     {
  790.             /* no standard Error givenfor Right Ascension , give it a large
  791.                value so CheckStarRecord() will use Tycho value if possible */
  792.             star.e_RA = 1000.0f;
  793.     }
  794.     else if (star.e_RA >1000.0)
  795.     {
  796.         cout << "Huge RA error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  797.     }
  798.     if (sscanf(buf + 112, "%f", &star.e_DE) != 1)
  799.     {
  800.             /* no standard Error given for Declination, give it a large value
  801.                so CheckStarRecord() will use the Tycho value if possible */
  802.             star.e_DE = 1000.0f;
  803.     }
  804.     else if (star.e_DE >1000.0)
  805.     {
  806.         cout << "Huge DE error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  807.     }
  809.     if (sscanf(buf + 224, "%f", &star.e_Mag) != 1)
  810.         // No error in HIPPARCOS for VMag, use error on BTmag instead.
  811.     {
  812.             /* no standard Error given give it a large value so CheckStarRecord() will use the Tycho value if possible */
  813.             star.e_Mag = 1000.0f;
  814.     }
  815.     else if (star.e_Mag >1000.0)
  816.     {
  817.         cout << "Huge BTmag error for HIP " << star.HIPCatalogNumber << " from Line " << lineno << " ." << endl;
  818.     }
  820.     stars.insert(stars.end(), star);
  822.     return true;
  823. }
  826. bool ReadComponentRecord(istream& in)
  827. {
  828.     HipparcosComponent component;
  829.     char buf[HipComponentRecordLength];
  831.     in.read(buf, HipComponentRecordLength);
  833.     uint32 hip;
  834.     if (sscanf(buf + 42, "%ud", &hip) != 1)
  835.     {
  836.         cout << "Missing HIP catalog number for component.n";
  837.         return false;
  838.     }
  840.     component.star = findStar(hip);
  841.     if (component.star == NULL)
  842.     {
  843.         cout << "Nonexistent HIP catalog number for component.n";
  844.         return false;
  845.     }
  847.     if (sscanf(buf + 40, "%c", &component.componentID) != 1)
  848.     {
  849.         cout << "Missing component identifier.n";
  850.         return false;
  851.     }
  853.     if (sscanf(buf + 175, "%c", &component.refComponentID) != 1)
  854.     {
  855.         cout << "Error reading reference component.n";
  856.         return false;
  857.     }
  858.     if (component.refComponentID == ' ')
  859.         component.refComponentID = component.componentID;
  861.     // Read astrometric information
  862.     if (sscanf(buf + 88, "%f", &component.ascension) != 1)
  863.     {
  864.         cout << "Missing ascension for component.n";
  865.         return false;
  866.     }
  867.     component.ascension = (float) (component.ascension * 24.0 / 360.0);
  869.     if (sscanf(buf + 101, "%f", &component.declination) != 1)
  870.     {
  871.         cout << "Missing declination for component.n";
  872.         return false;
  873.     }
  875.     // Read photometric information
  876.     if (sscanf(buf + 49, "%f", &component.appMag) != 1)
  877.     {
  878.         cout << "Missing magnitude for component.n";
  879.         return false;
  880.     }
  882.     // vMag and bMag may be necessary to guess the spectral type
  883.     if (sscanf(buf + 62, "%f", &component.bMag) != 1 ||
  884.         sscanf(buf + 75, "%f", &component.vMag) != 1)
  885.     {
  886.         component.bMag = component.vMag = component.appMag;
  887.     }
  888.     else
  889.     {
  890.         component.hasBV = true;
  891.     }
  893.     if (component.componentID != component.refComponentID)
  894.     {
  895.         if (sscanf(buf + 177, "%f", &component.positionAngle) != 1)
  896.         {
  897.             cout << "Missing position angle for component.n";
  898.             return false;
  899.         }
  900.         if (sscanf(buf + 185, "%f", &component.separation) != 1)
  901.         {
  902.             cout << "Missing separation for component.n";
  903.             return false;
  904.         }
  905.     }
  907.     components.insert(components.end(), component);
  909.     return true;
  910. };
  913. void BuildMultistarSystemCatalog()
  914. {
  915.     for (vector<HipparcosStar>::iterator iter = stars.begin();
  916.          iter != stars.end(); iter++)
  917.     {
  918.         if (iter->CCDMIdentifier != NullCCDMIdentifier)
  919.         {
  920.             MultistarSystemCatalog::iterator it =
  921.                 starSystems.find(iter->CCDMIdentifier);
  922.             if (it == starSystems.end())
  923.             {
  924.                 MultistarSystem* multiSystem = new MultistarSystem();
  925.                 multiSystem->nStars = 1;
  926.                 multiSystem->stars[0] = &*iter;
  927.                 starSystems.insert(MultistarSystemCatalog::value_type(iter->CCDMIdentifier, multiSystem));
  928.             }
  929.             else
  930.             {
  931.                 MultistarSystem* multiSystem = it->second;
  932.                 if (multiSystem->nStars == 4)
  933.                 {
  934.                     cout << "Number of stars in system exceeds 4n";
  935.                 }
  936.                 else
  937.                 {
  938.                     multiSystem->stars[multiSystem->nStars] = &*iter;
  939.                     multiSystem->nStars++;
  940.                 }
  941.             }
  942.         }
  943.     }
  944. }
  947. void ConstrainComponentParallaxes()
  948. {
  949.     for (MultistarSystemCatalog::iterator iter = starSystems.begin();
  950.          iter != starSystems.end(); iter++)
  951.     {
  952.         MultistarSystem* multiSystem = iter->second;
  953.         if (multiSystem->nStars > 1)
  954.         {
  955.             for (int i = 1; i < multiSystem->nStars; i++)
  956.                 multiSystem->stars[i]->parallax = multiSystem->stars[0]->parallax;
  957.         }
  958. #if 0
  959.         if (multiSystem->nStars > 2)
  960.         {
  961.             cout << multiSystem->nStars << ": ";
  962.             if (multiSystem->stars[0]->HDCatalogNumber != NullCatalogNumber)
  963.                 cout << "HD " << multiSystem->stars[0]->HDCatalogNumber;
  964.             else
  965.                 cout << "HIP " << multiSystem->stars[0]->HIPCatalogNumber;
  966.             cout << 'n';
  967.         }
  968. #endif
  969.     }
  970. }
  973. void CorrectErrors()
  974. {
  975.     for (vector<HipparcosStar>::iterator iter = stars.begin();
  976.          iter != stars.end(); iter++)
  977.     {
  978.         // Fix the spectral class of Capella, listed for some reason
  979.         // as M1 in the database.
  980.         if (iter->HDCatalogNumber == 34029)
  981.         {
  982.             iter->stellarClass = StellarClass(StellarClass::NormalStar,
  983.                                               StellarClass::Spectral_G, 0,
  984.                                               StellarClass::Lum_III);
  985.         }
  986.     }
  987. }
  990. // Process the vector of star components and insert those that are companions
  991. // of stars in the primary database into the companions vector.
  992. void CreateCompanionList()
  993. {
  994.     for (vector<HipparcosComponent>::iterator iter = components.begin();
  995.          iter != components.end(); iter++)
  996.     {
  997.         // Don't insert the reference component, as this star should already
  998.         // be in the primary database.
  999.         if (iter->componentID != iter->refComponentID)
  1000.         {
  1001.             int componentNumber = iter->componentID - 'A';
  1002.             if (componentNumber > 0 && componentNumber < 8)
  1003.             {
  1004.                 HipparcosStar star;
  1006.                 star.HDCatalogNumber = NullCatalogNumber;
  1007.                 star.HIPCatalogNumber = iter->star->HIPCatalogNumber |
  1008.                     (componentNumber << 25);
  1010.                 star.ascension = iter->ascension;
  1011.                 star.declination = iter->declination;
  1012.                 star.parallax = iter->star->parallax;
  1013.                 star.appMag = iter->appMag;
  1014.                 if (iter->hasBV)
  1015.                     star.stellarClass = guessSpectralType(iter->bMag - iter->vMag, 0.0f);
  1016.                 else
  1017.                     star.stellarClass = StellarClass(StellarClass::NormalStar,
  1018.                                                      StellarClass::Spectral_Unknown,
  1019.                                                      0, StellarClass::Lum_V);
  1021.                 star.CCDMIdentifier = iter->star->CCDMIdentifier;
  1022.                 star.parallaxError = iter->star->parallaxError;
  1024.                 companions.insert(companions.end(), star);
  1025.             }
  1026.         }
  1027.     }
  1028. }
  1031. void ShowStarsWithComponents()
  1032. {
  1033.     cout << "nStars with >2 componentsn";
  1034.     for (vector<HipparcosStar>::iterator iter = stars.begin();
  1035.          iter != stars.end(); iter++)
  1036.     {
  1037.         if (iter->nComponents > 2)
  1038.         {
  1039.             cout << (int) iter->nComponents << ": ";
  1040.             if (iter->HDCatalogNumber != NullCatalogNumber)
  1041.                 cout << "HD " << iter->HDCatalogNumber;
  1042.             else
  1043.                 cout << "HIP " << iter->HIPCatalogNumber;
  1044.             cout << 'n';
  1045.         }
  1046.     }
  1047. }
  1050. void CompareTycho()
  1051. {
  1052.     ifstream tycDatabase(TychoDatabaseFile.c_str(), ios::in | ios::binary);
  1053.     if (!tycDatabase.is_open())
  1054.     {
  1055.         cout << "Error opening " << TychoDatabaseFile << 'n';
  1056.         cout << "You may download this file from ftp://cdsarc.u-strasbg.fr/cats/I/239/n";
  1057.         return;
  1058.     }
  1060.     int recs=0;
  1061.     cout << "Comparing Tycho data set.n";
  1062.     okStars=0;
  1063.     lineno=0;
  1064.     changes=0;
  1065.     while (tycDatabase.good())
  1066.     {
  1067.         CheckStarRecord(tycDatabase);
  1068.         if (++recs % 10000 == 0)
  1069.             {
  1070.             if (verbose>=0)
  1071.                 cout << recs << " records.n";
  1072.             else
  1073.                 cout << ".";
  1074.                 cout.flush();
  1075.             }
  1076.     }
  1077.     if (verbose<0)
  1078.         cout << "n";
  1079.     else
  1080.         cout << recs  << " records checked, " << tested << " tested,  " << okStars << " checked out OK, and " << changes << " changes were made.n";
  1081. }
  1085. int main(int argc, char* argv[])
  1086. {
  1087.     assert(sizeof(StellarClass) == 2);
  1088.     verbose=0;
  1089.     dropstars=1;
  1090.     int c;
  1091.     while((c=getopt(argc,argv,"v::qd:"))>-1)
  1092.     {
  1093.         if (c=='?')
  1094.         {
  1095.             cout << "Usage: buildstardb [-v[<verbosity_level>] [-q] [-d <drop-level>n";
  1096.             exit(1);
  1097.         }
  1098.         else if (c=='v')
  1099.         {
  1100.             if (optarg)
  1101.                 {
  1102.                 verbose=(int)atol(optarg);
  1103.                 if (verbose<-1)
  1104.                     verbose=-1;
  1105.                 else if (verbose>3)
  1106.                     verbose=3;
  1107.                 }
  1108.             else
  1109.                 verbose=1;
  1110.         }
  1111.         else if (c=='d')
  1112.         {   /* Dropstar level. 0 = don't drop stars
  1113.                                1 = drop only non-naked eye visible stars
  1114.                                2 = drop all stars with strange values */
  1115.             dropstars=(int)atol(optarg);
  1116.             if (dropstars<0)
  1117.                 dropstars=0;
  1118.             else if (dropstars>2)
  1119.                 dropstars=2;
  1120.         }
  1121.         else if (c=='q')
  1122.         {
  1123.             verbose=-1;
  1124.         }
  1125.     }
  1127.     // Read star records from the primary HIPPARCOS catalog
  1128.     {
  1129.         ifstream mainDatabase(MainDatabaseFile.c_str(), ios::in | ios::binary);
  1130.         if (!mainDatabase.is_open())
  1131.         {
  1132.             cout << "Error opening " << MainDatabaseFile << 'n';
  1133.             cout << "You may download this file from ftp://cdsarc.u-strasbg.fr/cats/I/239/n";
  1134.             exit(1);
  1135.         }
  1137.         cout << "Reading HIPPARCOS data set.n";
  1138.         while (mainDatabase.good())
  1139.         {
  1140.             ReadStarRecord(mainDatabase);
  1141.             if (stars.size() % 10000 == 0)
  1142.                 {
  1143.                 if (verbose>=0)
  1144.                     cout << stars.size() << " records.n";
  1145.                 else
  1146.                     cout << ".";
  1147.                     cout.flush();
  1148.                 }
  1149.         }
  1150.         if (verbose<0)
  1151.             cout << "n";
  1152.     }
  1153.     if (verbose>=0)
  1154.     {
  1155.         cout << "Read " << stars.size() << " stars from main database.n";
  1156.         cout << "Adding the Sun...n";
  1157.     }
  1158.     stars.insert(stars.end(), TheSun());
  1160.     if (verbose>=0)
  1161.         cout << "Sorting stars...n";
  1162.     {
  1163.         starIndex.reserve(stars.size());
  1164.         for (vector<HipparcosStar>::iterator iter = stars.begin();
  1165.              iter != stars.end(); iter++)
  1166.         {
  1167.             starIndex.insert(starIndex.end(), &*iter);
  1168.         }
  1170.         HIPCatalogComparePredicate pred;
  1172.         // It may not even be necessary to sort the records, if the
  1173.         // HIPPARCOS catalog is strictly ordered by catalog number.  I'm not
  1174.         // sure about this however,
  1175.         random_shuffle(starIndex.begin(), starIndex.end());
  1176.         sort(starIndex.begin(), starIndex.end(), pred);
  1177.     }
  1179.     // Read component records
  1180.     {
  1181.         ifstream componentDatabase(ComponentDatabaseFile.c_str(),
  1182.                                    ios::in | ios::binary);
  1183.         if (!componentDatabase.is_open())
  1184.         {
  1185.             cout << "Error opening " << ComponentDatabaseFile << 'n';
  1186.             cout << "You may download this file from ftp://cdsarc.u-strasbg.fr/cats/I/239/n";
  1187.             exit(1);
  1188.         }
  1190.         if (verbose>=0)
  1191.             cout << "Reading HIPPARCOS component database.n";
  1192.         while (componentDatabase.good())
  1193.         {
  1194.             ReadComponentRecord(componentDatabase);
  1195.         }
  1196.     }
  1197.     if (verbose>=0)
  1198.         cout << "Read " << components.size() << " components.n";
  1200.     {
  1201.         int aComp = 0, bComp = 0, cComp = 0, dComp = 0, eComp = 0, otherComp = 0;
  1202.         int bvComp = 0;
  1203.         for (unsigned int i = 0; i < components.size(); i++)
  1204.         {
  1205.             switch (components[i].componentID)
  1206.             {
  1207.             case 'A':
  1208.                 aComp++; break;
  1209.             case 'B':
  1210.                 bComp++; break;
  1211.             case 'C':
  1212.                 cComp++; break;
  1213.             case 'D':
  1214.                 dComp++; break;
  1215.             case 'E':
  1216.                 eComp++; break;
  1217.             default:
  1218.                 otherComp++; break;
  1219.             }
  1220.             if (components[i].hasBV && components[i].componentID != 'A')
  1221.                 bvComp++;
  1222.         }
  1224.         if (verbose>=0)
  1225.         {
  1226.             cout << "A:" << aComp << "  B:" << bComp << "  C:" << cComp << "  D:" << dComp << "  E:" << eComp << 'n';
  1227.             cout << "Components with B-V mag: " << bvComp << 'n';
  1228.         }
  1229.     }
  1231.     if (verbose>=0)
  1232.         cout << "Building catalog of multiple star systems.n";
  1233.     BuildMultistarSystemCatalog();
  1235.     int nMultipleSystems = starSystems.size();
  1236.     if (verbose>=0)
  1237.         cout << "Stars in multiple star systems: " << nMultipleSystems << 'n';
  1239.     ConstrainComponentParallaxes();
  1241.     CorrectErrors();
  1243.     CompareTycho();
  1245.     // CreateCompanionList();
  1246.     if (verbose>=0)
  1247.     {
  1248.         cout << "Companion stars: " << companions.size() << 'n';
  1249.         cout << "Total stars: " << stars.size() + companions.size() << 'n';
  1250.     }
  1252.     if (verbose>0)
  1253.         ShowStarsWithComponents();
  1255.     const char* outputFile = "stars.dat";
  1256.     if (argv[optind])
  1257.         outputFile = argv[optind];
  1259.     cout << "Writing processed star records to " << outputFile << 'n';
  1260.     ofstream out(outputFile, ios::out | ios::binary);
  1261.     if (!out.good())
  1262.     {
  1263.         cout << "Error opening " << outputFile << 'n';
  1264.         exit(1);
  1265.     }
  1267.     s_er=0.0; // Zero the statistics values
  1268.     s_erq=0.0;
  1269.     n_er=0;
  1270.     s_er=0.0;
  1271.     s_erq=0.0;
  1272.     n_er=0;
  1273.     n_drop=0;
  1274.     n_dub=0;
  1275.     {
  1276.         vector<HipparcosStar>::iterator iter;
  1277.         for (iter = stars.begin(); iter != stars.end(); iter++)
  1278.             iter->analyze();
  1279.         for (iter = companions.begin(); iter != companions.end(); iter++)
  1280.             iter->analyze();
  1281.         binwrite(out, stars.size() + companions.size() - n_drop);
  1282.         float av_r,av_d; // average Right Ascension/Declination
  1283.         av_r=s_er/((float)n_er);
  1284.         av_d=s_ed/((float)n_ed);
  1285.         if (verbose>=0)
  1286.         {
  1287.             cout << "RA Error average: " << av_r << " with Standard Error: " << sqrt((s_erq+(square(s_er)/n_er) - (2*av_r*s_er))/(n_er-1)) << " .n";
  1288.             cout << "DE Error average: " << av_d << " with Standard Error: " << sqrt((s_edq+(square(s_ed)/n_ed) - (2*av_d*s_ed))/(n_ed-1)) << " .n";
  1289.         }
  1290.         for (iter = stars.begin(); iter != stars.end(); iter++)
  1291.             iter->write(out);
  1292.         for (iter = companions.begin(); iter != companions.end(); iter++)
  1293.             iter->write(out);
  1294.     }
  1295.     cout << "Stars processed: " << stars.size() + companions.size() << "   Number dropped: " << n_drop << "  number dubious: " << n_dub << " .n";
  1297. #if 0
  1298.     char* hdOutputFile = "hdxref.dat";
  1300.     cout << "Writing out HD cross reference to " << hdOutputFile << 'n';
  1301.     ofstream hdout(hdOutputFile, ios::out | ios::binary);
  1302.     if (!out.good())
  1303.     {
  1304.         cout << "Error opening " << hdOutputFile << 'n';
  1305.         exit(1);
  1306.     }
  1308.     {
  1309.         int nHD = 0;
  1310.         vector<HipparcosStar>::iterator iter;
  1311.         for (iter = stars.begin(); iter != stars.end(); iter++)
  1312.         {
  1313.             if (iter->HDCatalogNumber != NullCatalogNumber)
  1314.                 nHD++;
  1315.         }
  1316.         binwrite(hdout, nHD);
  1318.         cout << nHD << " stars have HD numbers.n";
  1320.         for (iter = stars.begin(); iter != stars.end(); iter++)
  1321.         {
  1322.             if (iter->HDCatalogNumber != NullCatalogNumber)
  1323.             {
  1324.                 binwrite(hdout, iter->HDCatalogNumber);
  1325.                 binwrite(hdout, iter->HIPCatalogNumber);
  1326.             }
  1327.         }
  1328.     }
  1329. #endif
  1331.     return 0;
  1332. }