Merge lp:~stellarium/stellarium/catalog-cleanup into lp:stellarium

[Alexander Wolf]

This branch contains re-packed stars catalogues without a many artifacts.

I prefer merge it by myself.

[treaves (treaves)]

Looks O.K.

I t would be nice to wrap this in a gui at some point.

review: Approve

=== modified file 'src/core/modules/StarMgr.cpp'
--- src/core/modules/StarMgr.cpp	2013-10-18 04:42:17 +0000
+++ src/core/modules/StarMgr.cpp	2013-11-22 14:47:49 +0000
@@ -68,7 +68,7 @@
 // This number must be incremented each time the content or file format of the stars catalogs change
 // It can also be incremented when the defaultStarsConfig.json file change.
 // It should always matchs the version field of the defaultStarsConfig.json file
-static const int StarCatalogFormatVersion = 4;
+static const int StarCatalogFormatVersion = 5;
 
 // Initialise statics
 bool StarMgr::flagSciNames = true;

=== modified file 'stars/default/CMakeLists.txt'
--- stars/default/CMakeLists.txt	2013-05-08 15:48:03 +0000
+++ stars/default/CMakeLists.txt	2013-11-22 14:47:49 +0000
@@ -1,5 +1,5 @@
 
 ########### install files ###############
 
-INSTALL(FILES defaultStarsConfig.json name.fab stars_0_0v0_3.cat stars_1_0v0_3.cat stars_2_0v0_3.cat stars_3_1v0_2.cat stars_hip_cids_0v0_0.cat stars_hip_sp_0v0_0.cat gcvs_hip_part.dat DESTINATION share/${PACKAGE}/stars/default)
+INSTALL(FILES defaultStarsConfig.json name.fab stars_0_0v0_4.cat stars_1_0v0_4.cat stars_2_0v0_4.cat stars_3_1v0_3.cat stars_hip_cids_0v0_0.cat stars_hip_sp_0v0_0.cat gcvs_hip_part.dat DESTINATION share/${PACKAGE}/stars/default)
 

