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