1 // Author and Copyright: Johannes Gajdosik, 2006
2 // Modifications: Alexander Wolf, 2013
3 //
4 // This program is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU General Public License
6 // as published by the Free Software Foundation; either version 2
7 // of the License, or (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program; if not, write to the Free Software
16 // Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
17 //
18 // g++ -O2 MakeCombinedCatalogue.C -o MakeCombinedCatalogue
19
20
21 // change catalogue locations according to your needs:
22
23 const char *nomad_names[]={
24 "/storage/1T/nomad.sml/Nomad00.sml",
25 "/storage/1T/nomad.sml/Nomad01.sml",
26 "/storage/1T/nomad.sml/Nomad02.sml",
27 "/storage/1T/nomad.sml/Nomad03.sml",
28 "/storage/1T/nomad.sml/Nomad04.sml",
29 "/storage/1T/nomad.sml/Nomad05.sml",
30 "/storage/1T/nomad.sml/Nomad06.sml",
31 "/storage/1T/nomad.sml/Nomad07.sml",
32 "/storage/1T/nomad.sml/Nomad08.sml",
33 "/storage/1T/nomad.sml/Nomad09.sml",
34 "/storage/1T/nomad.sml/Nomad10.sml",
35 "/storage/1T/nomad.sml/Nomad11.sml",
36 "/storage/1T/nomad.sml/Nomad12.sml",
37 "/storage/1T/nomad.sml/Nomad13.sml",
38 "/storage/1T/nomad.sml/Nomad14.sml",
39 "/storage/1T/nomad.sml/Nomad15.sml",
40 "/storage/1T/nomad.sml/Nomad16.sml",
41 0
42 };
43
44
45
46 #define NR_OF_HIP 120416
47
48 #define MAX_HIP_LEVEL 2
49 #define MAX_TYC_LEVEL 4
50 #define MAX_LEVEL 9
51
52 static
53 const char *output_file_names[MAX_LEVEL+1] = {
54 "stars0.cat",
55 "stars1.cat",
56 "stars2.cat",
57 "stars3.cat",
58 "stars4.cat",
59 "stars5.cat",
60 "stars6.cat",
61 "stars7.cat",
62 "stars8.cat",
63 "stars9.cat"
64 };
65
66 static
67 const char *component_ids_filename = "stars_hip_component_ids.cat";
68
69 static
70 const char *sp_filename = "stars_hip_sp.cat";
71
72 static
73 const double level_mag_limit[MAX_LEVEL+1] = {
74 6,7.5,9,
75 10.5,12,
76 13.5,15,16.5,18,19.5
77 };
78
79
80
81
82
83
84
85
86
87 /*
88
89 Vector (simple 3d-Vector)
90
91 Author and Copyright: Johannes Gajdosik, 2006
92
93 License: GPL
94
95 */
96
97
98 #include <math.h>
99
100 struct Vector {
101 double x[3];
102 inline double length(void) const;
103 inline void normalize(void);
104 const Vector &operator+=(const Vector &a) {
105 x[0] += a.x[0];
106 x[1] += a.x[1];
107 x[2] += a.x[2];
108 }
109 const Vector &operator-=(const Vector &a) {
110 x[0] -= a.x[0];
111 x[1] -= a.x[1];
112 x[2] -= a.x[2];
113 }
114 };
115
116 static inline
117 Vector operator^(const Vector &a,const Vector &b) {
118 Vector c;
119 c.x[0] = a.x[1]*b.x[2] - a.x[2]*b.x[1];
120 c.x[1] = a.x[2]*b.x[0] - a.x[0]*b.x[2];
121 c.x[2] = a.x[0]*b.x[1] - a.x[1]*b.x[0];
122 return c;
123 }
124
125 static inline
126 Vector operator-(const Vector &a,const Vector &b) {
127 Vector c;
128 c.x[0] = a.x[0]-b.x[0];
129 c.x[1] = a.x[1]-b.x[1];
130 c.x[2] = a.x[2]-b.x[2];
131 return c;
132 }
133
134 static inline
135 Vector operator+(const Vector &a,const Vector &b) {
136 Vector c;
137 c.x[0] = a.x[0]+b.x[0];
138 c.x[1] = a.x[1]+b.x[1];
139 c.x[2] = a.x[2]+b.x[2];
140 return c;
141 }
142
143 static inline
144 Vector operator*(double a,const Vector &b) {
145 Vector c;
146 c.x[0] = a*b.x[0];
147 c.x[1] = a*b.x[1];
148 c.x[2] = a*b.x[2];
149 return c;
150 }
151
152 static inline
153 Vector operator*(const Vector &b,double a) {
154 Vector c;
155 c.x[0] = a*b.x[0];
156 c.x[1] = a*b.x[1];
157 c.x[2] = a*b.x[2];
158 return c;
159 }
160
161 static inline
162 double operator*(const Vector &a,const Vector &b) {
163 return a.x[0]*b.x[0]+a.x[1]*b.x[1]+a.x[2]*b.x[2];
164 }
165
length(void)166 double Vector::length(void) const {
167 return sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
168 }
169
normalize(void)170 void Vector::normalize(void) {
171 const double l = 1.0/length();
172 x[0] *= l;
173 x[1] *= l;
174 x[2] *= l;
175 }
176
177
178
179
ChangeEpoch(double delta_years,double & ra,double & dec,double pm_ra,double pm_dec)180 void ChangeEpoch(double delta_years,
181 double &ra,double &dec, // degrees
182 double pm_ra,double pm_dec // mas/yr
183 ) {
184 ra *= (M_PI/180);
185 dec *= (M_PI/180);
186
187 const double cdec = cos(dec);
188 Vector x = {cos(ra)*cdec,sin(ra)*cdec,sin(dec)};
189 const Vector north = {0.0,0.0,1.0};
190
191 Vector axis0 = north ^ x;
192 axis0.normalize();
193 const Vector axis1 = x ^ axis0;
194
195 const double f = delta_years*(0.001/3600)*(M_PI/180);
196
197 x += pm_ra*f*axis1 + pm_dec*f*axis0;
198
199 ra = atan2(x.x[1],x.x[0]);
200 if (ra < 0.0) ra += 2*M_PI;
201 dec = atan2(x.x[2],sqrt(x.x[0]*x.x[0]+x.x[1]*x.x[1]));
202
203 ra *= (180/M_PI);
204 dec *= (180/M_PI);
205 }
206
207
208
209
210
211
212 //------------------------------------------------
213
214 /*
215
216 GeodesicGrid: a library for dividing the sphere into triangle zones
217 by subdividing the icosahedron
218
219 Author and Copyright: Johannes Gajdosik, 2006
220
221 This library requires a simple Vector library,
222 which may have different copyright and license.
223
224 In the moment I choose to distribute the library under the GPL,
225 later I may choose to additionally distribute it under a more
226 relaxed license like the LGPL. If you want to have the library
227 under another license, please ask me.
228
229
230
231 This library is free software; you can redistribute it and/or
232 modify it under the terms of the GNU General Public License
233 as published by the Free Software Foundation; either version 2
234 of the License, or (at your option) any later version.
235
236 This library is distributed in the hope that it will be useful,
237 but WITHOUT ANY WARRANTY; without even the implied warranty of
238 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
239 GNU General Public License for more details.
240
241 You should have received a copy of the GNU General Public License
242 along with this library; if not, write to the Free Software
243 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
244
245 */
246
247 using namespace std;
248
249 struct HalfSpace {
250 // all x with x*n-d >= 0
251 Vector n;
252 double d;
insideHalfSpace253 bool inside(const Vector &x) const {return (x*n>=d);}
254 };
255
256 class GeodesicGrid {
257 // Grid of triangles (zones) on the sphere with radius 1,
258 // generated by subdividing the icosahedron.
259 // level 0: just the icosahedron, 20 zones
260 // level 1: 80 zones, level 2: 320 zones, ...
261 // Each zone has a unique integer number in the range
262 // [0,getNrOfZones()-1].
263 public:
264 GeodesicGrid(int max_level);
265 ~GeodesicGrid(void);
getMaxLevel(void)266 int getMaxLevel(void) const {return max_level;}
nrOfZones(int level)267 static int nrOfZones(int level)
268 {return (20<<(level<<1));} // 20*4^level
getNrOfZones(void)269 int getNrOfZones(void) const {return nrOfZones(max_level);}
270 void getTriangleCorners(int lev, int index, Vector& c0, Vector& c1, Vector& c2) const;
271 // Return the position of the 3 corners for the triangle at the given level and index
272 int getPartnerTriangle(int lev, int index) const;
273 // Return the index of the partner triangle with which to form a paralelogram
274 typedef void (VisitFunc)(int lev,int index,
275 const Vector &c0,
276 const Vector &c1,
277 const Vector &c2,
278 void *context);
279 void visitTriangles(int max_visit_level,
280 VisitFunc *func,void *context) const;
281
282 int searchZone(const Vector &v,int search_level) const;
283 // find the zone number in which a given point lies.
284 // prerequisite: v*v==1
285 // When the point lies on the border of two or more zones,
286 // one such zone is returned (always the same one,
287 // because the algorithm is deterministic).
288
289 void searchZones(const HalfSpace *half_spaces,int nr_of_half_spaces,
290 int **inside,int **border,int max_search_level) const;
291 // find all zones that lie fully(inside) or partly(border)
292 // in the intersection of the given half spaces.
293 // The result is accurate when (0,0,0) lies on the border of
294 // each half space. If this is not the case,
295 // the result may be inaccurate, because it is assumed, that
296 // a zone lies in a half space when its 3 corners lie in
297 // this half space.
298 // inside[l] points to the begin of an integer array of size nrOfZones(l),
299 // border[l] points to one after the end of the same integer array.
300 // The array will be filled from the beginning with the inside zone numbers
301 // and from the end with the border zone numbers of the given level l
302 // for 0<=l<=getMaxLevel().
303 // inside[l] will not contain zones that are already contained
304 // in inside[l1] for some l1 < l.
305 // In order to restrict search depth set max_search_level < max_level,
306 // for full search depth set max_search_level = max_level,
307
308 private:
309 const Vector& getTriangleCorner(int lev, int index, int cornerNumber) const;
310 void initTriangle(int lev,int index,
311 const Vector &c0,
312 const Vector &c1,
313 const Vector &c2);
314 void visitTriangles(int lev,int index,
315 const Vector &c0,
316 const Vector &c1,
317 const Vector &c2,
318 int max_visit_level,
319 VisitFunc *func,
320 void *context) const;
321 void searchZones(int lev,int index,
322 const HalfSpace *half_spaces,
323 const int nr_of_half_spaces,
324 const int *index_of_used_half_spaces,
325 const int half_spaces_used,
326 const bool *corner0_inside,
327 const bool *corner1_inside,
328 const bool *corner2_inside,
329 int **inside,int **border,int max_search_level) const;
330 private:
331 const int max_level;
332 struct Triangle {
333 Vector e0,e1,e2; // Seitenmittelpunkte
334 };
335 Triangle **triangles;
336 // 20*(4^0+4^1+...+4^n)=20*(4*(4^n)-1)/3 triangles total
337 // 2+10*4^n corners
338 };
339
340 class GeodesicSearchResult {
341 public:
342 GeodesicSearchResult(const GeodesicGrid &grid);
343 ~GeodesicSearchResult(void);
344 void search(const HalfSpace *half_spaces,
345 const int nr_of_half_spaces,
346 int max_search_level);
347 void search(const Vector &e0,const Vector &e1,const Vector &e2,const Vector &e3,
348 int max_search_level);
349 void print(void) const;
350 private:
351 friend class GeodesicSearchInsideIterator;
352 friend class GeodesicSearchBorderIterator;
353 const GeodesicGrid &grid;
354 int **const zones;
355 int **const inside;
356 int **const border;
357 };
358
359 class GeodesicSearchBorderIterator {
360 public:
GeodesicSearchBorderIterator(const GeodesicSearchResult & r,int level)361 GeodesicSearchBorderIterator(const GeodesicSearchResult &r,int level)
362 : r(r),level((level<0)?0:(level>r.grid.getMaxLevel())
363 ?r.grid.getMaxLevel():level),
364 end(r.zones[GeodesicSearchBorderIterator::level]+
365 GeodesicGrid::nrOfZones(GeodesicSearchBorderIterator::level))
366 {reset();}
reset(void)367 void reset(void) {index = r.border[level];}
next(void)368 int next(void) // returns -1 when finished
369 {if (index < end) {return *index++;} return -1;}
370 private:
371 const GeodesicSearchResult &r;
372 const int level;
373 const int *const end;
374 const int *index;
375 };
376
377
378 class GeodesicSearchInsideIterator {
379 public:
GeodesicSearchInsideIterator(const GeodesicSearchResult & r,int level)380 GeodesicSearchInsideIterator(const GeodesicSearchResult &r,int level)
381 : r(r),max_level((level<0)?0:(level>r.grid.getMaxLevel())
382 ?r.grid.getMaxLevel():level)
383 {reset();}
384 void reset(void);
385 int next(void); // returns -1 when finished
386 private:
387 const GeodesicSearchResult &r;
388 const int max_level;
389 int level;
390 int max_count;
391 int *index_p;
392 int *end_p;
393 int index;
394 int count;
395 };
396
397 #include <stdlib.h>
398 #include <cassert>
399 #include <iostream>
400
401 using namespace std;
402
403 static const double icosahedron_G = 0.5*(1.0+sqrt(5.0));
404 static const double icosahedron_b = 1.0/sqrt(1.0+icosahedron_G*icosahedron_G);
405 static const double icosahedron_a = icosahedron_b*icosahedron_G;
406
407 static const Vector icosahedron_corners[12] = {
408 { icosahedron_a, -icosahedron_b, 0.0},
409 { icosahedron_a, icosahedron_b, 0.0},
410 {-icosahedron_a, icosahedron_b, 0.0},
411 {-icosahedron_a, -icosahedron_b, 0.0},
412 { 0.0, icosahedron_a, -icosahedron_b},
413 { 0.0, icosahedron_a, icosahedron_b},
414 { 0.0, -icosahedron_a, icosahedron_b},
415 { 0.0, -icosahedron_a, -icosahedron_b},
416 {-icosahedron_b, 0.0, icosahedron_a},
417 { icosahedron_b, 0.0, icosahedron_a},
418 { icosahedron_b, 0.0, -icosahedron_a},
419 {-icosahedron_b, 0.0, -icosahedron_a}
420 };
421
422 struct TopLevelTriangle {
423 int corners[3]; // index der Ecken
424 };
425
426
427 static const TopLevelTriangle icosahedron_triangles[20] = {
428 { 1, 0,10}, // 1
429 { 0, 1, 9}, // 0
430 { 0, 9, 6}, // 12
431 { 9, 8, 6}, // 9
432 { 0, 7,10}, // 16
433 { 6, 7, 0}, // 6
434 { 7, 6, 3}, // 7
435 { 6, 8, 3}, // 14
436 {11,10, 7}, // 11
437 { 7, 3,11}, // 18
438 { 3, 2,11}, // 3
439 { 2, 3, 8}, // 2
440 {10,11, 4}, // 10
441 { 2, 4,11}, // 19
442 { 5, 4, 2}, // 5
443 { 2, 8, 5}, // 15
444 { 4, 1,10}, // 17
445 { 4, 5, 1}, // 4
446 { 5, 9, 1}, // 13
447 { 8, 9, 5} // 8
448 };
449
450
451
GeodesicGrid(const int lev)452 GeodesicGrid::GeodesicGrid(const int lev)
453 :max_level(lev<0?0:lev) {
454 if (max_level > 0) {
455 triangles = new Triangle*[max_level];
456 int nr_of_triangles = 20;
457 for (int i=0;i<max_level;i++) {
458 triangles[i] = new Triangle[nr_of_triangles];
459 nr_of_triangles *= 4;
460 }
461 for (int i=0;i<20;i++) {
462 const int *const corners = icosahedron_triangles[i].corners;
463 initTriangle(0,i,
464 icosahedron_corners[corners[0]],
465 icosahedron_corners[corners[1]],
466 icosahedron_corners[corners[2]]);
467 }
468 } else {
469 triangles = 0;
470 }
471 }
472
~GeodesicGrid(void)473 GeodesicGrid::~GeodesicGrid(void) {
474 if (max_level > 0) {
475 for (int i=max_level-1;i>=0;i--) delete[] triangles[i];
476 delete triangles;
477 }
478 }
479
getTriangleCorners(int lev,int index,Vector & h0,Vector & h1,Vector & h2)480 void GeodesicGrid::getTriangleCorners(int lev,int index,
481 Vector &h0,
482 Vector &h1,
483 Vector &h2) const {
484 if (lev <= 0) {
485 const int *const corners = icosahedron_triangles[index].corners;
486 h0 = icosahedron_corners[corners[0]];
487 h1 = icosahedron_corners[corners[1]];
488 h2 = icosahedron_corners[corners[2]];
489 } else {
490 lev--;
491 const int i = index>>2;
492 Triangle &t(triangles[lev][i]);
493 switch (index&3) {
494 case 0: {
495 Vector c0,c1,c2;
496 getTriangleCorners(lev,i,c0,c1,c2);
497 h0 = c0;
498 h1 = t.e2;
499 h2 = t.e1;
500 } break;
501 case 1: {
502 Vector c0,c1,c2;
503 getTriangleCorners(lev,i,c0,c1,c2);
504 h0 = t.e2;
505 h1 = c1;
506 h2 = t.e0;
507 } break;
508 case 2: {
509 Vector c0,c1,c2;
510 getTriangleCorners(lev,i,c0,c1,c2);
511 h0 = t.e1;
512 h1 = t.e0;
513 h2 = c2;
514 } break;
515 case 3:
516 h0 = t.e0;
517 h1 = t.e1;
518 h2 = t.e2;
519 break;
520 }
521 }
522 }
523
getPartnerTriangle(int lev,int index)524 int GeodesicGrid::getPartnerTriangle(int lev, int index) const
525 {
526 if (lev==0)
527 {
528 assert(index<20);
529 return (index&1) ? index+1 : index-1;
530 }
531 switch(index&7)
532 {
533 case 2:
534 case 6:
535 return index+1;
536 case 3:
537 case 7:
538 return index-1;
539 case 0:
540 return (lev==1) ? index+5 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
541 case 1:
542 return (lev==1) ? index+3 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
543 case 4:
544 return (lev==1) ? index-3 : (getPartnerTriangle(lev-1, index>>2)<<2)+1;
545 case 5:
546 return (lev==1) ? index-5 : (getPartnerTriangle(lev-1, index>>2)<<2)+0;
547 default:
548 assert(0);
549 }
550 }
551
initTriangle(int lev,int index,const Vector & c0,const Vector & c1,const Vector & c2)552 void GeodesicGrid::initTriangle(int lev,int index,
553 const Vector &c0,
554 const Vector &c1,
555 const Vector &c2) {
556 Triangle &t(triangles[lev][index]);
557 t.e0 = c1+c2;
558 t.e0.normalize();
559 t.e1 = c2+c0;
560 t.e1.normalize();
561 t.e2 = c0+c1;
562 t.e2.normalize();
563 lev++;
564 if (lev < max_level) {
565 index *= 4;
566 initTriangle(lev,index+0,c0,t.e2,t.e1);
567 initTriangle(lev,index+1,t.e2,c1,t.e0);
568 initTriangle(lev,index+2,t.e1,t.e0,c2);
569 initTriangle(lev,index+3,t.e0,t.e1,t.e2);
570 }
571 }
572
573
visitTriangles(int max_visit_level,VisitFunc * func,void * context)574 void GeodesicGrid::visitTriangles(int max_visit_level,
575 VisitFunc *func,
576 void *context) const {
577 if (func && max_visit_level >= 0) {
578 if (max_visit_level > max_level) max_visit_level = max_level;
579 for (int i=0;i<20;i++) {
580 const int *const corners = icosahedron_triangles[i].corners;
581 visitTriangles(0,i,
582 icosahedron_corners[corners[0]],
583 icosahedron_corners[corners[1]],
584 icosahedron_corners[corners[2]],
585 max_visit_level,func,context);
586 }
587 }
588 }
589
visitTriangles(int lev,int index,const Vector & c0,const Vector & c1,const Vector & c2,int max_visit_level,VisitFunc * func,void * context)590 void GeodesicGrid::visitTriangles(int lev,int index,
591 const Vector &c0,
592 const Vector &c1,
593 const Vector &c2,
594 int max_visit_level,
595 VisitFunc *func,
596 void *context) const {
597 (*func)(lev,index,c0,c1,c2,context);
598 Triangle &t(triangles[lev][index]);
599 lev++;
600 if (lev <= max_visit_level) {
601 index *= 4;
602 visitTriangles(lev,index+0,c0,t.e2,t.e1,max_visit_level,func,context);
603 visitTriangles(lev,index+1,t.e2,c1,t.e0,max_visit_level,func,context);
604 visitTriangles(lev,index+2,t.e1,t.e0,c2,max_visit_level,func,context);
605 visitTriangles(lev,index+3,t.e0,t.e1,t.e2,max_visit_level,func,context);
606 }
607 }
608
609
searchZone(const Vector & v,int search_level)610 int GeodesicGrid::searchZone(const Vector &v,int search_level) const {
611 for (int i=0;i<20;i++) {
612 const int *const corners = icosahedron_triangles[i].corners;
613 const Vector &c0(icosahedron_corners[corners[0]]);
614 const Vector &c1(icosahedron_corners[corners[1]]);
615 const Vector &c2(icosahedron_corners[corners[2]]);
616 if (((c0^c1)*v >= 0.0) && ((c1^c2)*v >= 0.0) && ((c2^c0)*v >= 0.0)) {
617 // v lies inside this icosahedron triangle
618 for (int lev=0;lev<search_level;lev++) {
619 Triangle &t(triangles[lev][i]);
620 i <<= 2;
621 if ((t.e1^t.e2)*v <= 0.0) {
622 // i += 0;
623 } else if ((t.e2^t.e0)*v <= 0.0) {
624 i += 1;
625 } else if ((t.e0^t.e1)*v <= 0.0) {
626 i += 2;
627 } else {
628 i += 3;
629 }
630 }
631 return i;
632 }
633 }
634 cerr << "software error: not found" << endl;
635 exit(-1);
636 // shut up the compiler
637 return -1;
638 }
639
searchZones(const HalfSpace * half_spaces,const int nr_of_half_spaces,int ** inside_list,int ** border_list,int max_search_level)640 void GeodesicGrid::searchZones(const HalfSpace *half_spaces,
641 const int nr_of_half_spaces,
642 int **inside_list,int **border_list,
643 int max_search_level) const {
644 if (max_search_level < 0) max_search_level = 0;
645 else if (max_search_level > max_level) max_search_level = max_level;
646 int halfs_used[nr_of_half_spaces];
647 for (int h=0;h<nr_of_half_spaces;h++) {halfs_used[h] = h;}
648 bool corner_inside[12][nr_of_half_spaces];
649 for (int h=0;h<nr_of_half_spaces;h++) {
650 const HalfSpace &half_space(half_spaces[h]);
651 for (int i=0;i<12;i++) {
652 corner_inside[i][h] = half_space.inside(icosahedron_corners[i]);
653 }
654 }
655 for (int i=0;i<20;i++) {
656 searchZones(0,i,
657 half_spaces,nr_of_half_spaces,halfs_used,nr_of_half_spaces,
658 corner_inside[icosahedron_triangles[i].corners[0]],
659 corner_inside[icosahedron_triangles[i].corners[1]],
660 corner_inside[icosahedron_triangles[i].corners[2]],
661 inside_list,border_list,max_search_level);
662 }
663 }
664
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)665 void GeodesicGrid::searchZones(int lev,int index,
666 const HalfSpace *half_spaces,
667 const int nr_of_half_spaces,
668 const int *index_of_used_half_spaces,
669 const int half_spaces_used,
670 const bool *corner0_inside,
671 const bool *corner1_inside,
672 const bool *corner2_inside,
673 int **inside_list,int **border_list,
674 const int max_search_level) const {
675 int halfs_used[half_spaces_used];
676 int halfs_used_count = 0;
677 for (int h=0;h<half_spaces_used;h++) {
678 const int i = index_of_used_half_spaces[h];
679 if (!corner0_inside[i] && !corner1_inside[i] && !corner2_inside[i]) {
680 // totally outside this HalfSpace
681 return;
682 } else if (corner0_inside[i] && corner1_inside[i] && corner2_inside[i]) {
683 // totally inside this HalfSpace
684 } else {
685 // on the border of this HalfSpace
686 halfs_used[halfs_used_count++] = i;
687 }
688 }
689 if (halfs_used_count == 0) {
690 // this triangle(lev,index) lies inside all halfspaces
691 **inside_list = index;
692 (*inside_list)++;
693 } else {
694 (*border_list)--;
695 **border_list = index;
696 if (lev < max_search_level) {
697 Triangle &t(triangles[lev][index]);
698 lev++;
699 index <<= 2;
700 inside_list++;
701 border_list++;
702 bool edge0_inside[nr_of_half_spaces];
703 bool edge1_inside[nr_of_half_spaces];
704 bool edge2_inside[nr_of_half_spaces];
705 for (int h=0;h<halfs_used_count;h++) {
706 const int i = halfs_used[h];
707 const HalfSpace &half_space(half_spaces[i]);
708 edge0_inside[i] = half_space.inside(t.e0);
709 edge1_inside[i] = half_space.inside(t.e1);
710 edge2_inside[i] = half_space.inside(t.e2);
711 }
712 searchZones(lev,index+0,
713 half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
714 corner0_inside,edge2_inside,edge1_inside,
715 inside_list,border_list,max_search_level);
716 searchZones(lev,index+1,
717 half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
718 edge2_inside,corner1_inside,edge0_inside,
719 inside_list,border_list,max_search_level);
720 searchZones(lev,index+2,
721 half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
722 edge1_inside,edge0_inside,corner2_inside,
723 inside_list,border_list,max_search_level);
724 searchZones(lev,index+3,
725 half_spaces,nr_of_half_spaces,halfs_used,halfs_used_count,
726 edge0_inside,edge1_inside,edge2_inside,
727 inside_list,border_list,max_search_level);
728 }
729 }
730 }
731
732
733
GeodesicSearchResult(const GeodesicGrid & grid)734 GeodesicSearchResult::GeodesicSearchResult(const GeodesicGrid &grid)
735 :grid(grid),
736 zones(new int*[grid.getMaxLevel()+1]),
737 inside(new int*[grid.getMaxLevel()+1]),
738 border(new int*[grid.getMaxLevel()+1]) {
739 for (int i=0;i<=grid.getMaxLevel();i++) {
740 zones[i] = new int[GeodesicGrid::nrOfZones(i)];
741 }
742 }
743
~GeodesicSearchResult(void)744 GeodesicSearchResult::~GeodesicSearchResult(void) {
745 for (int i=grid.getMaxLevel();i>=0;i--) {
746 delete zones[i];
747 }
748 delete border;
749 delete inside;
750 delete zones;
751 }
752
search(const HalfSpace * half_spaces,const int nr_of_half_spaces,int max_search_level)753 void GeodesicSearchResult::search(const HalfSpace *half_spaces,
754 const int nr_of_half_spaces,
755 int max_search_level) {
756 for (int i=grid.getMaxLevel();i>=0;i--) {
757 inside[i] = zones[i];
758 border[i] = zones[i]+GeodesicGrid::nrOfZones(i);
759 }
760 grid.searchZones(half_spaces,nr_of_half_spaces,inside,border,max_search_level);
761 }
762
search(const Vector & e0,const Vector & e1,const Vector & e2,const Vector & e3,int max_search_level)763 void GeodesicSearchResult::search(const Vector &e0,const Vector &e1,
764 const Vector &e2,const Vector &e3,
765 int max_search_level) {
766 HalfSpace half_spaces[4];
767 half_spaces[0].d = e0*((e1-e0)^(e2-e0));
768 if (half_spaces[0].d > 0) {
769 half_spaces[0].n = e0^e1;
770 half_spaces[1].n = e1^e2;
771 half_spaces[2].n = e2^e3;
772 half_spaces[3].n = e3^e0;
773 half_spaces[0].d = 0.0;
774 half_spaces[1].d = 0.0;
775 half_spaces[2].d = 0.0;
776 half_spaces[3].d = 0.0;
777 search(half_spaces,4,max_search_level);
778 } else {
779 half_spaces[0].n = (e1-e0)^(e2-e0);
780 search(half_spaces,1,max_search_level);
781 }
782 }
783
784
785
reset(void)786 void GeodesicSearchInsideIterator::reset(void) {
787 level = 0;
788 max_count = 1<<(max_level<<1); // 4^max_level
789 index_p = r.zones[0];
790 end_p = r.inside[0];
791 index = (*index_p) * max_count;
792 count = (index_p < end_p) ? 0 : max_count;
793 }
794
next(void)795 int GeodesicSearchInsideIterator::next(void) {
796 if (count < max_count) return index+(count++);
797 index_p++;
798 if (index_p < end_p) {
799 index = (*index_p) * max_count;
800 count = 1;
801 return index;
802 }
803 while (level < max_level) {
804 level++;
805 max_count >>= 2;
806 index_p = r.zones[level];
807 end_p = r.inside[level];
808 if (index_p < end_p) {
809 index = (*index_p) * max_count;
810 count = 1;
811 return index;
812 }
813 }
814 return -1;
815 }
816
817 //------------------------------------------------
818
819
820 #include <stdio.h>
821 #include <stdlib.h>
822 #include <string.h>
823 #include <math.h>
824
825 #include <map>
826 #include <list>
827 #include <string>
828 #include <iostream>
829
830 using namespace std;
831
832
833
834 static int restrict_output_level_min = 0;
835 static int restrict_output_level_max = MAX_LEVEL;
836
837
838
839 struct FaintStar { // 6 byte
840 int x0:18;
841 int x1:18;
842 unsigned int b_v:7;
843 unsigned int mag:5;
844 bool operator<(const FaintStar &s) const {return (mag < s.mag);}
845 } __attribute__ ((__packed__)) ;
846
847
848 struct TycStar { // 10 byte
849 int x0:20;
850 int x1:20;
851 int dx0:14;
852 int dx1:14;
853 unsigned int b_v:7;
854 unsigned int mag:5;
855 bool operator<(const TycStar &s) const {return (mag < s.mag);}
856 } __attribute__ ((__packed__));
857
858
859 struct HipStarCompressed {
860 int hip:24; // 17 bits needed
861 unsigned char comp_ids_int; // 5 bits needed
862 int x0; // 32 bits needed
863 int x1; // 32 bits needed
864 unsigned char b_v; // 8 bits needed
865 unsigned char mag; // 5 bits needed
866 unsigned short sp_int; // 14 bits needed
867 int dx0,dx1,plx;
868 } __attribute__ ((__packed__));
869
870 struct HipStar : public HipStarCompressed {
871 string component_ids,sp;
setPlxSpHipStar872 void setPlxSp(double plx,const string &sp) {
873 HipStar::plx = (int)floor(0.5+100.0*plx);
874 if (HipStar::plx < 0) HipStar::plx = 0;
875 HipStar::sp = sp;
876 }
877 bool operator<(const HipStar &s) const {return (mag < s.mag);}
878 };
879
880
881
882 struct ZoneData {
~ZoneDataZoneData883 virtual ~ZoneData(void) {}
884 virtual int getArraySize(void) const = 0;
885 virtual HipStar *add(int level,
886 int tyc1,int tyc2,int tyc3,
887 int hip,const char *component_ids,
888 double x0,double x1,double dx0,double dx1,
889 double mag,double b_v,
890 double plx,const char *sp,
891 bool &does_not_fit) = 0;
892 void writeInfoToOutput(FILE *f) const;
893 virtual void writeStarsToOutput(FILE *f) = 0;
894 Vector center;
895 Vector axis0;
896 Vector axis1;
897 };
898
899 struct HipZoneData : public ZoneData {
900 list<HipStar> stars;
getArraySizeHipZoneData901 int getArraySize(void) const {return stars.size();}
902 HipStar *add(int level,
903 int tyc1,int tyc2,int tyc3,
904 int hip,const char *component_ids,
905 double x0,double x1,double dx0,double dx1,
906 double mag,double b_v,double plx,const char *sp,
907 bool &does_not_fit);
908 void writeStarsToOutput(FILE *f);
909 };
910
911
912 static list<HipStar*> hip_index[NR_OF_HIP+1];
913
SqueezeHip(void)914 static void SqueezeHip(void) {
915 // build lookup maps
916 map<string,int> component_ids_map,sp_map;
917 for (int i=0;i<=NR_OF_HIP;i++) {
918 for (list<HipStar*>::const_iterator it(hip_index[i].begin());
919 it!=hip_index[i].end();it++) {
920 HipStar *h = *it;
921 component_ids_map[h->component_ids]++;
922 sp_map[h->sp]++;
923 }
924 }
925 int component_ids_size = 0;
926 for (map<string,int>::iterator it(component_ids_map.begin());
927 it!=component_ids_map.end();it++,component_ids_size++) {
928 it->second = component_ids_size;
929 }
930 if (component_ids_size >= 32) {
931 cerr << "SqueezeHip: too many component_ids: "
932 << component_ids_size << endl;
933 }
934 int sp_size = 0;
935 for (map<string,int>::iterator it(sp_map.begin());
936 it!=sp_map.end();it++,sp_size++) {
937 it->second = sp_size;
938 }
939 if (sp_size >= 16384) {
940 cerr << "SqueezeHip: too many sp: " << sp_size << endl;
941 }
942 // translate the strings for the hip stars into integers:
943 for (int i=0;i<=NR_OF_HIP;i++) {
944 for (list<HipStar*>::const_iterator it(hip_index[i].begin());
945 it!=hip_index[i].end();it++) {
946 HipStar *h = *it;
947 h->comp_ids_int = component_ids_map[h->component_ids];
948 h->sp_int = sp_map[h->sp];
949 }
950 }
951 // write output files for component_ids and sp
952 FILE *f = fopen(component_ids_filename,"wb");
953 if (f==0) {
954 cout << "fopen(" << component_ids_filename << ") failed" << endl;
955 exit(-1);
956 }
957 for (map<string,int>::const_iterator it(component_ids_map.begin());
958 it!=component_ids_map.end();it++,component_ids_size++) {
959 fprintf(f,"%s\n",it->first.c_str());
960 }
961 fclose(f);
962 f = fopen(sp_filename,"wb");
963 if (f==0) {
964 cout << "fopen(" << sp_filename << ") failed" << endl;
965 exit(-1);
966 }
967 for (map<string,int>::const_iterator it(sp_map.begin());
968 it!=sp_map.end();it++,sp_size++) {
969 fprintf(f,"%s\n",it->first.c_str());
970 }
971 fclose(f);
972 }
973
974 struct TycZoneData : public ZoneData {
975 list<TycStar> stars;
getArraySizeTycZoneData976 int getArraySize(void) const {return stars.size();}
977 HipStar *add(int level,
978 int tyc1,int tyc2,int tyc3,
979 int hip,const char *component_ids,
980 double x0,double x1,double dx0,double dx1,
981 double mag,double b_v,double plx,const char *sp,
982 bool &does_not_fit);
983 void writeStarsToOutput(FILE *f);
984 };
985
986 struct FaintZoneData : public ZoneData {
987 list<FaintStar> stars;
getArraySizeFaintZoneData988 int getArraySize(void) const {return stars.size();}
989 HipStar *add(int level,
990 int tyc1,int tyc2,int tyc3,
991 int hip,const char *component_ids,
992 double x0,double x1,double dx0,double dx1,
993 double mag,double b_v,double plx,const char *sp,
994 bool &does_not_fit);
995 void writeStarsToOutput(FILE *f);
996 };
997
998
999 class Accumulator {
1000 public:
1001 Accumulator(void);
1002 ~Accumulator(void);
1003 int addStar(int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
1004 double ra, // degrees
1005 double dec, // degrees
1006 double pma,double pmd,
1007 double mag,double b_v,
1008 double plx,const char *sp);
1009 void writeOutput(const char *fnames[]);
1010 private:
1011 GeodesicGrid grid;
1012 struct ZoneArray {
ZoneArrayZoneArray1013 ZoneArray(int level,GeodesicGrid &grid)
1014 : level(level),nr_of_zones(GeodesicGrid::nrOfZones(level)),grid(grid)
1015 {scale = 0.0;nr_of_stars = 0;}
~ZoneArrayZoneArray1016 virtual ~ZoneArray(void) {}
1017 virtual ZoneData *getZone(int index) const = 0;
1018 virtual void writeHeaderToOutput(FILE *f) const = 0;
1019 const int level;
1020 const int nr_of_zones;
1021 const GeodesicGrid &grid;
1022 double scale;
1023 // int scale_int;
1024 int nr_of_stars;
1025 HipStar *addStar(int tyc1,int tyc2,int tyc3,
1026 int hip,const char *component_ids,
1027 double ra, // degrees
1028 double dec, // degrees
1029 double pma,double pmd,
1030 double mag,double b_v,
1031 double plx,const char *sp,
1032 bool &does_not_fit);
1033 void writeOutput(const char *fname);
1034 };
1035 struct HipZoneArray : public ZoneArray {
HipZoneArrayHipZoneArray1036 HipZoneArray(int level,GeodesicGrid &grid)
1037 : ZoneArray(level,grid),zones(new HipZoneData[nr_of_zones]) {}
~HipZoneArrayHipZoneArray1038 ~HipZoneArray(void) {delete[] zones;}
1039 void writeHeaderToOutput(FILE *f) const;
1040 HipZoneData *const zones;
getZoneHipZoneArray1041 ZoneData *getZone(int index) const {
1042 if (index<0 || index>=nr_of_zones) {
1043 cerr << "getZone: bad index" << endl;
1044 exit(-1);
1045 }
1046 return zones+index;
1047 }
1048 };
1049 struct TycZoneArray : public ZoneArray {
TycZoneArrayTycZoneArray1050 TycZoneArray(int level,GeodesicGrid &grid)
1051 : ZoneArray(level,grid),zones(new TycZoneData[nr_of_zones]) {}
~TycZoneArrayTycZoneArray1052 ~TycZoneArray(void) {delete[] zones;}
1053 void writeHeaderToOutput(FILE *f) const;
1054 TycZoneData *zones;
getZoneTycZoneArray1055 ZoneData *getZone(int index) const {
1056 if (index<0 || index>=nr_of_zones) {
1057 cerr << "getZone: bad index" << endl;
1058 exit(-1);
1059 }
1060 return zones+index;
1061 }
1062 };
1063 struct FaintZoneArray : public ZoneArray {
FaintZoneArrayFaintZoneArray1064 FaintZoneArray(int level,GeodesicGrid &grid)
1065 : ZoneArray(level,grid),zones(new FaintZoneData[nr_of_zones]) {}
~FaintZoneArrayFaintZoneArray1066 ~FaintZoneArray(void) {delete[] zones;}
1067 void writeHeaderToOutput(FILE *f) const;
1068 FaintZoneData *zones;
getZoneFaintZoneArray1069 ZoneData *getZone(int index) const {
1070 if (index<0 || index>=nr_of_zones) {
1071 cerr << "getZone: bad index" << endl;
1072 exit(-1);
1073 }
1074 return zones+index;
1075 }
1076 };
1077 ZoneArray *zone_array[MAX_LEVEL+1];
1078
1079 private:
initTriangleFunc(int lev,int index,const Vector & c0,const Vector & c1,const Vector & c2,void * user_data)1080 static void initTriangleFunc(int lev,int index,
1081 const Vector &c0,
1082 const Vector &c1,
1083 const Vector &c2,
1084 void *user_data) {
1085 reinterpret_cast<Accumulator*>(user_data)
1086 ->initTriangle(lev,index,c0,c1,c2);
1087 }
1088 void initTriangle(int lev,int index,
1089 const Vector &c0,
1090 const Vector &c1,
1091 const Vector &c2);
1092 };
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102 static
WriteByte(FILE * f,int x)1103 void WriteByte(FILE *f,int x) {
1104 unsigned char c = (unsigned int)x;
1105 if (1!=fwrite(&c,1,1,f)) {
1106 cerr << "WriteByte: fwrite failed" << endl;
1107 exit(-1);
1108 }
1109 }
1110
1111 static
WriteInt(FILE * f,int x)1112 void WriteInt(FILE *f,int x) {
1113 if (4!=fwrite(&x,1,4,f)) {
1114 cerr << "WriteInt: fwrite failed" << endl;
1115 exit(-1);
1116 }
1117 }
1118
1119 //static
1120 //int DoubleToInt(double x) {
1121 // return (int)floor(0.5+x*0x7FFFFFFF);
1122 //}
1123
1124 //static
1125 //double IntToDouble(int x) {
1126 // double rval = x;
1127 // rval /= 0x7FFFFFFF;
1128 // return rval;
1129 //}
1130
writeInfoToOutput(FILE * f)1131 void ZoneData::writeInfoToOutput(FILE *f) const {
1132 // WriteInt(f,DoubleToInt(center.x[0]));
1133 // WriteInt(f,DoubleToInt(center.x[1]));
1134 // WriteInt(f,DoubleToInt(center.x[2]));
1135 WriteInt(f,getArraySize());
1136 }
1137
writeStarsToOutput(FILE * f)1138 void HipZoneData::writeStarsToOutput(FILE *f) {
1139 stars.sort();
1140 for (list<HipStar>::const_iterator it(stars.begin());
1141 it!=stars.end();it++) {
1142 if (sizeof(HipStarCompressed)!=
1143 fwrite(&(*it),1,sizeof(HipStarCompressed),f)) {
1144 cerr << "HipZoneData::writeStarsToOutput: fwrite failed" << endl;
1145 exit(-1);
1146 }
1147 }
1148 }
1149
writeStarsToOutput(FILE * f)1150 void TycZoneData::writeStarsToOutput(FILE *f) {
1151 stars.sort();
1152 for (list<TycStar>::const_iterator it(stars.begin());
1153 it!=stars.end();it++) {
1154 if (sizeof(TycStar)!=fwrite(&(*it),1,sizeof(TycStar),f)) {
1155 cerr << "TycZoneData::writeStarsToOutput: fwrite failed" << endl;
1156 exit(-1);
1157 }
1158 }
1159 }
1160
1161
writeStarsToOutput(FILE * f)1162 void FaintZoneData::writeStarsToOutput(FILE *f) {
1163 stars.sort();
1164 for (list<FaintStar>::const_iterator it(stars.begin());
1165 it!=stars.end();it++) {
1166 if (sizeof(FaintStar)!=fwrite(&(*it),1,sizeof(FaintStar),f)) {
1167 cerr << "FaintZoneData::writeStarsToOutput: fwrite failed" << endl;
1168 exit(-1);
1169 }
1170 }
1171 }
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
Accumulator(void)1191 Accumulator::Accumulator(void)
1192 :grid(MAX_LEVEL) {
1193 int l=0;
1194 for (;l<=MAX_HIP_LEVEL;l++) {
1195 zone_array[l] = new HipZoneArray(l,grid);
1196 }
1197 for (;l<=MAX_TYC_LEVEL;l++) {
1198 zone_array[l] = new TycZoneArray(l,grid);
1199 }
1200 for (;l<=MAX_LEVEL;l++) {
1201 zone_array[l] = new FaintZoneArray(l,grid);
1202 }
1203 grid.visitTriangles(MAX_LEVEL,initTriangleFunc,this);
1204 for (int l=0;l<=MAX_LEVEL;l++) {
1205 // zone_array[l]->scale_int = DoubleToInt(zone_array[l]->scale);
1206 // while (IntToDouble(zone_array[l]->scale_int) > zone_array[l]->scale) {
1207 // zone_array[l]->scale_int--;
1208 // }
1209 // zone_array[l]->scale = IntToDouble(zone_array[l]->scale_int);
1210 cout << "Accumulator::Accumulator: level " << zone_array[l]->level
1211 << ", scale: " << zone_array[l]->scale << endl;
1212 }
1213 }
1214
~Accumulator(void)1215 Accumulator::~Accumulator(void) {
1216 for (int l=0;l<=MAX_LEVEL;l++) {
1217 delete zone_array[l];
1218 }
1219 }
1220
1221 static const Vector north = {0.0,0.0,1.0};
1222
initTriangle(int lev,int index,const Vector & c0,const Vector & c1,const Vector & c2)1223 void Accumulator::initTriangle(int lev,int index,
1224 const Vector &c0,
1225 const Vector &c1,
1226 const Vector &c2) {
1227 ZoneData &z(*zone_array[lev]->getZone(index));
1228 double &scale(zone_array[lev]->scale);
1229 z.center = c0+c1+c2;
1230 z.center.normalize();
1231 z.axis0 = north ^ z.center;
1232 z.axis0.normalize();
1233 z.axis1 = z.center ^ z.axis0;
1234
1235 double mu0,mu1,f,h;
1236
1237 mu0 = (c0-z.center)*z.axis0;
1238 mu1 = (c0-z.center)*z.axis1;
1239 f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
1240 h = fabs(mu0)*f;
1241 if (scale < h) scale = h;
1242 h = fabs(mu1)*f;
1243 if (scale < h) scale = h;
1244
1245 mu0 = (c1-z.center)*z.axis0;
1246 mu1 = (c1-z.center)*z.axis1;
1247 f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
1248 h = fabs(mu0)*f;
1249 if (scale < h) scale = h;
1250 h = fabs(mu1)*f;
1251 if (scale < h) scale = h;
1252
1253 mu0 = (c2-z.center)*z.axis0;
1254 mu1 = (c2-z.center)*z.axis1;
1255 f = 1.0/sqrt(1.0-mu0*mu0-mu1*mu1);
1256 h = fabs(mu0)*f;
1257 if (scale < h) scale = h;
1258 h = fabs(mu1)*f;
1259 if (scale < h) scale = h;
1260 }
1261
1262 static
BvToColorIndex(double b_v)1263 unsigned char BvToColorIndex(double b_v) {
1264 b_v *= 1000.0;
1265 if (b_v < -500) {
1266 b_v = -500;
1267 } else if (b_v > 3499) {
1268 b_v = 3499;
1269 }
1270 return (unsigned int)floor(0.5+127.0*((500.0+b_v)/4000.0));
1271 }
1272
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)1273 HipStar *FaintZoneData::add(int level,
1274 int tyc1,int tyc2,int tyc3,
1275 int hip,const char *component_ids,
1276 double x0,double x1,double dx0,double dx1,
1277 double mag,double b_v,double plx,const char *sp,
1278 bool &does_not_fit) {
1279 if (mag>=level_mag_limit[level]) {
1280 cout << "too faint" << endl;
1281 exit(-1);
1282 }
1283 FaintStar s;
1284 s.x0 = ((int)floor(0.5+x0*((1<<17)-1))); // 18 bit signed
1285 s.b_v = BvToColorIndex(b_v); // 0..127: 7 Bit unsigned
1286 s.x1 = ((int)floor(0.5+x1*((1<<17)-1))); // 18 bit signed
1287 s.mag = (unsigned int)floor(30*(mag-level_mag_limit[level-1])/
1288 (level_mag_limit[level]-level_mag_limit[level-1]));
1289 // steps(accuracy) of 0.05mag, 5 bit unsigned
1290 stars.push_back(s);
1291 does_not_fit = false;
1292 return 0;
1293 }
1294
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)1295 HipStar *TycZoneData::add(int level,
1296 int tyc1,int tyc2,int tyc3,
1297 int hip,const char *component_ids,
1298 double x0,double x1,double dx0,double dx1,
1299 double mag,double b_v,double plx,const char *sp,
1300 bool &does_not_fit) {
1301 TycStar s;
1302 s.x0 = ((int)floor(0.5+x0*((1<<19)-1))); // 20 bit signed
1303 s.b_v = BvToColorIndex(b_v); // 0..127: 7 Bit unsigned
1304 s.x1 = ((int)floor(0.5+x1*((1<<19)-1))); // 20 bit signed
1305 s.mag = (unsigned int)floor(30*(mag-level_mag_limit[level-1])/
1306 (level_mag_limit[level]-level_mag_limit[level-1]));
1307 // steps(accuracy) of 0.05mag, 5 bit unsigned
1308 const int dx0_int = (int)floor(0.5+dx0*10.0); // 14 bit signed
1309 const int dx1_int = (int)floor(0.5+dx1*10.0); // 14 bit signed
1310 if (dx0_int >= 8192 || dx0_int<-8192 ||
1311 dx1_int >= 8192 || dx1_int<-8192) {
1312 // does not fit, must store in Hip format
1313 does_not_fit = true;
1314 return 0;
1315 }
1316 s.dx0 = dx0_int;
1317 s.dx1 = dx1_int;
1318 stars.push_back(s);
1319 does_not_fit = false;
1320 return 0;
1321 }
1322
1323
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)1324 HipStar *HipZoneData::add(int level,
1325 int tyc1,int tyc2,int tyc3,
1326 int hip,const char *component_ids,
1327 double x0,double x1,double dx0,double dx1,
1328 double mag,double b_v,double plx,const char *sp,
1329 bool &does_not_fit) {
1330 HipStar s;
1331 s.x0 = ((int)floor(0.5+x0*((1<<31)-1))); // 32 bit signed
1332 s.b_v = BvToColorIndex(b_v); // 0..127: 7 Bit unsigned
1333 s.x1 = ((int)floor(0.5+x1*((1<<31)-1))); // 32 bit signed
1334 const int mag_int = (int)floor(20*
1335 (mag-((level==0)?-2.0:level_mag_limit[level-1])));
1336 if (mag_int < 0 || mag_int > 255) {
1337 cerr << "HipZoneData(" << level << ")::add: bad mag: " << mag << endl;
1338 exit(-1);
1339 }
1340 s.mag = mag_int;
1341 s.dx0 = (int)floor(0.5+dx0*10.0);
1342 s.dx1 = (int)floor(0.5+dx1*10.0);
1343 s.hip = hip;
1344 s.component_ids = component_ids;
1345 s.comp_ids_int = 0;
1346 s.sp_int = 0;
1347 s.plx = (int)floor(0.5+100.0*plx);
1348 s.sp = sp;
1349 stars.push_back(s);
1350 does_not_fit = false;
1351 return &(stars.back());
1352 }
1353
1354
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)1355 HipStar *Accumulator::ZoneArray::addStar(
1356 int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
1357 double ra,double dec,double pma,double pmd,
1358 double mag,double b_v,double plx,const char *sp,
1359 bool &does_not_fit) {
1360 ra *= (M_PI/180.0);
1361 dec *= (M_PI/180.0);
1362 if (ra < 0 || ra >= 2*M_PI || dec < -0.5*M_PI || dec > 0.5*M_PI) {
1363 cerr << "ZoneArray(l=" << level << ")::addStar: "
1364 "bad ra/dec: " << ra << ',' << dec << endl;
1365 exit (-1);
1366 }
1367 const double cdec = cos(dec);
1368 Vector pos;
1369 pos.x[0] = cos(ra)*cdec;
1370 pos.x[1] = sin(ra)*cdec;
1371 pos.x[2] = sin(dec);
1372 const int zone = grid.searchZone(pos,level);
1373
1374 ZoneData &z(*getZone(zone));
1375 const double mu0 = (pos-z.center) * z.axis0;
1376 const double mu1 = (pos-z.center) * z.axis1;
1377 const double d = 1.0 - mu0*mu0 - mu1*mu1;
1378 const double sd = sqrt(d);
1379 const double x0 = mu0/(sd*scale);
1380 const double x1 = mu1/(sd*scale);
1381 const double h = mu0*pma + mu1*pmd;
1382 const double dx0 = (pma*d + mu0*h) / (sd*d);
1383 const double dx1 = (pmd*d + mu1*h) / (sd*d);
1384
1385 //cout << "Accumulator::ZoneArray(l=" << level << ")::addStar: " << zone << endl;
1386 HipStar *rval = z.add(level,tyc1,tyc2,tyc3,hip,component_ids,
1387 x0,x1,dx0,dx1,mag,b_v,plx,sp,does_not_fit);
1388 if (!does_not_fit) nr_of_stars++;
1389 //cout << "Accumulator::ZoneArray(l=" << level << ")::addStar: 999";
1390 return rval;
1391 }
1392
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)1393 int Accumulator::addStar(int tyc1,int tyc2,int tyc3,int hip,const char *component_ids,
1394 double ra,double dec,double pma,double pmd,
1395 double mag,double b_v,double plx,const char *sp) {
1396 // const int packed_hip = PackHip(hip,component_ids);
1397 // const int packed_tyc = PackTyc(tyc1,tyc2,tyc3);
1398
1399 int l=0;
1400 while (l<MAX_LEVEL && mag>=level_mag_limit[l]) l++;
1401 if (hip>0 && l>MAX_HIP_LEVEL) l = MAX_HIP_LEVEL;
1402 // store the faint tyc2 stars as faint stars
1403 // else if (tyc1>0 && l>MAX_TYC_LEVEL) l = MAX_TYC_LEVEL;
1404
1405
1406 if (restrict_output_level_max < l ||
1407 restrict_output_level_min > l) {
1408 // you may want to restrict the output
1409 // in order to restrict the amount of main memory used.
1410 // Calling the program twice may be faster than calling once
1411 // without enough main memory (avoid swapping).
1412 return 0;
1413 }
1414
1415 bool does_not_fit;
1416 HipStar *s = zone_array[l]->addStar(tyc1,tyc2,tyc3,hip,component_ids,
1417 ra,dec,pma,pmd,
1418 mag,b_v,plx,sp,does_not_fit);
1419 if (does_not_fit) s = zone_array[MAX_HIP_LEVEL]
1420 ->addStar(tyc1,tyc2,tyc3,hip,component_ids,
1421 ra,dec,pma,pmd,
1422 mag,b_v,plx,sp,does_not_fit);
1423 if (s) hip_index[hip].push_back(s);
1424 return 0;
1425 }
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440 #define FILE_MAGIC 0x835f040a
1441
writeHeaderToOutput(FILE * f)1442 void Accumulator::HipZoneArray::writeHeaderToOutput(FILE *f) const {
1443 cout << "Accumulator::HipZoneArray(" << level << ")::writeHeaderToOutput: "
1444 << nr_of_stars << endl;
1445 WriteInt(f,FILE_MAGIC);
1446 WriteInt(f,0); // type
1447 WriteInt(f,0); // major version
1448 WriteInt(f,4); // minor version
1449 WriteInt(f,level);
1450 // WriteInt(f,scale_int);
1451 if (level == 0) {
1452 WriteInt(f,-2000); // min_mag
1453 } else {
1454 WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
1455 }
1456 WriteInt(f,12800); // mag_range
1457 WriteInt(f,256); // mag_steps
1458 }
1459
writeHeaderToOutput(FILE * f)1460 void Accumulator::TycZoneArray::writeHeaderToOutput(FILE *f) const {
1461 cout << "Accumulator::TycZoneArray(" << level << ")::writeHeaderToOutput: "
1462 << nr_of_stars << endl;
1463 WriteInt(f,FILE_MAGIC);
1464 WriteInt(f,1); // type
1465 WriteInt(f,0); // major version
1466 WriteInt(f,3); // minor version
1467 WriteInt(f,level);
1468 // WriteInt(f,scale_int);
1469 WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
1470 WriteInt(f,(int)floor(0.5+1000.0*(level_mag_limit[level]
1471 -level_mag_limit[level-1]))); // mag_range
1472 WriteInt(f,30); // mag_steps
1473 }
1474
writeHeaderToOutput(FILE * f)1475 void Accumulator::FaintZoneArray::writeHeaderToOutput(FILE *f) const {
1476 cout << "Accumulator::FaintZoneArray(" << level << ")::writeHeaderToOutput: "
1477 << nr_of_stars << endl;
1478 WriteInt(f,FILE_MAGIC);
1479 WriteInt(f,2); // type
1480 WriteInt(f,0); // major version
1481 WriteInt(f,1); // minor version
1482 WriteInt(f,level);
1483 // WriteInt(f,scale_int);
1484 WriteInt(f,(int)floor(0.5+1000.0*level_mag_limit[level-1])); // min_mag
1485 WriteInt(f,(int)floor(0.5+1000.0*(level_mag_limit[level]
1486 -level_mag_limit[level-1]))); // mag_range
1487 WriteInt(f,30); // mag_steps
1488 }
1489
1490
writeOutput(const char * fname)1491 void Accumulator::ZoneArray::writeOutput(const char *fname) {
1492 if (nr_of_stars <= 0) return;
1493 FILE *f = fopen(fname,"wb");
1494 if (f==0) {
1495 cerr << "Accumulator::writeOutput(" << fname << "): "
1496 "fopen failed" << endl;
1497 return;
1498 }
1499 writeHeaderToOutput(f);
1500 for (int zone=0;zone<nr_of_zones;zone++) {
1501 getZone(zone)->writeInfoToOutput(f);
1502 }
1503 int permille = 0;
1504 for (int zone=0;zone<nr_of_zones;zone++) {
1505 getZone(zone)->writeStarsToOutput(f);
1506 int p = (1000*zone)/nr_of_zones;
1507 if (p != permille) {
1508 permille = p;
1509 cout << "\rAccumulator::writeOutput(" << fname << "): "
1510 << permille << "permille written" << flush;
1511 }
1512 }
1513 cout << endl;
1514 fclose(f);
1515 }
1516
writeOutput(const char * fnames[])1517 void Accumulator::writeOutput(const char *fnames[]) {
1518 for (int l=0;l<=MAX_LEVEL;l++) {
1519 zone_array[l]->writeOutput(fnames[l]);
1520 }
1521 }
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
ReadHipTycFile(Accumulator & accu)1537 int ReadHipTycFile(Accumulator &accu) {
1538 int count = 0;
1539 FILE *f;
1540 const char *fname = "HipTyc";
1541 f = fopen(fname,"r");
1542 if (f == 0) {
1543 fprintf(stderr,"Could not open file \"%s\".\n",fname);
1544 exit(-1);
1545 }
1546
1547 int hip,tyc1,tyc2,tyc3;
1548 char cids[32];
1549 char sp[256];
1550 int mag,b_v,VarFlag;
1551 double ra,dec,Plx,pm_ra,pm_dec;
1552 while (14==fscanf(f,"%d%d%d%d%s%d%lf%lf%lf%lf%lf%d%d%s",
1553 &hip,&tyc1,&tyc2,&tyc3,cids,&VarFlag,
1554 &ra,&dec,&Plx,&pm_ra,&pm_dec,&mag,&b_v,sp)) {
1555 if (b_v>-500 && b_v<3450) {
1556 const int rc = accu.addStar(tyc1,tyc2,tyc3,hip,
1557 cids[0]=='?'?"":cids,
1558 ra, // degrees
1559 dec, // degrees
1560 pm_ra,pm_dec,0.001*mag,0.001*b_v,
1561 Plx,sp[0]=='?'?"":sp);
1562 if (rc < 0) {
1563 // never mind: propably no magnitude for Hiparcos star
1564 // fprintf(stderr,"File \"%s\", record %d: Error 13 %d %d \"%s\"\n",
1565 // fname,count,rc,hip,sp);
1566 // exit(-1);
1567 }
1568 count++;
1569 }
1570 }
1571 fclose(f);
1572 return count;
1573 }
1574
1575
1576 const unsigned short int SHORT_ASTSRCBIT0 = 0x0001; // Astrometry source bit 0
1577 const unsigned short int SHORT_ASTSRCBIT1 = 0x0002; // Astrometry source bit 1
1578 const unsigned short int SHORT_ASTSRCBIT2 = 0x0004; // Astrometry source bit 2
1579 const unsigned short int SHORT_UBBIT = 0x0008;
1580 const unsigned short int SHORT_TMBIT = 0x0010;
1581 const unsigned short int SHORT_XRBIT = 0x0020;
1582 const unsigned short int SHORT_IUCBIT = 0x0040;
1583 const unsigned short int SHORT_ITYBIT = 0x0080;
1584 const unsigned short int SHORT_OMAGBIT = 0x0100;
1585 const unsigned short int SHORT_EMAGBIT = 0x0200;
1586 const unsigned short int SHORT_TMONLY = 0x0400;
1587 const unsigned short int SHORT_SPIKE = 0x0800;
1588 const unsigned short int SHORT_TYCONF = 0x1000;
1589 const unsigned short int SHORT_BSCONF = 0x2000;
1590 const unsigned short int SHORT_BSART = 0x4000;
1591 const unsigned short int SHORT_USEME = 0x8000;
1592
1593 struct Short_NOMAD_Record {
1594 int ra,spd,pm_ra,pm_spd;
1595 short int b,v,r;
1596 unsigned short int flags;
1597 };
1598
1599
1600
1601 #define READ_SIZE 10000
1602 static Short_NOMAD_Record buff[READ_SIZE];
1603
ReadNOMADFile(const char * fname,Accumulator & accu)1604 void ReadNOMADFile(const char *fname,Accumulator &accu) {
1605 int count;
1606 FILE *f;
1607 f = fopen(fname,"r");
1608 if (f == 0) {
1609 fprintf(stderr,"Could not open file \"%s\".\n",fname);
1610 exit(-1);
1611 }
1612 int total = 0;
1613 do {
1614 count = fread(buff,sizeof(Short_NOMAD_Record),READ_SIZE,f);
1615 total += count;
1616 printf("\rfread(%s,...) returned %6d, total = %8d",
1617 fname,count,total);
1618 fflush(stdout);
1619 int i;
1620 for (i=0;i<count;i++) {
1621 int mag = 30000;
1622 int b_v;
1623 if (buff[i].v >= 30000) {
1624 if (buff[i].b >= 30000) {
1625 if (buff[i].r >= 30000) {
1626 // cerr << "no magnitude at all" << endl;
1627 } else {
1628 mag = buff[i].r + 3499; // just an assumption
1629 b_v = 3499;
1630 }
1631 } else {
1632 mag = buff[i].b + 500; // just an assumption
1633 b_v = -500;
1634 }
1635 } else {
1636 mag = buff[i].v;
1637 if (buff[i].b >= 30000) {
1638 if (buff[i].r >= 30000) {
1639 b_v = 0;
1640 } else {
1641 b_v = buff[i].v-buff[i].r; // desperate
1642 }
1643 } else {
1644 b_v = buff[i].b-buff[i].v;
1645 }
1646 }
1647
1648 int nr_of_measurements = 0;
1649 if (buff[i].b < 30000) nr_of_measurements++;
1650 if (buff[i].v < 30000) nr_of_measurements++;
1651 if (buff[i].r < 30000) nr_of_measurements++;
1652
1653 if (mag < 19500 && b_v>-500 && b_v<3450 &&
1654 ((buff[i].flags&SHORT_USEME) ||
1655 (
1656 ((buff[i].flags&(SHORT_ASTSRCBIT0|SHORT_ASTSRCBIT1|SHORT_ASTSRCBIT2))!=1 ||
1657 (((buff[i].flags&SHORT_UBBIT)==0)
1658 &&
1659 (mag>14000 || nr_of_measurements>1) &&
1660 (mag>13000 || nr_of_measurements>2)
1661 )) &&
1662 ((buff[i].flags&SHORT_SPIKE)==0) &&
1663 ((buff[i].flags&SHORT_BSART)==0) &&
1664 ((buff[i].flags&SHORT_TMONLY)==0)))) {
1665 if (accu.addStar(0,0,0,0,"",
1666 buff[i].ra/(3600.0*1000.0), // degrees
1667 (buff[i].spd-90*3600*1000)/(3600.0*1000.0), // degrees
1668 0.1*buff[i].pm_ra,0.1*buff[i].pm_spd,
1669 0.001*mag,0.001*b_v,
1670 0,"") < 0) {
1671 fprintf(stderr,"File \"%s\", record %d: Error 16\n",fname,count);
1672 exit(-1);
1673 }
1674 }
1675 }
1676 } while (count == READ_SIZE);
1677 printf("\n");
1678 fclose(f);
1679 }
1680
1681
main(int argc,char * argv[])1682 int main(int argc,char *argv[]) {
1683 if (argc != 3 ||
1684 1 != sscanf(argv[1],"%d",&restrict_output_level_min) ||
1685 1 != sscanf(argv[2],"%d",&restrict_output_level_max)) {
1686 cerr << "Usage: " << argv[0] << " level_min level_max" << endl
1687 << " (like " << argv[0] << " 0 6)" << endl << endl;
1688 return -1;
1689 }
1690
1691 Accumulator accu;
1692 int n=0;
1693 n = ReadHipTycFile(accu);
1694 printf("HipTyc: %d records processed.\n",n);
1695 SqueezeHip();
1696
1697 for (int c=0;nomad_names[c];c++) {
1698 ReadNOMADFile(nomad_names[c],accu);
1699 }
1700
1701 accu.writeOutput(output_file_names);
1702
1703 return 0;
1704 }
1705