1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <unistd.h>
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <math.h>
10 #include <string.h>
11
12 #include "healpix.h"
13 #include "starutil.h"
14 #include "rdlist.h"
15 #include "bl.h"
16
17 char* OPTIONS = "hqf:N:";
18
printHelp(char * progname)19 void printHelp(char* progname) {
20 fprintf(stderr, "Usage: %s [options]\n"
21 " -f <rdls-file>\n"
22 " [-N nside] (default 1)\n"
23 " [-h] print help msg\n"
24 " [-q] quiet mode\n",
25 progname);
26 }
27
28
main(int argc,char ** args)29 int main(int argc, char** args) {
30 char* filename = NULL;
31 int npoints;
32 int i, j;
33 int* healpixes;
34 int argchar;
35 char* progname = args[0];
36 il** lists;
37 anbool quiet = FALSE;
38 rdlist* rdls;
39 int Nside = 1;
40 int N;
41
42 while ((argchar = getopt (argc, args, OPTIONS)) != -1)
43 switch (argchar) {
44 case 'N':
45 Nside = atoi(optarg);
46 break;
47 case 'f':
48 filename = optarg;
49 break;
50 case 'h':
51 printHelp(progname);
52 exit(0);
53 case 'q':
54 quiet = TRUE;
55 break;
56 case '?':
57 fprintf(stderr, "Unknown option `-%c'.\n", optopt);
58 default:
59 exit(-1);
60 }
61
62 if (!filename) {
63 printHelp(progname);
64 exit(-1);
65 }
66
67 fprintf(stderr, "Opening RDLS file %s...\n", filename);
68 rdls = rdlist_open(filename);
69 if (!rdls) {
70 fprintf(stderr, "Failed to open RDLS file.\n");
71 exit(-1);
72 }
73
74 N = 12 * Nside * Nside;
75
76 healpixes = malloc(N * sizeof(int));
77 lists = calloc(N, sizeof(il*));
78
79 /*
80 for (i=0; i<N; i++) {
81 lists[i] = il_new(256);
82 }
83 */
84
85 for (j=1; j<=rdls_n_fields(rdls); j++) {
86 rd* points;
87
88 points = rdlist_get_field(rdls, j);
89 if (!points) {
90 fprintf(stderr, "error reading field %i\n", j);
91 break;
92 }
93
94 memset(healpixes, 0, N * sizeof(int));
95
96 npoints = rd_size(points);
97
98 for (i=0; i<npoints; i++) {
99 double ra, dec;
100 int hp;
101
102 ra = deg2rad(rd_refra (points, i));
103 dec = deg2rad(rd_refdec(points, i));
104
105 if (Nside > 1)
106 hp = radectohealpix_nside(ra, dec, Nside);
107 else
108 hp = radectohealpix(ra, dec);
109 if ((hp < 0) || (hp >= N)) {
110 printf("hp=%i\n", hp);
111 continue;
112 }
113 healpixes[hp] = 1;
114 }
115 if (!quiet) {
116 printf("Field %i: healpixes ", j);
117 for (i=0; i<N; i++) {
118 if (healpixes[i])
119 printf("%i ", i);
120 }
121 printf("\n");
122 fflush(stdout);
123 }
124
125 for (i=0; i<N; i++)
126 if (healpixes[i]) {
127 if (!lists[i])
128 lists[i] = il_new(256);
129 il_append(lists[i], j);
130 }
131
132 free_rd(points);
133 }
134
135 for (i=0; i<N; i++) {
136 int N;
137 if (!lists[i])
138 continue;
139 printf("HP %i: ", i);
140 N = il_size(lists[i]);
141 for (j=0; j<N; j++)
142 printf("%i ", il_get(lists[i], j));
143 il_free(lists[i]);
144 printf("\n");
145 }
146
147 free(lists);
148 free(healpixes);
149
150 rdlist_close(rdls);
151 return 0;
152 }
153