=== modified file 'stars/default/defaultStarsConfig.json'
--- stars/default/defaultStarsConfig.json	2013-05-03 17:35:56 +0000
+++ stars/default/defaultStarsConfig.json	2013-11-22 14:47:49 +0000
@@ -1,93 +1,93 @@
 {
-	"version": 4,
+	"version": 5,
 	"hipSpectralFile": "stars_hip_sp_0v0_0.cat",
 	"hipComponentsIdsFile": "stars_hip_cids_0v0_0.cat",
 	"catalogs":
 	[
 		{
 			"id": "stars0",
-			"fileName": "stars_0_0v0_3.cat",
+			"fileName": "stars_0_0v0_4.cat",
 			"count": 0.005,
 			"magRange": [-2, 6],
 			"sizeMb": 0.1,
-			"checksum": "89b6833280270f0846f2c7dae1e310e0",
+			"checksum": "d18dafce604d850efeab79aa8402222c",
 			"checked": true
 		},
 		{
 			"id": "stars1",
-			"fileName": "stars_1_0v0_3.cat",
+			"fileName": "stars_1_0v0_4.cat",
 			"count": 0.022,
 			"magRange": [6, 7.5],
 			"sizeMb": 0.6,
-			"checksum": "59eeeb88a0939930a7ae1042fb542bdd",
+			"checksum": "cdf19a423ed3e129bb5cd6fa0d3a09e2",
 			"checked": true
 		},
 		{
 			"id": "stars2",
-			"fileName": "stars_2_0v0_3.cat",
+			"fileName": "stars_2_0v0_4.cat",
 			"count": 0.15,
 			"magRange": [7.5, 9],
 			"sizeMb": 4.1,
-			"checksum": "43e4a39a7198761e675b0f490de8a9ae",
+			"checksum": "5a4c185f79c31c13b0e27d761e407f97",
 			"checked": true
 		},
 		{
 			"id": "stars3",
-			"fileName": "stars_3_1v0_2.cat",
+			"fileName": "stars_3_1v0_3.cat",
 			"count": 0.43,
 			"magRange": [9, 10.5],
-			"sizeMb": 4.2,
-			"checksum": "7bfebdea2d6fd7eae8958f36edd82a77",
+			"sizeMb": 4.1,
+			"checksum": "0b2481f67b0af30dbb027d62761cc061",
 			"checked": true
 		},
 		{
 			"id": "stars4",
-			"fileName": "stars_4_1v0_0.cat",
+			"fileName": "stars_4_1v0_1.cat",
 			"count": 1.7,
 			"magRange": [10.5, 12],
 			"sizeMb": 17,
-			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_4_1v0_0.cat",
-			"checksum": "3abe460ce558418dbdab9a251c173cc4",
+			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_4_1v0_1.cat",
+			"checksum": "7a3b8b3471a755661b6d3686eb4c7cc2",
 			"checked": false
 		},
 		{
 			"id": "stars5",
-			"fileName": "stars_5_2v0_0.cat",
-			"count": 7.7,
+			"fileName": "stars_5_2v0_1.cat",
+			"count": 7.1,
 			"magRange": [12, 13.5],
-			"sizeMb": 44,
-			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_5_2v0_0.cat",
-			"checksum": "67278f69d584b95710eba3e322e22f2e",
+			"sizeMb": 41,
+			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_5_2v0_1.cat",
+			"checksum": "50e798428d9402766805da5cf0e8e819",
 			"checked": false
 		},
 		{
 			"id": "stars6",
-			"fileName": "stars_6_2v0_0.cat",
-			"count": 26.6,
+			"fileName": "stars_6_2v0_1.cat",
+			"count": 24.7,
 			"magRange": [13.5, 15],
-			"sizeMb": 153,
-			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_6_2v0_0.cat",
-			"checksum": "2f9a908226a5b7732319b950729b658b",
+			"sizeMb": 142,
+			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_6_2v0_1.cat",
+			"checksum": "3005959bf7baca36d266def2fecf36d3",
 			"checked": false
 		},
 		{
 			"id": "stars7",
-			"fileName": "stars_7_2v0_0.cat",
-			"count": 57.8,
+			"fileName": "stars_7_2v0_1.cat",
+			"count": 50.8,
 			"magRange": [15, 16.5],
-			"sizeMb": 333,
-			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_7_2v0_0.cat",
-			"checksum": "e05a444abb23232267b10567901b547b",
+			"sizeMb": 292,
+			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_7_2v0_1.cat",
+			"checksum": "15f856fb4c7830a71cb8d60e1decc387",
 			"checked": false
 		},
 		{
 			"id": "stars8",
-			"fileName": "stars_8_2v0_0.cat",
-			"count": 116.9,
+			"fileName": "stars_8_2v0_1.cat",
+			"count": 92.3,
 			"magRange": [16.5, 18],
-			"sizeMb": 674,
-			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_8_2v0_0.cat",
-			"checksum": "b804136f5c5584e9508f21a1f7aa1a28",
+			"sizeMb": 534,
+			"url": "http://downloads.sourceforge.net/sourceforge/stellarium/stars_8_2v0_1.cat",
+			"checksum": "9e2e362022824c60d7e4d94ef8c3af12",
 			"checked": false
 		}
 	]

=== removed file 'stars/default/stars_0_0v0_3.cat'
Binary files stars/default/stars_0_0v0_3.cat	2013-05-03 17:35:56 +0000 and stars/default/stars_0_0v0_3.cat	1970-01-01 00:00:00 +0000 differ
=== added file 'stars/default/stars_0_0v0_4.cat'
Binary files stars/default/stars_0_0v0_4.cat	1970-01-01 00:00:00 +0000 and stars/default/stars_0_0v0_4.cat	2013-11-22 14:47:49 +0000 differ
=== removed file 'stars/default/stars_1_0v0_3.cat'
Binary files stars/default/stars_1_0v0_3.cat	2013-05-03 17:35:56 +0000 and stars/default/stars_1_0v0_3.cat	1970-01-01 00:00:00 +0000 differ
=== added file 'stars/default/stars_1_0v0_4.cat'
Binary files stars/default/stars_1_0v0_4.cat	1970-01-01 00:00:00 +0000 and stars/default/stars_1_0v0_4.cat	2013-11-22 14:47:49 +0000 differ
=== removed file 'stars/default/stars_2_0v0_3.cat'
Binary files stars/default/stars_2_0v0_3.cat	2013-05-03 17:35:56 +0000 and stars/default/stars_2_0v0_3.cat	1970-01-01 00:00:00 +0000 differ
=== added file 'stars/default/stars_2_0v0_4.cat'
Binary files stars/default/stars_2_0v0_4.cat	1970-01-01 00:00:00 +0000 and stars/default/stars_2_0v0_4.cat	2013-11-22 14:47:49 +0000 differ
=== removed file 'stars/default/stars_3_1v0_2.cat'
Binary files stars/default/stars_3_1v0_2.cat	2013-05-03 17:35:56 +0000 and stars/default/stars_3_1v0_2.cat	1970-01-01 00:00:00 +0000 differ
=== added file 'stars/default/stars_3_1v0_3.cat'
Binary files stars/default/stars_3_1v0_3.cat	1970-01-01 00:00:00 +0000 and stars/default/stars_3_1v0_3.cat	2013-11-22 14:47:49 +0000 differ
=== added file 'util/MakeCombinedCatalogue.v2.C'
--- util/MakeCombinedCatalogue.v2.C	1970-01-01 00:00:00 +0000
+++ util/MakeCombinedCatalogue.v2.C	2013-11-22 14:47:49 +0000
@@ -0,0 +1,1690 @@
+// Author and Copyright: Johannes Gajdosik, 2006
+// Modifications: Alexander Wolf, 2013
+// License: GPL
+// g++ -O2 MakeCombinedCatalogue.C -o MakeCombinedCatalogue
+
+
+// change catalogue locations according to your needs:
+
+const char *nomad_names[]={
+  "/storage/1T/nomad.sml/Nomad00.sml",
+  "/storage/1T/nomad.sml/Nomad01.sml",
+  "/storage/1T/nomad.sml/Nomad02.sml",
+  "/storage/1T/nomad.sml/Nomad03.sml",
+  "/storage/1T/nomad.sml/Nomad04.sml",
+  "/storage/1T/nomad.sml/Nomad05.sml",
+  "/storage/1T/nomad.sml/Nomad06.sml",
+  "/storage/1T/nomad.sml/Nomad07.sml",
+  "/storage/1T/nomad.sml/Nomad08.sml",
+  "/storage/1T/nomad.sml/Nomad09.sml",
+  "/storage/1T/nomad.sml/Nomad10.sml",
+  "/storage/1T/nomad.sml/Nomad11.sml",
+  "/storage/1T/nomad.sml/Nomad12.sml",
+  "/storage/1T/nomad.sml/Nomad13.sml",
+  "/storage/1T/nomad.sml/Nomad14.sml",
+  "/storage/1T/nomad.sml/Nomad15.sml",
+  "/storage/1T/nomad.sml/Nomad16.sml",
+  0
+};
+
+
+
+#define NR_OF_HIP 120416
+
+#define MAX_HIP_LEVEL 2
+#define MAX_TYC_LEVEL 4
+#define MAX_LEVEL 9
+
+static
+const char *output_file_names[MAX_LEVEL+1] = {
+  "stars0.cat",
+  "stars1.cat",
+  "stars2.cat",
+  "stars3.cat",
+  "stars4.cat",
+  "stars5.cat",
+  "stars6.cat",
+  "stars7.cat",
+  "stars8.cat",
+  "stars9.cat"
+};
+
+static
+const char *component_ids_filename = "stars_hip_component_ids.cat";
+
+static
+const char *sp_filename = "stars_hip_sp.cat";
+
+static
+const double level_mag_limit[MAX_LEVEL+1] = {
+  6,7.5,9,
+  10.5,12,
+  13.5,15,16.5,18,19.5
+};
+
+
+
+
+
+
+
+
+
+/*
+
+Vector (simple 3d-Vector)
+
+Author and Copyright: Johannes Gajdosik, 2006
+
+License: GPL
+
+*/
+
+
+#include <math.h>
+
+struct Vector {
+  double x[3];
+  inline double length(void) const;
+  inline void normalize(void);
+  const Vector &operator+=(const Vector &a) {
+    x[0] += a.x[0];
+    x[1] += a.x[1];
+    x[2] += a.x[2];
+  }
+  const Vector &operator-=(const Vector &a) {
+    x[0] -= a.x[0];
+    x[1] -= a.x[1];
+    x[2] -= a.x[2];
+  }
+};
+
+static inline
+Vector operator^(const Vector &a,const Vector &b) {
+  Vector c;
+  c.x[0] = a.x[1]*b.x[2] - a.x[2]*b.x[1];
+  c.x[1] = a.x[2]*b.x[0] - a.x[0]*b.x[2];
+  c.x[2] = a.x[0]*b.x[1] - a.x[1]*b.x[0];
+  return c;
+}
+
+static inline
+Vector operator-(const Vector &a,const Vector &b) {
+  Vector c;
+  c.x[0] = a.x[0]-b.x[0];
+  c.x[1] = a.x[1]-b.x[1];
+  c.x[2] = a.x[2]-b.x[2];
+  return c;
+}
+
+static inline
+Vector operator+(const Vector &a,const Vector &b) {
+  Vector c;
+  c.x[0] = a.x[0]+b.x[0];
+  c.x[1] = a.x[1]+b.x[1];
+  c.x[2] = a.x[2]+b.x[2];
+  return c;
+}
+
+static inline
+Vector operator*(double a,const Vector &b) {
+  Vector c;
+  c.x[0] = a*b.x[0];
+  c.x[1] = a*b.x[1];
+  c.x[2] = a*b.x[2];
+  return c;
+}
+
+static inline
+Vector operator*(const Vector &b,double a) {
+  Vector c;
+  c.x[0] = a*b.x[0];
+  c.x[1] = a*b.x[1];
+  c.x[2] = a*b.x[2];
+  return c;
+}
+
+static inline
+double operator*(const Vector &a,const Vector &b) {
+  return a.x[0]*b.x[0]+a.x[1]*b.x[1]+a.x[2]*b.x[2];
+}
+
+double Vector::length(void) const {
+  return sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
+}
+
+void Vector::normalize(void) {
+  const double l = 1.0/length();
+  x[0] *= l;
+  x[1] *= l;
+  x[2] *= l;
+}
+
+
+
+
+void ChangeEpoch(double delta_years,
+                 double &ra,double &dec, // degrees
+                 double pm_ra,double pm_dec // mas/yr
+                ) {
+  ra *= (M_PI/180);
+  dec *= (M_PI/180);
+
+  const double cdec = cos(dec);
+  Vector x = {cos(ra)*cdec,sin(ra)*cdec,sin(dec)};
+  const Vector north = {0.0,0.0,1.0};
+
+  Vector axis0 = north ^ x;
+  axis0.normalize();
+  const Vector axis1 = x ^ axis0;
+
+  const double f = delta_years*(0.001/3600)*(M_PI/180);
+
+  x += pm_ra*f*axis1 + pm_dec*f*axis0;
+
+  ra = atan2(x.x[1],x.x[0]);
+  if (ra < 0.0) ra += 2*M_PI;
+  dec = atan2(x.x[2],sqrt(x.x[0]*x.x[0]+x.x[1]*x.x[1]));
+
+  ra *= (180/M_PI);
+  dec *= (180/M_PI);
+}
+
+
+
+
+
+
+//------------------------------------------------
+
+/*
+
+GeodesicGrid: a library for dividing the sphere into triangle zones
+by subdividing the icosahedron
+
+Author and Copyright: Johannes Gajdosik, 2006
+
+This library requires a simple Vector library,
+which may have different copyright and license.
+
+In the moment I choose to distribute the library under the GPL,
+later I may choose to additionally distribute it under a more
+relaxed license like the LGPL. If you want to have the library
+under another license, please ask me.
+
+
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+using namespace std;
+
+struct HalfSpace {
+    // all x with x*n-d >= 0
+  Vector n;
+  double d;
+  bool inside(const Vector &x) const {return (x*n>=d);}
+};
+
+class GeodesicGrid {
+    // Grid of triangles (zones) on the sphere with radius 1,
+    // generated by subdividing the icosahedron.
+    // level 0: just the icosahedron, 20 zones
+    // level 1: 80 zones, level 2: 320 zones, ...
+    // Each zone has a unique integer number in the range
+    // [0,getNrOfZones()-1].
+public:
+  GeodesicGrid(int max_level);
+  ~GeodesicGrid(void);
+  int getMaxLevel(void) const {return max_level;}
+  static int nrOfZones(int level)
+    {return (20<<(level<<1));} // 20*4^level
+  int getNrOfZones(void) const {return nrOfZones(max_level);}
+  void getTriangleCorners(int lev, int index, Vector& c0, Vector& c1, Vector& c2) const;
+    // Return the position of the 3 corners for the triangle at the given level and index
+  int getPartnerTriangle(int lev, int index) const;
+    // Return the index of the partner triangle with which to form a paralelogram
+  typedef void (VisitFunc)(int lev,int index,
+                           const Vector &c0,
+                           const Vector &c1,
+                           const Vector &c2,
+                           void *context);
+  void visitTriangles(int max_visit_level,
+                      VisitFunc *func,void *context) const;
+
+  int searchZone(const Vector &v,int search_level) const;
+    // find the zone number in which a given point lies.
+    // prerequisite: v*v==1
+    // When the point lies on the border of two or more zones,
+    // one such zone is returned (always the same one,
+    // because the algorithm is deterministic).
+
+  void searchZones(const HalfSpace *half_spaces,int nr_of_half_spaces,
+                   int **inside,int **border,int max_search_level) const;
+    // find all zones that lie fully(inside) or partly(border)
+    // in the intersection of the given half spaces.
+    // The result is accurate when (0,0,0) lies on the border of
+    // each half space. If this is not the case,
+    // the result may be inaccurate, because it is assumed, that
+    // a zone lies in a half space when its 3 corners lie in
+    // this half space.
+    // inside[l] points to the begin of an integer array of size nrOfZones(l),
+    // border[l] points to one after the end of the same integer array.
+    // The array will be filled from the beginning with the inside zone numbers
+    // and from the end with the border zone numbers of the given level l
+    // for 0<=l<=getMaxLevel().
+    // inside[l] will not contain zones that are already contained
+    // in inside[l1] for some l1 < l.
+    // In order to restrict search depth set max_search_level < max_level,
+    // for full search depth set max_search_level = max_level,
+
+private:
+  const Vector& getTriangleCorner(int lev, int index, int cornerNumber) const;
+  void initTriangle(int lev,int index,
+                    const Vector &c0,
+                    const Vector &c1,
+                    const Vector &c2);
+  void visitTriangles(int lev,int index,
+                      const Vector &c0,
+                      const Vector &c1,
+                      const Vector &c2,
+                      int max_visit_level,
+                      VisitFunc *func,
+                      void *context) const;
+  void searchZones(int lev,int index,
+                   const HalfSpace *half_spaces,
+                   const int nr_of_half_spaces,
+                   const int *index_of_used_half_spaces,
+                   const int half_spaces_used,
+                   const bool *corner0_inside,
+                   const bool *corner1_inside,
+                   const bool *corner2_inside,
+                   int **inside,int **border,int max_search_level) const;
+private:
+  const int max_level;
+  struct Triangle {
+    Vector e0,e1,e2;   // Seitenmittelpunkte
+  };
+  Triangle **triangles;
+  // 20*(4^0+4^1+...+4^n)=20*(4*(4^n)-1)/3 triangles total
+  // 2+10*4^n corners
+};
+
+class GeodesicSearchResult {
+public:
+  GeodesicSearchResult(const GeodesicGrid &grid);
+  ~GeodesicSearchResult(void);
+  void search(const HalfSpace *half_spaces,
+              const int nr_of_half_spaces,
+              int max_search_level);
+  void search(const Vector &e0,const Vector &e1,const Vector &e2,const Vector &e3,
+              int max_search_level);
+  void print(void) const;
+private:
+  friend class GeodesicSearchInsideIterator;
+  friend class GeodesicSearchBorderIterator;
+  const GeodesicGrid &grid;
+  int **const zones;
+  int **const inside;
+  int **const border;
+};
+
+class GeodesicSearchBorderIterator {
+public:
+  GeodesicSearchBorderIterator(const GeodesicSearchResult &r,int level)
+    : r(r),level((level<0)?0:(level>r.grid.getMaxLevel())
+                          ?r.grid.getMaxLevel():level),
+      end(r.zones[GeodesicSearchBorderIterator::level]+
+                  GeodesicGrid::nrOfZones(GeodesicSearchBorderIterator::level))
+      {reset();}
+  void reset(void) {index = r.border[level];}
+  int next(void) // returns -1 when finished
+    {if (index < end) {return *index++;} return -1;}
+private:
+  const GeodesicSearchResult &r;
+  const int level;
+  const int *const end;
+  const int *index;
+};
+
+
+class GeodesicSearchInsideIterator {
+public:
+  GeodesicSearchInsideIterator(const GeodesicSearchResult &r,int level)
+    : r(r),max_level((level<0)?0:(level>r.grid.getMaxLevel())
+                              ?r.grid.getMaxLevel():level)
+     {reset();}
+  void reset(void);
+  int next(void); // returns -1 when finished
+private:
+  const GeodesicSearchResult &r;
+  const int max_level;
+  int level;
+  int max_count;
+  int *index_p;
+  int *end_p;
+  int index;
+  int count;
+};
+
+#include <stdlib.h>
+#include <cassert>
+#include <iostream>
+
+using namespace std;
+
+static const double icosahedron_G = 0.5*(1.0+sqrt(5.0));
+static const double icosahedron_b = 1.0/sqrt(1.0+icosahedron_G*icosahedron_G);
+static const double icosahedron_a = icosahedron_b*icosahedron_G;
+
+static const Vector icosahedron_corners[12] = {
+  { icosahedron_a, -icosahedron_b,            0.0},
+  { icosahedron_a,  icosahedron_b,            0.0},
+  {-icosahedron_a,  icosahedron_b,            0.0},
+  {-icosahedron_a, -icosahedron_b,            0.0},
+  {           0.0,  icosahedron_a, -icosahedron_b},
+  {           0.0,  icosahedron_a,  icosahedron_b},
+  {           0.0, -icosahedron_a,  icosahedron_b},
+  {           0.0, -icosahedron_a, -icosahedron_b},
+  {-icosahedron_b,            0.0,  icosahedron_a},
+  { icosahedron_b,            0.0,  icosahedron_a},
+  { icosahedron_b,            0.0, -icosahedron_a},
+  {-icosahedron_b,            0.0, -icosahedron_a}
+};
+
+struct TopLevelTriangle {
+  int corners[3];   // index der Ecken
+};
+
+
+static const TopLevelTriangle icosahedron_triangles[20] = {
+  { 1, 0,10}, //  1
+  { 0, 1, 9}, //  0
+  { 0, 9, 6}, // 12
+  { 9, 8, 6}, //  9
+  { 0, 7,10}, // 16
+  { 6, 7, 0}, //  6
+  { 7, 6, 3}, //  7
+  { 6, 8, 3}, // 14
+  {11,10, 7}, // 11
+  { 7, 3,11}, // 18
+  { 3, 2,11}, //  3
+  { 2, 3, 8}, //  2
+  {10,11, 4}, // 10
+  { 2, 4,11}, // 19
+  { 5, 4, 2}, //  5
+  { 2, 8, 5}, // 15
+  { 4, 1,10}, // 17
+  { 4, 5, 1}, // 4
+  { 5, 9, 1}, // 13
+  { 8, 9, 5}  //  8
+};
+
+
+
+GeodesicGrid::GeodesicGrid(const int lev)
+             :max_level(lev<0?0:lev) {
+  if (max_level > 0) {
+    triangles = new Triangle*[max_level];
+    int nr_of_triangles = 20;
+    for (int i=0;i<max_level;i++) {
+      triangles[i] = new Triangle[nr_of_triangles];
+      nr_of_triangles *= 4;
+    }
+    for (int i=0;i<20;i++) {
+      const int *const corners = icosahedron_triangles[i].corners;
+      initTriangle(0,i,
+                   icosahedron_corners[corners[0]],
+                   icosahedron_corners[corners[1]],
+                   icosahedron_corners[corners[2]]);
+    }
+  } else {
+    triangles = 0;
+  }
+}
+
+GeodesicGrid::~GeodesicGrid(void) {
+  if (max_level > 0) {
+    for (int i=max_level-1;i>=0;i--) delete[] triangles[i];
+    delete triangles;
+  }
+}
+
+void GeodesicGrid::getTriangleCorners(int lev,int index,
+                                      Vector &h0,
+                                      Vector &h1,
+                                      Vector &h2) const {
+  if (lev <= 0) {
+    const int *const corners = icosahedron_triangles[index].corners;
+    h0 = icosahedron_corners[corners[0]];
+    h1 = icosahedron_corners[corners[1]];
+    h2 = icosahedron_corners[corners[2]];
+  } else {
+    lev--;
+    const int i = index>>2;
+    Triangle &t(triangles[lev][i]);
+    switch (index&3) {
+      case 0: {
+        Vector c0,c1,c2;
+        getTriangleCorners(lev,i,c0,c1,c2);
+        h0 = c0;
+        h1 = t.e2;
+        h2 = t.e1;
+      } break;
+      case 1: {
+        Vector c0,c1,c2;
+        getTriangleCorners(lev,i,c0,c1,c2);
+        h0 = t.e2;
+        h1 = c1;
+        h2 = t.e0;
+      } break;
+      case 2: {
+        Vector c0,c1,c2;
+        getTriangleCorners(lev,i,c0,c1,c2);
+        h0 = t.e1;
+        h1 = t.e0;
+        h2 = c2;
+      } break;
+      case 3:
+        h0 = t.e0;
+        h1 = t.e1;
+        h2 = t.e2;
+        break;
+    }
+  }
+}
+
+int GeodesicGrid::getPartnerTriangle(int lev, int index) const
+{
+	if (lev==0)
+	{
+		assert(index<20);
+		return (index&1) ? index+1 : index-1;
+	}
+	switch(index&7)
+	{
+		case 2:
+		case 6:
+			return index+1;
+		case 3:
+		case 7:
+			return index-1;
+		case 0:
+			return (lev==1) ? index+5 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
+		case 1:
+			return (lev==1) ? index+3 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
+		case 4:
+			return (lev==1) ? index-3 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
+		case 5:
+			return (lev==1) ? index-5 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
+		default:
+			assert(0);
+	}
+}
+
+void GeodesicGrid::initTriangle(int lev,int index,
+                                const Vector &c0,
+                                const Vector &c1,
+                                const Vector &c2) {
+  Triangle &t(triangles[lev][index]);
+  t.e0 = c1+c2;
+  t.e0.normalize();
+  t.e1 = c2+c0;
+  t.e1.normalize();
+  t.e2 = c0+c1;
+  t.e2.normalize();
+  lev++;
+  if (lev < max_level) {
+    index *= 4;
+    initTriangle(lev,index+0,c0,t.e2,t.e1);
+    initTriangle(lev,index+1,t.e2,c1,t.e0);
+    initTriangle(lev,index+2,t.e1,t.e0,c2);
+    initTriangle(lev,index+3,t.e0,t.e1,t.e2);
+  }
+}
+
+
+void GeodesicGrid::visitTriangles(int max_visit_level,
+                                  VisitFunc *func,
+                                  void *context) const {
+  if (func && max_visit_level >= 0) {
+    if (max_visit_level > max_level) max_visit_level = max_level;
+    for (int i=0;i<20;i++) {
+      const int *const corners = icosahedron_triangles[i].corners;
+      visitTriangles(0,i,
+                     icosahedron_corners[corners[0]],
+                     icosahedron_corners[corners[1]],
+                     icosahedron_corners[corners[2]],
+                     max_visit_level,func,context);
+    }
+  }
+}
+
+void GeodesicGrid::visitTriangles(int lev,int index,
+                                  const Vector &c0,
+                                  const Vector &c1,
+                                  const Vector &c2,
+                                  int max_visit_level,
+                                  VisitFunc *func,
+                                  void *context) const {
+  (*func)(lev,index,c0,c1,c2,context);
+  Triangle &t(triangles[lev][index]);
+  lev++;
+  if (lev <= max_visit_level) {
+    index *= 4;
+    visitTriangles(lev,index+0,c0,t.e2,t.e1,max_visit_level,func,context);
+    visitTriangles(lev,index+1,t.e2,c1,t.e0,max_visit_level,func,context);
+    visitTriangles(lev,index+2,t.e1,t.e0,c2,max_visit_level,func,context);
+    visitTriangles(lev,index+3,t.e0,t.e1,t.e2,max_visit_level,func,context);
+  }
+}
+
+
+int GeodesicGrid::searchZone(const Vector &v,int search_level) const {
+  for (int i=0;i<20;i++) {
+    const int *const corners = icosahedron_triangles[i].corners;
+    const Vector &c0(icosahedron_corners[corners[0]]);
+    const Vector &c1(icosahedron_corners[corners[1]]);
+    const Vector &c2(icosahedron_corners[corners[2]]);
+    if (((c0^c1)*v >= 0.0) && ((c1^c2)*v >= 0.0) && ((c2^c0)*v >= 0.0)) {
+        // v lies inside this icosahedron triangle
+      for (int lev=0;lev<search_level;lev++) {
+        Triangle &t(triangles[lev][i]);
+        i <<= 2;
+        if ((t.e1^t.e2)*v <= 0.0) {
+          // i += 0;
+        } else if ((t.e2^t.e0)*v <= 0.0) {
+          i += 1;
+        } else if ((t.e0^t.e1)*v <= 0.0) {
+          i += 2;
+        } else {
+          i += 3;
+        }
+      }
+      return i;
+    }
+  }
+  cerr << "software error: not found" << endl;
+  exit(-1);
+    // shut up the compiler
+  return -1;
+}
+
+void GeodesicGrid::searchZones(const HalfSpace *half_spaces,
+                               const int nr_of_half_spaces,
+                               int **inside_list,int **border_list,
+                               int max_search_level) const {
+  if (max_search_level < 0) max_search_level = 0;
+  else if (max_search_level > max_level) max_search_level = max_level;
+  int halfs_used[nr_of_half_spaces];
+  for (int h=0;h<nr_of_half_spaces;h++) {halfs_used[h] = h;}
+  bool corner_inside[12][nr_of_half_spaces];
+  for (int h=0;h<nr_of_half_spaces;h++) {
+    const HalfSpace &half_space(half_spaces[h]);
+    for (int i=0;i<12;i++) {
+      corner_inside[i][h] = half_space.inside(icosahedron_corners[i]);
+    }
+  }
+  for (int i=0;i<20;i++) {
+    searchZones(0,i,
+                half_spaces,nr_of_half_spaces,halfs_used,nr_of_half_spaces,
+                corner_inside[icosahedron_triangles[i].corners[0]],
+                corner_inside[icosahedron_triangles[i].corners[1]],
+                corner_inside[icosahedron_triangles[i].corners[2]],
+                inside_list,border_list,max_search_level);
+  }
+}
+
+void GeodesicGrid::searchZones(int lev,int index,
+                               const HalfSpace *half_spaces,
+                               const int nr_of_half_spaces,
+                               const int *index_of_used_half_spaces,
+                               const int half_spaces_used,
+                               const bool *corner0_inside,
+                               const bool *corner1_inside,
+                               const bool *corner2_inside,
+                               int **inside_list,int **border_list,
+                               const int max_search_level) const {
+  int halfs_used[half_spaces_used];
+  int halfs_used_count = 0;
+  for (int h=0;h<half_spaces_used;h++) {
+    const int i = index_of_used_half_spaces[h];
+    if (!corner0_inside[i] && !corner1_inside[i] && !corner2_inside[i]) {
+        // totally outside this HalfSpace
+      return;
+    } else if (corner0_inside[i] && corner1_inside[i] && corner2_inside[i]) {
+        // totally inside this HalfSpace
+    } else {
+        // on the border of this HalfSpace
+      halfs_used[halfs_used_count++] = i;
+    }
+  }
+  if (halfs_used_count == 0) {
+      // this triangle(lev,index) lies inside all halfspaces
+    **inside_list = index;
+    (*inside_list)++;
+  } else {
+    (*border_list)--;
+    **border_list = index;
+    if (lev < max_search_level) {
+      Triangle &t(triangles[lev][index]);
+      lev++;
+      index <<= 2;
+      inside_list++;
+      border_list++;
+      bool edge0_inside[nr_of_half_spaces];
+      bool edge1_inside[nr_of_half_spaces];
+      bool edge2_inside[nr_of_half_spaces];
+      for (int h=0;h<halfs_used_count;h++) {
+        const int i = halfs_used[h];
+        const HalfSpace &half_space(half_spaces[i]);
+        edge0_inside[i] = half_space.inside(t.e0);
+        edge1_inside[i] = half_space.inside(t.e1);
+        edge2_inside[i] = half_space.inside(t.e2);
+      }
+      searchZones(lev,index+0,
+                  half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
+                  corner0_inside,edge2_inside,edge1_inside,
+                  inside_list,border_list,max_search_level);
+      searchZones(lev,index+1,
+                  half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
+                  edge2_inside,corner1_inside,edge0_inside,
+                  inside_list,border_list,max_search_level);
+      searchZones(lev,index+2,
+                  half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
+                  edge1_inside,edge0_inside,corner2_inside,
+                  inside_list,border_list,max_search_level);
+      searchZones(lev,index+3,
+                  half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
+                  edge0_inside,edge1_inside,edge2_inside,
+                  inside_list,border_list,max_search_level);
+    }
+  }
+}
+
+
+
+GeodesicSearchResult::GeodesicSearchResult(const GeodesicGrid &grid)
+                     :grid(grid),
+                      zones(new int*[grid.getMaxLevel()+1]),
+                      inside(new int*[grid.getMaxLevel()+1]),
+                      border(new int*[grid.getMaxLevel()+1]) {
+  for (int i=0;i<=grid.getMaxLevel();i++) {
+    zones[i] = new int[GeodesicGrid::nrOfZones(i)];
+  }
+}
+
+GeodesicSearchResult::~GeodesicSearchResult(void) {
+  for (int i=grid.getMaxLevel();i>=0;i--) {
+    delete zones[i];
+  }
+  delete border;
+  delete inside;
+  delete zones;
+}
+
+void GeodesicSearchResult::search(const HalfSpace *half_spaces,
+                                  const int nr_of_half_spaces,
+                                  int max_search_level) {
+  for (int i=grid.getMaxLevel();i>=0;i--) {
+    inside[i] = zones[i];
+    border[i] = zones[i]+GeodesicGrid::nrOfZones(i);
+  }
+  grid.searchZones(half_spaces,nr_of_half_spaces,inside,border,max_search_level);
+}
+
+void GeodesicSearchResult::search(const Vector &e0,const Vector &e1,
+                                  const Vector &e2,const Vector &e3,
+                                  int max_search_level) {
+  HalfSpace half_spaces[4];
+  half_spaces[0].d = e0*((e1-e0)^(e2-e0));
+  if (half_spaces[0].d > 0) {
+    half_spaces[0].n = e0^e1;
+    half_spaces[1].n = e1^e2;
+    half_spaces[2].n = e2^e3;
+    half_spaces[3].n = e3^e0;
+    half_spaces[0].d = 0.0;
+    half_spaces[1].d = 0.0;
+    half_spaces[2].d = 0.0;
+    half_spaces[3].d = 0.0;
+    search(half_spaces,4,max_search_level);
+  } else {
+    half_spaces[0].n = (e1-e0)^(e2-e0);
+    search(half_spaces,1,max_search_level);
+  }
+}
+
+
+
+void GeodesicSearchInsideIterator::reset(void) {
+  level = 0;
+  max_count = 1<<(max_level<<1); // 4^max_level
+  index_p = r.zones[0];
+  end_p = r.inside[0];
+  index = (*index_p) * max_count;
+  count = (index_p < end_p) ? 0 : max_count;
+}
+
+int GeodesicSearchInsideIterator::next(void) {
+  if (count < max_count) return index+(count++);
+  index_p++;
+  if (index_p < end_p) {
+    index = (*index_p) * max_count;
+    count = 1;
+    return index;
+  }
+  while (level < max_level) {
+    level++;
+    max_count >>= 2;
+    index_p = r.zones[level];
+    end_p = r.inside[level];
+    if (index_p < end_p) {
+      index = (*index_p) * max_count;
+      count = 1;
+      return index;
+    }
+  }
+  return -1;
+}
+
+//------------------------------------------------
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <map>
+#include <list>
+#include <string>
+#include <iostream>
+
+using namespace std;
+
+
+
+static int restrict_output_level_min = 0;
+static int restrict_output_level_max = MAX_LEVEL;
+
+
+
+struct FaintStar {  // 6 byte
+  int x0:18;
+  int x1:18;
+  unsigned int b_v:7;
+  unsigned int mag:5;
+  bool operator<(const FaintStar &s) const {return (mag < s.mag);}
+} __attribute__ ((__packed__)) ;
+
+
+struct TycStar {  // 10 byte
+  int x0:20;
+  int x1:20;
+  int dx0:14;
+  int dx1:14;
+  unsigned int b_v:7;
+  unsigned int mag:5;
+  bool operator<(const TycStar &s) const {return (mag < s.mag);}
+} __attribute__ ((__packed__));
+
+
+struct HipStarCompressed {
+  int hip:24;                  // 17 bits needed
+  unsigned char comp_ids_int;  //  5 bits needed
+  int x0;                      // 32 bits needed
+  int x1;                      // 32 bits needed
+  unsigned char b_v;           //  8 bits needed
+  unsigned char mag;           //  5 bits needed
+  unsigned short sp_int;       // 14 bits needed
+  int dx0,dx1,plx;
+} __attribute__ ((__packed__));
+
+struct HipStar : public HipStarCompressed {
+  string component_ids,sp;
+  void setPlxSp(double plx,const string &sp) {
+    HipStar::plx = (int)floor(0.5+100.0*plx);
+    if (HipStar::plx < 0) HipStar::plx = 0;
+    HipStar::sp = sp;
+  }
+  bool operator<(const HipStar &s) const {return (mag < s.mag);}
+};
+
+
+
+struct ZoneData {
+  virtual ~ZoneData(void) {}
+  virtual int getArraySize(void) const = 0;
+  virtual HipStar *add(int level,
+                       int tyc1,int tyc2,int tyc3,
+                       int hip,const char *component_ids,
+                       double x0,double x1,double dx0,double dx1,
+                       double mag,double b_v,
+                       double plx,const char *sp,
+                       bool &does_not_fit) = 0;
+  void writeInfoToOutput(FILE *f) const;
+  virtual void writeStarsToOutput(FILE *f) = 0;
+  Vector center;
+  Vector axis0;
+  Vector axis1;
+};
+
+struct HipZoneData : public ZoneData {
+  list<HipStar> stars;
+  int getArraySize(void) const {return stars.size();}
+  HipStar *add(int level,
+                       int tyc1,int tyc2,int tyc3,
+                       int hip,const char *component_ids,
+                       double x0,double x1,double dx0,double dx1,
+                       double mag,double b_v,double plx,const char *sp,
+                       bool &does_not_fit);
+  void writeStarsToOutput(FILE *f);
+};
+
+
+static list<HipStar*> hip_index[NR_OF_HIP+1];
+
+static void SqueezeHip(void) {
+    // build lookup maps
+  map<string,int> component_ids_map,sp_map;
+  for (int i=0;i<=NR_OF_HIP;i++) {
+    for (list<HipStar*>::const_iterator it(hip_index[i].begin());
+         it!=hip_index[i].end();it++) {
+      HipStar *h = *it;
+      component_ids_map[h->component_ids]++;
+      sp_map[h->sp]++;
+    }
+  }
+  int component_ids_size = 0;
+  for (map<string,int>::iterator it(component_ids_map.begin());
+       it!=component_ids_map.end();it++,component_ids_size++) {
+    it->second = component_ids_size;
+  }
+  if (component_ids_size >= 32) {
+    cerr << "SqueezeHip: too many component_ids: "
+         << component_ids_size << endl;
+  }
+  int sp_size = 0;
+  for (map<string,int>::iterator it(sp_map.begin());
+       it!=sp_map.end();it++,sp_size++) {
+    it->second = sp_size;
+  }
+  if (sp_size >= 16384) {
+    cerr << "SqueezeHip: too many sp: " << sp_size << endl;
+  }
+    // translate the strings for the hip stars into integers:
+  for (int i=0;i<=NR_OF_HIP;i++) {
+    for (list<HipStar*>::const_iterator it(hip_index[i].begin());
+         it!=hip_index[i].end();it++) {
+      HipStar *h = *it;
+      h->comp_ids_int = component_ids_map[h->component_ids];
+      h->sp_int = sp_map[h->sp];
+    }
+  }
+    // write output files for component_ids and sp
+  FILE *f = fopen(component_ids_filename,"wb");
+  if (f==0) {
+    cout << "fopen(" << component_ids_filename << ") failed" << endl;
+    exit(-1);
+  }
+  for (map<string,int>::const_iterator it(component_ids_map.begin());
+       it!=component_ids_map.end();it++,component_ids_size++) {
+    fprintf(f,"%s\n",it->first.c_str());
+  }
+  fclose(f);
+  f = fopen(sp_filename,"wb");
+  if (f==0) {
+    cout << "fopen(" << sp_filename << ") failed" << endl;
+    exit(-1);
+  }
+  for (map<string,int>::const_iterator it(sp_map.begin());
+       it!=sp_map.end();it++,sp_size++) {
+    fprintf(f,"%s\n",it->first.c_str());
+  }
+  fclose(f);
+}
+
+struct TycZoneData : public ZoneData {
+  list<TycStar> stars;
+  int getArraySize(void) const {return stars.size();}
+  HipStar *add(int level,
+                       int tyc1,int tyc2,int tyc3,
+                       int hip,const char *component_ids,
+                       double x0,double x1,double dx0,double dx1,
+                       double mag,double b_v,double plx,const char *sp,
+                       bool &does_not_fit);
+  void writeStarsToOutput(FILE *f);
+};
+
+struct FaintZoneData : public ZoneData {
+  list<FaintStar> stars;
+  int getArraySize(void) const {return stars.size();}
+  HipStar *add(int level,
+                       int tyc1,int tyc2,int tyc3,
+                       int hip,const char *component_ids,
+                       double x0,double x1,double dx0,double dx1,
+                       double mag,double b_v,double plx,const char *sp,
+                       bool &does_not_fit);
+  void writeStarsToOutput(FILE *f);
+};
+
+
+class Accumulator {
+public:
+  Accumulator(void);
+  ~Accumulator(void);
+  int addStar(int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
+              double ra, // degrees
+              double dec, // degrees
+              double pma,double pmd,
+              double mag,double b_v,
+              double plx,const char *sp);
+  void writeOutput(const char *fnames[]);
+private:
+  GeodesicGrid grid;
+  struct ZoneArray {
+    ZoneArray(int level,GeodesicGrid &grid)
+      : level(level),nr_of_zones(GeodesicGrid::nrOfZones(level)),grid(grid)
+        {scale = 0.0;nr_of_stars = 0;}
+    virtual ~ZoneArray(void) {}
+    virtual ZoneData *getZone(int index) const = 0;
+    virtual void writeHeaderToOutput(FILE *f) const = 0;
+    const int level;
+    const int nr_of_zones;
+    const GeodesicGrid &grid;
+    double scale;
+//    int scale_int;
+    int nr_of_stars;
+    HipStar *addStar(int tyc1,int tyc2,int tyc3,
+                     int hip,const char *component_ids,
+                     double ra, // degrees
+                     double dec, // degrees
+                     double pma,double pmd,
+                     double mag,double b_v,
+                     double plx,const char *sp,
+                     bool &does_not_fit);
+    void writeOutput(const char *fname);
+  };
+  struct HipZoneArray : public ZoneArray {
+    HipZoneArray(int level,GeodesicGrid &grid)
+      : ZoneArray(level,grid),zones(new HipZoneData[nr_of_zones]) {}
+    ~HipZoneArray(void) {delete[] zones;}
+    void writeHeaderToOutput(FILE *f) const;
+    HipZoneData *const zones;
+    ZoneData *getZone(int index) const {
+      if (index<0 || index>=nr_of_zones) {
+        cerr << "getZone: bad index" << endl;
+        exit(-1);
+      }
+      return zones+index;
+    }
+  };
+  struct TycZoneArray : public ZoneArray {
+    TycZoneArray(int level,GeodesicGrid &grid)
+      : ZoneArray(level,grid),zones(new TycZoneData[nr_of_zones]) {}
+    ~TycZoneArray(void) {delete[] zones;}
+    void writeHeaderToOutput(FILE *f) const;
+    TycZoneData *zones;
+    ZoneData *getZone(int index) const {
+      if (index<0 || index>=nr_of_zones) {
+        cerr << "getZone: bad index" << endl;
+        exit(-1);
+      }
+      return zones+index;
+    }
+  };
+  struct FaintZoneArray : public ZoneArray {
+    FaintZoneArray(int level,GeodesicGrid &grid)
+      : ZoneArray(level,grid),zones(new FaintZoneData[nr_of_zones]) {}
+    ~FaintZoneArray(void) {delete[] zones;}
+    void writeHeaderToOutput(FILE *f) const;
+    FaintZoneData *zones;
+    ZoneData *getZone(int index) const {
+      if (index<0 || index>=nr_of_zones) {
+        cerr << "getZone: bad index" << endl;
+        exit(-1);
+      }
+      return zones+index;
+    }
+  };
+  ZoneArray *zone_array[MAX_LEVEL+1];
+  
+private:
+  static void initTriangleFunc(int lev,int index,
+                    const Vector &c0,
+                    const Vector &c1,
+                    const Vector &c2,
+                    void *user_data) {
+    reinterpret_cast<Accumulator*>(user_data)
+      ->initTriangle(lev,index,c0,c1,c2);
+  }
+  void initTriangle(int lev,int index,
+                    const Vector &c0,
+                    const Vector &c1,
+                    const Vector &c2);
+};
+
+
+
+
+
+
+
+
+
+static
+void WriteByte(FILE *f,int x) {
+  unsigned char c = (unsigned int)x;
+  if (1!=fwrite(&c,1,1,f)) {
+    cerr << "WriteByte: fwrite failed" << endl;
+    exit(-1);
+  }
+}
+
+static
+void WriteInt(FILE *f,int x) {
+  if (4!=fwrite(&x,1,4,f)) {
+    cerr << "WriteInt: fwrite failed" << endl;
+    exit(-1);
+  }
+}
+
+//static
+//int DoubleToInt(double x) {
+//  return (int)floor(0.5+x*0x7FFFFFFF);
+//}
+
+//static
+//double IntToDouble(int x) {
+//  double rval = x;
+//  rval /= 0x7FFFFFFF;
+//  return rval;
+//}
+
+void ZoneData::writeInfoToOutput(FILE *f) const {
+//  WriteInt(f,DoubleToInt(center.x[0]));
+//  WriteInt(f,DoubleToInt(center.x[1]));
+//  WriteInt(f,DoubleToInt(center.x[2]));
+  WriteInt(f,getArraySize());
+}
+
+void HipZoneData::writeStarsToOutput(FILE *f) {
+  stars.sort();
+  for (list<HipStar>::const_iterator it(stars.begin());
+       it!=stars.end();it++) {
+    if (sizeof(HipStarCompressed)!=
+        fwrite(&(*it),1,sizeof(HipStarCompressed),f)) {
+      cerr << "HipZoneData::writeStarsToOutput: fwrite failed" << endl;
+      exit(-1);
+    }
+  }
+}
+
+void TycZoneData::writeStarsToOutput(FILE *f) {
+  stars.sort();
+  for (list<TycStar>::const_iterator it(stars.begin());
+       it!=stars.end();it++) {
+    if (sizeof(TycStar)!=fwrite(&(*it),1,sizeof(TycStar),f)) {
+      cerr << "TycZoneData::writeStarsToOutput: fwrite failed" << endl;
+      exit(-1);
+    }
+  }
+}
+
+
+void FaintZoneData::writeStarsToOutput(FILE *f) {
+  stars.sort();
+  for (list<FaintStar>::const_iterator it(stars.begin());
+       it!=stars.end();it++) {
+    if (sizeof(FaintStar)!=fwrite(&(*it),1,sizeof(FaintStar),f)) {
+      cerr << "FaintZoneData::writeStarsToOutput: fwrite failed" << endl;
+      exit(-1);
+    }
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Accumulator::Accumulator(void)
+            :grid(MAX_LEVEL) {
+  int l=0;
+  for (;l<=MAX_HIP_LEVEL;l++) {
+    zone_array[l] = new HipZoneArray(l,grid);
+  }
+  for (;l<=MAX_TYC_LEVEL;l++) {
+    zone_array[l] = new TycZoneArray(l,grid);
+  }
+  for (;l<=MAX_LEVEL;l++) {
+    zone_array[l] = new FaintZoneArray(l,grid);
+  }
+  grid.visitTriangles(MAX_LEVEL,initTriangleFunc,this);
+  for (int l=0;l<=MAX_LEVEL;l++) {
+//    zone_array[l]->scale_int = DoubleToInt(zone_array[l]->scale);
+//    while (IntToDouble(zone_array[l]->scale_int) > zone_array[l]->scale) {
+//      zone_array[l]->scale_int--;
+//    }
+//    zone_array[l]->scale = IntToDouble(zone_array[l]->scale_int);
+    cout << "Accumulator::Accumulator: level " << zone_array[l]->level
+         << ", scale: " << zone_array[l]->scale << endl;
+  }
+}
+
+Accumulator::~Accumulator(void) {
+  for (int l=0;l<=MAX_LEVEL;l++) {
+    delete zone_array[l];
+  }
+}
+
+static const Vector north = {0.0,0.0,1.0};
+
+void Accumulator::initTriangle(int lev,int index,
+                               const Vector &c0,
+                               const Vector &c1,
+                               const Vector &c2) {
+  ZoneData &z(*zone_array[lev]->getZone(index));
+  double &scale(zone_array[lev]->scale);
+  z.center = c0+c1+c2;
+  z.center.normalize();
+  z.axis0 = north ^ z.center;
+  z.axis0.normalize();
+  z.axis1 = z.center ^ z.axis0;
+
+  double mu0,mu1,f,h;
+
+  mu0 = (c0-z.center)*z.axis0;
+  mu1 = (c0-z.center)*z.axis1;
+  f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
+  h = fabs(mu0)*f;
+  if (scale < h) scale = h;
+  h = fabs(mu1)*f;
+  if (scale < h) scale = h;
+
+  mu0 = (c1-z.center)*z.axis0;
+  mu1 = (c1-z.center)*z.axis1;
+  f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
+  h = fabs(mu0)*f;
+  if (scale < h) scale = h;
+  h = fabs(mu1)*f;
+  if (scale < h) scale = h;
+
+  mu0 = (c2-z.center)*z.axis0;
+  mu1 = (c2-z.center)*z.axis1;
+  f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
+  h = fabs(mu0)*f;
+  if (scale < h) scale = h;
+  h = fabs(mu1)*f;
+  if (scale < h) scale = h;
+}
+
+static
+unsigned char BvToColorIndex(double b_v) {
+  b_v *= 1000.0;
+  if (b_v < -500) {
+    b_v = -500;
+  } else if (b_v > 3499) {
+    b_v = 3499;
+  }
+  return (unsigned int)floor(0.5+127.0*((500.0+b_v)/4000.0));
+}
+
+HipStar *FaintZoneData::add(int level,
+                       int tyc1,int tyc2,int tyc3,
+                       int hip,const char *component_ids,
+                       double x0,double x1,double dx0,double dx1,
+                       double mag,double b_v,double plx,const char *sp,
+                       bool &does_not_fit) {
+  if (mag>=level_mag_limit[level]) {
+    cout << "too faint" << endl;
+    exit(-1);
+  }
+  FaintStar s;
+  s.x0  = ((int)floor(0.5+x0*((1<<17)-1))); // 18 bit signed
+  s.b_v = BvToColorIndex(b_v);       // 0..127: 7 Bit unsigned
+  s.x1  = ((int)floor(0.5+x1*((1<<17)-1))); // 18 bit signed
+  s.mag = (unsigned int)floor(30*(mag-level_mag_limit[level-1])/
+                        (level_mag_limit[level]-level_mag_limit[level-1]));
+                 // steps(accuracy) of 0.05mag, 5 bit unsigned
+  stars.push_back(s);
+  does_not_fit = false;
+  return 0;
+}
+
+HipStar *TycZoneData::add(int level,
+                          int tyc1,int tyc2,int tyc3,
+                          int hip,const char *component_ids,
+                          double x0,double x1,double dx0,double dx1,
+                          double mag,double b_v,double plx,const char *sp,
+                          bool &does_not_fit) {
+  TycStar s;
+  s.x0  = ((int)floor(0.5+x0*((1<<19)-1))); // 20 bit signed
+  s.b_v = BvToColorIndex(b_v);       // 0..127: 7 Bit unsigned
+  s.x1  = ((int)floor(0.5+x1*((1<<19)-1))); // 20 bit signed
+  s.mag = (unsigned int)floor(30*(mag-level_mag_limit[level-1])/
+                        (level_mag_limit[level]-level_mag_limit[level-1]));
+                 // steps(accuracy) of 0.05mag, 5 bit unsigned
+  const int dx0_int = (int)floor(0.5+dx0*10.0); // 14 bit signed
+  const int dx1_int = (int)floor(0.5+dx1*10.0); // 14 bit signed
+  if (dx0_int >= 8192 || dx0_int<-8192 ||
+      dx1_int >= 8192 || dx1_int<-8192) {
+      // does not fit, must store in Hip format
+    does_not_fit = true;
+    return 0;
+  }
+  s.dx0 = dx0_int;
+  s.dx1 = dx1_int;
+  stars.push_back(s);
+  does_not_fit = false;
+  return 0;
+}
+
+
+HipStar *HipZoneData::add(int level,
+                          int tyc1,int tyc2,int tyc3,
+                          int hip,const char *component_ids,
+                          double x0,double x1,double dx0,double dx1,
+                          double mag,double b_v,double plx,const char *sp,
+                          bool &does_not_fit) {
+  HipStar s;
+  s.x0  = ((int)floor(0.5+x0*((1<<31)-1))); // 32 bit signed
+  s.b_v = BvToColorIndex(b_v);       // 0..127: 7 Bit unsigned
+  s.x1  = ((int)floor(0.5+x1*((1<<31)-1))); // 32 bit signed
+  const int mag_int = (int)floor(20*
+             (mag-((level==0)?-2.0:level_mag_limit[level-1])));
+  if (mag_int < 0 || mag_int > 255) {
+    cerr << "HipZoneData(" << level << ")::add: bad mag: " << mag << endl;
+    exit(-1);
+  }
+  s.mag = mag_int;
+  s.dx0 = (int)floor(0.5+dx0*10.0);
+  s.dx1 = (int)floor(0.5+dx1*10.0);
+  s.hip = hip;
+  s.component_ids = component_ids;
+  s.comp_ids_int = 0;
+  s.sp_int = 0;
+  s.plx = (int)floor(0.5+100.0*plx);
+  s.sp = sp;
+  stars.push_back(s);
+  does_not_fit = false;
+  return &(stars.back());
+}
+
+
+HipStar *Accumulator::ZoneArray::addStar(
+                           int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
+                           double ra,double dec,double pma,double pmd,
+                           double mag,double b_v,double plx,const char *sp,
+                           bool &does_not_fit) {
+  ra *= (M_PI/180.0);
+  dec *= (M_PI/180.0);
+  if (ra < 0 || ra >= 2*M_PI || dec < -0.5*M_PI || dec > 0.5*M_PI) {
+    cerr << "ZoneArray(l=" << level << ")::addStar: "
+            "bad ra/dec: " << ra << ',' << dec << endl;
+    exit (-1);
+  }
+  const double cdec = cos(dec);
+  Vector pos;
+  pos.x[0] = cos(ra)*cdec;
+  pos.x[1] = sin(ra)*cdec;
+  pos.x[2] = sin(dec);
+  const int zone = grid.searchZone(pos,level);
+
+  ZoneData &z(*getZone(zone));
+  const double mu0 = (pos-z.center) * z.axis0;
+  const double mu1 = (pos-z.center) * z.axis1;
+  const double d = 1.0 - mu0*mu0 - mu1*mu1;
+  const double sd = sqrt(d);
+  const double x0 = mu0/(sd*scale);
+  const double x1 = mu1/(sd*scale);
+  const double h = mu0*pma + mu1*pmd;
+  const double dx0 = (pma*d + mu0*h) / (sd*d);
+  const double dx1 = (pmd*d + mu1*h) / (sd*d);
+
+//cout << "Accumulator::ZoneArray(l=" << level << ")::addStar: " << zone << endl;
+  HipStar *rval = z.add(level,tyc1,tyc2,tyc3,hip,component_ids,
+                        x0,x1,dx0,dx1,mag,b_v,plx,sp,does_not_fit);
+  if (!does_not_fit) nr_of_stars++;
+//cout << "Accumulator::ZoneArray(l=" << level << ")::addStar: 999";
+  return rval;
+}
+
+int Accumulator::addStar(int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
+                         double ra,double dec,double pma,double pmd,
+                         double mag,double b_v,double plx,const char *sp) {
+//  const int packed_hip = PackHip(hip,component_ids);
+//  const int packed_tyc = PackTyc(tyc1,tyc2,tyc3);
+
+  int l=0;
+  while (l<MAX_LEVEL && mag>=level_mag_limit[l]) l++;
+  if (hip>0 && l>MAX_HIP_LEVEL) l = MAX_HIP_LEVEL;
+//  store the faint tyc2 stars as faint stars
+//  else if (tyc1>0 && l>MAX_TYC_LEVEL) l = MAX_TYC_LEVEL;
+
+
+  if (restrict_output_level_max < l ||
+      restrict_output_level_min > l) {
+      // you may want to restrict the output
+      // in order to restrict the amount of main memory used.
+      // Calling the program twice may be faster than calling once
+      // without enough main memory (avoid swapping).
+    return 0;
+  }
+
+  bool does_not_fit;
+  HipStar *s = zone_array[l]->addStar(tyc1,tyc2,tyc3,hip,component_ids,
+                                      ra,dec,pma,pmd,
+                                      mag,b_v,plx,sp,does_not_fit);
+  if (does_not_fit) s = zone_array[MAX_HIP_LEVEL]
+                            ->addStar(tyc1,tyc2,tyc3,hip,component_ids,
+                                      ra,dec,pma,pmd,
+                                      mag,b_v,plx,sp,does_not_fit);
+  if (s) hip_index[hip].push_back(s);
+  return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#define FILE_MAGIC 0x835f040a
+
+void Accumulator::HipZoneArray::writeHeaderToOutput(FILE *f) const {
+  cout << "Accumulator::HipZoneArray(" << level << ")::writeHeaderToOutput: "
+       << nr_of_stars << endl;
+  WriteInt(f,FILE_MAGIC);
+  WriteInt(f,0); // type
+  WriteInt(f,0); // major version
+  WriteInt(f,4); // minor version
+  WriteInt(f,level);
+//  WriteInt(f,scale_int);
+  if (level == 0) {
+    WriteInt(f,-2000); // min_mag
+  } else {
+  WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
+  }
+  WriteInt(f,12800); // mag_range
+  WriteInt(f,256);   // mag_steps
+}
+
+void Accumulator::TycZoneArray::writeHeaderToOutput(FILE *f) const {
+  cout << "Accumulator::TycZoneArray(" << level << ")::writeHeaderToOutput: "
+       << nr_of_stars << endl;
+  WriteInt(f,FILE_MAGIC);
+  WriteInt(f,1); // type
+  WriteInt(f,0); // major version
+  WriteInt(f,3); // minor version
+  WriteInt(f,level);
+//  WriteInt(f,scale_int);
+  WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
+  WriteInt(f,(int)floor(0.5+1000.0*(level_mag_limit[level]
+                                   -level_mag_limit[level-1]))); // mag_range
+  WriteInt(f,30);   // mag_steps
+}
+
+void Accumulator::FaintZoneArray::writeHeaderToOutput(FILE *f) const {
+  cout << "Accumulator::FaintZoneArray(" << level << ")::writeHeaderToOutput: "
+       << nr_of_stars << endl;
+  WriteInt(f,FILE_MAGIC);
+  WriteInt(f,2); // type
+  WriteInt(f,0); // major version
+  WriteInt(f,1); // minor version
+  WriteInt(f,level);
+//  WriteInt(f,scale_int);
+  WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
+  WriteInt(f,(int)floor(0.5+1000.0*(level_mag_limit[level]
+                                   -level_mag_limit[level-1]))); // mag_range
+  WriteInt(f,30);   // mag_steps
+}
+
+
+void Accumulator::ZoneArray::writeOutput(const char *fname) {
+  if (nr_of_stars <= 0) return;
+  FILE *f = fopen(fname,"wb");
+  if (f==0) {
+    cerr << "Accumulator::writeOutput(" << fname << "): "
+            "fopen failed" << endl;
+    return;
+  }
+  writeHeaderToOutput(f);
+  for (int zone=0;zone<nr_of_zones;zone++) {
+    getZone(zone)->writeInfoToOutput(f);
+  }
+  int permille = 0;
+  for (int zone=0;zone<nr_of_zones;zone++) {
+    getZone(zone)->writeStarsToOutput(f);
+    int p = (1000*zone)/nr_of_zones;
+    if (p != permille) {
+      permille = p;
+      cout << "\rAccumulator::writeOutput(" << fname << "): "
+           << permille << "permille written" << flush;
+    }
+  }
+  cout << endl;
+  fclose(f);
+}
+
+void Accumulator::writeOutput(const char *fnames[]) {
+  for (int l=0;l<=MAX_LEVEL;l++) {
+    zone_array[l]->writeOutput(fnames[l]);
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+int ReadHipTycFile(Accumulator &accu) {
+  int count = 0;
+  FILE *f;
+  const char *fname = "HipTyc";
+  f = fopen(fname,"r");
+  if (f == 0) {
+    fprintf(stderr,"Could not open file \"%s\".\n",fname);
+    exit(-1);
+  }
+
+  int hip,tyc1,tyc2,tyc3;
+  char cids[32];
+  char sp[256];
+  int mag,b_v,VarFlag;
+  double ra,dec,Plx,pm_ra,pm_dec;
+  while (14==fscanf(f,"%d%d%d%d%s%d%lf%lf%lf%lf%lf%d%d%s",
+                    &hip,&tyc1,&tyc2,&tyc3,cids,&VarFlag,
+                    &ra,&dec,&Plx,&pm_ra,&pm_dec,&mag,&b_v,sp)) {
+      if (b_v>-500 && b_v<3450) {
+      const int rc = accu.addStar(tyc1,tyc2,tyc3,hip,
+                                  cids[0]=='?'?"":cids,
+                                  ra, // degrees
+                                  dec, // degrees
+                                  pm_ra,pm_dec,0.001*mag,0.001*b_v,
+                                  Plx,sp[0]=='?'?"":sp);
+      if (rc < 0) {
+          // never mind: propably no magnitude for Hiparcos star
+//        fprintf(stderr,"File \"%s\", record %d: Error 13 %d %d \"%s\"\n",
+//                fname,count,rc,hip,sp);
+//        exit(-1);
+      }
+    count++;
+    }
+  }
+  fclose(f);
+  return count;
+}
+
+
+const unsigned short int SHORT_ASTSRCBIT0 = 0x0001; // Astrometry source bit 0
+const unsigned short int SHORT_ASTSRCBIT1 = 0x0002; // Astrometry source bit 1
+const unsigned short int SHORT_ASTSRCBIT2 = 0x0004; // Astrometry source bit 2
+const unsigned short int SHORT_UBBIT   = 0x0008;
+const unsigned short int SHORT_TMBIT   = 0x0010;
+const unsigned short int SHORT_XRBIT   = 0x0020;
+const unsigned short int SHORT_IUCBIT  = 0x0040;
+const unsigned short int SHORT_ITYBIT  = 0x0080;
+const unsigned short int SHORT_OMAGBIT = 0x0100;
+const unsigned short int SHORT_EMAGBIT = 0x0200;
+const unsigned short int SHORT_TMONLY  = 0x0400;
+const unsigned short int SHORT_SPIKE   = 0x0800;
+const unsigned short int SHORT_TYCONF  = 0x1000;
+const unsigned short int SHORT_BSCONF  = 0x2000;
+const unsigned short int SHORT_BSART   = 0x4000;
+const unsigned short int SHORT_USEME   = 0x8000;
+
+struct Short_NOMAD_Record {
+  int ra,spd,pm_ra,pm_spd;
+  short int b,v,r;
+  unsigned short int flags;
+};
+
+
+
+#define READ_SIZE 10000
+static Short_NOMAD_Record buff[READ_SIZE];
+
+void ReadNOMADFile(const char *fname,Accumulator &accu) {
+  int count;
+  FILE *f;
+  f = fopen(fname,"r");
+  if (f == 0) {
+    fprintf(stderr,"Could not open file \"%s\".\n",fname);
+    exit(-1);
+  }
+  int total = 0;
+  do {
+    count = fread(buff,sizeof(Short_NOMAD_Record),READ_SIZE,f);
+    total += count;
+    printf("\rfread(%s,...) returned %6d, total = %8d",
+           fname,count,total);
+    fflush(stdout);
+    int i;
+    for (i=0;i<count;i++) {
+      int mag = 30000;
+      int b_v;
+      if (buff[i].v >= 30000) {
+        if (buff[i].b >= 30000) {
+          if (buff[i].r >= 30000) {
+//          cerr << "no magnitude at all" << endl;
+          } else {
+            mag = buff[i].r + 3499; // just an assumption
+            b_v = 3499;
+          }
+        } else {
+          mag = buff[i].b + 500; // just an assumption
+          b_v = -500;
+        }
+      } else {
+        mag = buff[i].v;
+        if (buff[i].b >= 30000) {
+          if (buff[i].r >= 30000) {
+            b_v = 0;
+          } else {
+            b_v = buff[i].v-buff[i].r; // desperate
+          }
+        } else {
+          b_v = buff[i].b-buff[i].v;
+        }
+      }
+
+      int nr_of_measurements = 0;
+      if (buff[i].b < 30000) nr_of_measurements++;
+      if (buff[i].v < 30000) nr_of_measurements++;
+      if (buff[i].r < 30000) nr_of_measurements++;
+
+      if (mag < 19500 && b_v>-500 && b_v<3450 &&
+          ((buff[i].flags&SHORT_USEME) ||
+           (
+           ((buff[i].flags&(SHORT_ASTSRCBIT0|SHORT_ASTSRCBIT1|SHORT_ASTSRCBIT2))!=1 ||
+              (((buff[i].flags&SHORT_UBBIT)==0)
+               &&
+                  (mag>14000 || nr_of_measurements>1) &&
+                  (mag>13000 || nr_of_measurements>2)
+                   )) &&
+            ((buff[i].flags&SHORT_SPIKE)==0) &&
+            ((buff[i].flags&SHORT_BSART)==0) &&
+            ((buff[i].flags&SHORT_TMONLY)==0)))) {
+        if (accu.addStar(0,0,0,0,"",
+                         buff[i].ra/(3600.0*1000.0), // degrees
+                         (buff[i].spd-90*3600*1000)/(3600.0*1000.0), // degrees
+                         0.1*buff[i].pm_ra,0.1*buff[i].pm_spd,
+                         0.001*mag,0.001*b_v,
+                         0,"") < 0) {
+          fprintf(stderr,"File \"%s\", record %d: Error 16\n",fname,count);
+          exit(-1);
+        }
+      }
+    }
+  } while (count == READ_SIZE);
+  printf("\n");
+  fclose(f);
+}
+
+
+int main(int argc,char *argv[]) {
+  if (argc != 3 ||
+      1 != sscanf(argv[1],"%d",&restrict_output_level_min) ||
+      1 != sscanf(argv[2],"%d",&restrict_output_level_max)) {
+    cerr << "Usage: " << argv[0] << " level_min level_max" << endl
+         << " (like " << argv[0] << " 0 6)" << endl << endl;
+    return -1;
+  }
+
+  Accumulator accu;
+  int n=0;
+  n = ReadHipTycFile(accu);
+  printf("HipTyc: %d records processed.\n",n);
+  SqueezeHip();
+  
+  for (int c=0;nomad_names[c];c++) {
+    ReadNOMADFile(nomad_names[c],accu);
+  }
+
+  accu.writeOutput(output_file_names);
+
+  return 0;
+}