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 <math.h>
9
10 #include "os-features.h"
11 #include "scamp-catalog.h"
12 #include "usnob-fits.h"
13 //#include "healpix-utils.h"
14 #include "healpix.h"
15 #include "starutil.h"
16 #include "errors.h"
17 #include "log.h"
18
19 const char* OPTIONS = "hu:o:A:D:r:n:RBNv";
20
print_help(char * progname)21 void print_help(char* progname) {
22 printf("Usage: %s\n"
23 " -u <usnob-fits-filename-pattern>: eg /path/to/usnob/usnob_hp%%03i.fits\n"
24 " -n <usnob-fits-nside> : Nside of USNOB healpixelization.\n"
25 " -o <scamp-reference-catalog> : output filename\n"
26 " -A <RA center>\n"
27 " -D <Dec center>\n"
28 " -r <radius in arcmin>\n"
29 " [-R] : use Red mags\n"
30 " [-B] : use Blue mags\n"
31 " [-N] : use Infrared mags\n"
32 " [-v]: verbose\n"
33 "\n", progname);
34 }
35
36
main(int argc,char ** args)37 int main(int argc, char** args) {
38 int c;
39 char* usnobpath = NULL;
40 char* scampref = NULL;
41 double ra = 0.0;
42 double dec = 0.0;
43 double radius = 0.0;
44 double xyz[3];
45 double range;
46 int healpixes[9];
47 int nhp;
48 int nside;
49 int i;
50 scamp_cat_t* scamp;
51 anbool red, blue, infrared;
52 int loglvl = LOG_MSG;
53
54 red = blue = infrared = FALSE;
55
56 while ((c = getopt(argc, args, OPTIONS)) != -1) {
57 switch (c) {
58 case 'h':
59 print_help(args[0]);
60 exit(0);
61 case 'u':
62 usnobpath = optarg;
63 break;
64 case 'n':
65 nside = atoi(optarg);
66 break;
67 case 'o':
68 scampref = optarg;
69 break;
70 case 'A':
71 ra = atora(optarg);
72 break;
73 case 'D':
74 dec = atodec(optarg);
75 break;
76 case 'r':
77 radius = atof(optarg);
78 break;
79 case 'R':
80 red = TRUE;
81 break;
82 case 'B':
83 blue = TRUE;
84 break;
85 case 'N':
86 infrared = TRUE;
87 break;
88 case 'v':
89 loglvl++;
90 break;
91 }
92 }
93 log_init(loglvl);
94
95 if (ra == HUGE_VAL || dec == HUGE_VAL || !usnobpath || !scampref || radius == 0.0 || !nside) {
96 print_help(args[0]);
97 printf("\n\nNeed RA, Dec, USNOB path, Nside, Scamp output file, and radius.\n");
98 exit(-1);
99 }
100
101 if ((red ? 1:0) + (blue ? 1:0) + (infrared ? 1:0) != 1) {
102 print_help(args[0]);
103 printf("Must select exactly one of Red, Blue and Infrared (-R, -B, -N)\n");
104 exit(-1);
105 }
106
107 logverb("(RA,Dec) center (%g, %g) degrees\n", ra, dec);
108 logverb("Search radius %g arcmin\n", radius);
109
110 scamp = scamp_catalog_open_for_writing(scampref, TRUE);
111 if (!scamp ||
112 scamp_catalog_write_field_header(scamp, NULL)) {
113 ERROR("Failed to open SCAMP reference catalog for writing: \"%s\"", scampref);
114 exit(-1);
115 }
116
117 radecdeg2xyzarr(ra, dec, xyz);
118 range = arcmin2dist(radius);
119 nhp = healpix_get_neighbours_within_range(xyz, range, healpixes, nside);
120
121 logverb("Found %i healpixes within range.\n", nhp);
122
123 for (i=0; i<nhp; i++) {
124 int hp = healpixes[i];
125 char* path;
126 usnob_fits* usnob;
127 int j, N;
128 int nspikes = 0;
129 int nanspikes = 0;
130 int ntycho = 0;
131 int noutofrange = 0;
132 int nnomag = 0;
133 int nwritten = 0;
134
135 asprintf(&path, usnobpath, hp);
136 logmsg("Opening USNOB file \"%s\"\n", path);
137 usnob = usnob_fits_open(path);
138 if (!usnob) {
139 ERROR("Failed to open USNOB file \"%s\"", path);
140 exit(-1);
141 }
142 N = usnob_fits_count_entries(usnob);
143 logmsg("Reading %i entries.\n", N);
144 for (j=0; j<N; j++) {
145 usnob_entry* entry;
146 scamp_ref_t ref;
147 float mag;
148
149 entry = usnob_fits_read_entry(usnob);
150
151 if (!usnob_is_usnob_star(entry)) {
152 ntycho++;
153 continue;
154 }
155 if (distsq_between_radecdeg(ra, dec, entry->ra, entry->dec) > range*range) {
156 noutofrange++;
157 continue;
158 }
159 if (entry->diffraction_spike) {
160 nspikes++;
161 continue;
162 }
163 if (entry->an_diffraction_spike) {
164 nanspikes++;
165 continue;
166 }
167 if ((red && usnob_get_red_mag (entry, &mag)) ||
168 (blue && usnob_get_blue_mag (entry, &mag)) ||
169 (infrared && usnob_get_infrared_mag(entry, &mag))) {
170 nnomag++;
171 continue;
172 }
173 ref.ra = entry->ra;
174 ref.dec = entry->dec;
175 ref.err_a = entry->sigma_ra;
176 ref.err_b = entry->sigma_dec;
177 ref.mag = mag;
178 // from USNOB docs.
179 ref.err_mag = 0.25;
180 scamp_catalog_write_reference(scamp, &ref);
181 nwritten++;
182 }
183 usnob_fits_close(usnob);
184
185 logmsg("Read a total of %i USNOB entries.\n", N);
186 logmsg("Rejected %i Tycho-2 stars.\n", ntycho);
187 logmsg("Rejected %i stars that were out of range.\n", noutofrange);
188 logmsg("Rejected %i diffraction spikes (marked by USNOB).\n", nspikes);
189 logmsg("Rejected %i diffraction spikes (marked by Astrometry.net).\n", nanspikes);
190 logmsg("Rejected %i stars that were missing magnitude measurements.\n", nnomag);
191 logmsg("Wrote %i stars.\n", nwritten);
192 }
193
194 if (scamp_catalog_close(scamp)) {
195 ERROR("Failed to close SCAMP reference catalog \"%s\"", scampref);
196 exit(-1);
197 }
198
199 return 0;
200 }
201
202