1
2 /*
3 * match: a package to match lists of stars (or other items)
4 * Copyright (C) 2000 Michael William Richmond
5 *
6 * Contact: Michael William Richmond
7 * Physics Department
8 * Rochester Institute of Technology
9 * 85 Lomb Memorial Drive
10 * Rochester, NY 14623-5603
11 * E-mail: mwrsps@rit.edu
12 *
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
27 *
28 */
29
30 /*
31 * <AUTO>
32 * FILE: project_coords.c
33 *
34 * <HTML>
35 * Project (RA, Dec) coords of a list of stars around some central point,
36 * creating a list with "plate coordinates" xi and eta, corresponding
37 * to the positions of the stars on a tangent plane projection.
38 *
39 * Given
40 * - an ASCII file consisting of a list of stars, with one line per
41 * star and multiple columns of information separated by white space
42 * - the numbers of the columns containing RA and Dec coords of stars
43 * (in decimal degrees)
44 * - a central RA and Dec, each in decimal degrees
45 *
46 * run through the data file. For each entry, calculate the projected
47 * coords (xi, eta) of the star, then replace the (RA, Dec) values with
48 * these projected values. Leave all other information in the ASCII
49 * file as-is.
50 *
51 * Print the results to stdout, or place them into the file given
52 * by the optional "outfile" command-line argument.
53 *
54 * Usage: project_coords starfile1 xcol ycol ra dec [outfile=]
55 *
56 * 6/28/2001: added 10 missing "%s" in the sscanf format string
57 * used to read in data in "proc_star_file".
58 * MWR
59 *
60 * 10/25/2003: added new command-line options "asec" and "arcsec",
61 * which cause the output (xi, eta) values to be
62 * converted from radians to arcseconds before being
63 * printed to output. Thanks to John Blakeslee and
64 * the ACS team.
65 *
66 * </HTML>
67 * </AUTO>
68 */
69
70
71 #include <stdio.h>
72 #include <stdlib.h>
73 #include <string.h>
74 #include <math.h>
75 #include "misc.h"
76 #include "degtorad.h"
77 #define RADtoASEC (3600.*180./M_PI)
78
79 #undef DEBUG /* get some of diagnostic output */
80
81 static int
82 proc_star_file(char *file, int racol, int deccol, double ra, double dec,
83 FILE *out_fp, int doASEC);
84
85
86 #define USAGE "usage: project_coords starfile1 racol deccol ra dec [outfile=] [asec]"
87 char *progname = "project_coords";
88
89
90 int
main(int argc,char * argv[])91 main
92 (
93 int argc,
94 char *argv[]
95 )
96 {
97 char *fileA;
98 int i, doASEC=0;
99 int racol, deccol;
100 double ra, dec;
101 char outfile[100];
102 FILE *fp;
103 TRANS trans;
104
105 /* make special search for --version, --help, and help */
106 if ((argc < 2) ||
107 (strncmp(argv[1], "--help", 6) == 0) ||
108 (strncmp(argv[1], "help", 4) == 0)) {
109 printf("%s\n", USAGE);
110 exit(0);
111 }
112 if (strncmp(argv[1], "--version", 9) == 0) {
113 printf("%s: version %s\n", progname, VERSION);
114 exit(0);
115 }
116
117 if (argc < 5) {
118 fprintf(stderr, "%s\n", USAGE);
119 exit(1);
120 }
121
122 fp = stdout;
123
124 /* parse the arguments */
125 fileA = argv[1];
126 if (sscanf(argv[2], "%d", &racol) != 1) {
127 shFatal("invalid argument for column for RA values in file");
128 }
129 if (sscanf(argv[3], "%d", &deccol) != 1) {
130 shFatal("invalid argument for column for Dec values in file");
131 }
132 if (sscanf(argv[4], "%lf", &ra) != 1) {
133 shFatal("invalid argument for RA");
134 }
135 if (sscanf(argv[5], "%lf", &dec) != 1) {
136 shFatal("invalid argument for Dec");
137 }
138
139 /*
140 * check for optional arguments "outfile" and "asec"
141 */
142 for (i = 6; i < argc; i++) {
143 if (strncmp(argv[i], "outfile=", 8) == 0) {
144 if (sscanf(argv[i] + 8, "%s", outfile) != 1) {
145 shFatal("invalid argument for outfile argument");
146 }
147 if ((fp = fopen(outfile, "w")) == NULL) {
148 shFatal("can't open file %s for output", outfile);
149 }
150 }
151 else if ((strncmp(argv[i],"asec",4) == 0) ||
152 (strncmp(argv[i],"arcsec",6) == 0)) {
153 doASEC=1;
154 #ifdef DEBUG
155 printf("Will output delta coords in projected ARCSEC, not RAD.\n");
156 #endif
157 }
158 else {
159 /* this isn't any known argument. Complain and quit */
160 shFatal("Invalid argument: %s", argv[i]);
161 }
162 }
163
164 /* now walk through the file and do the dirty work */
165 if (proc_star_file(fileA, racol, deccol, ra, dec, fp, doASEC) != SH_SUCCESS) {
166 shError("can't process data from file %s", fileA);
167 exit(1);
168 }
169
170
171 if (fp != stdout) {
172 fclose(fp);
173 }
174
175 return(0);
176 }
177
178
179
180 /****************************************************************************
181 * ROUTINE: proc_star_file
182 *
183 * walk through the given file, one line at a time.
184 *
185 * If the line starts with COMMENT_CHAR, place it into the output stream.
186 * If the line is completely blank, place it into the output stream.
187 *
188 * Otherwise,
189 * - read in the entire line,
190 * - figure out the RA and Dec coords of the star
191 * - calculate the (xi, eta) coords of the star projected onto the
192 * plane tangent to the sky at (central_ra, central_dec)
193 * - if the "doASEC" arg is non-zero, convert the values of (xi, eta)
194 * from radians to arcseconds
195 * - print out the line, replacing the "RA" with xi, the "Dec" with eta
196 *
197 * RETURNS:
198 * SH_SUCCESS if all goes well
199 * SH_GENERIC_ERROR if not
200 */
201
202 static int
proc_star_file(char * file,int racol,int deccol,double central_ra,double central_dec,FILE * out_fp,int doASEC)203 proc_star_file
204 (
205 char *file, /* I: name of input file with star list */
206 int racol, /* I: position of column with RA positions */
207 int deccol, /* I: position of column with Dec positions */
208 double central_ra, /* I: central RA of tangent plane (degrees) */
209 double central_dec, /* I: central Dec of tangent plane (degrees) */
210 FILE *out_fp, /* I: place output into this stream */
211 int doASEC /* I: if > 0, write offsets in arcsec */
212 )
213 {
214 char line[LINELEN];
215 char col[MAX_DATA_COL + 1][MAX_COL_LENGTH + 1];
216 int i, ncol;
217 int last_column = -1;
218 double raval, decval;
219 double cent_ra_rad, cent_dec_rad;
220 double dec_rad;
221 double delta_ra;
222 double xx, yy, xi, eta;
223 FILE *in_fp;
224
225 if ((in_fp = fopen(file, "r")) == NULL) {
226 shError("proc_star_file: can't open file %s for input", file);
227 return(SH_GENERIC_ERROR);
228 }
229
230 last_column = (racol > deccol ? racol : deccol);
231 cent_ra_rad = central_ra*DEGTORAD;
232 cent_dec_rad = central_dec*DEGTORAD;
233
234 while (fgets(line, LINELEN, in_fp) != NULL) {
235
236 if (line[0] == COMMENT_CHAR) {
237 continue;
238 }
239 if (is_blank(line)) {
240 continue;
241 }
242
243 shAssert(MAX_DATA_COL >= 20);
244 ncol = sscanf(line, "%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
245 &(col[0][0]), &(col[1][0]), &(col[2][0]), &(col[3][0]),
246 &(col[4][0]), &(col[5][0]), &(col[6][0]), &(col[7][0]),
247 &(col[8][0]), &(col[9][0]),
248 &(col[10][0]), &(col[11][0]), &(col[12][0]), &(col[13][0]),
249 &(col[14][0]), &(col[15][0]), &(col[16][0]), &(col[17][0]),
250 &(col[18][0]), &(col[19][0]),
251 &(col[20][0]), &(col[21][0]), &(col[22][0]), &(col[23][0]),
252 &(col[24][0]), &(col[25][0]), &(col[26][0]), &(col[27][0]),
253 &(col[28][0]), &(col[29][0]));
254 if (last_column > ncol) {
255 shError("proc_star_file: not enough entries in following line; skipping");
256 shError(" %s", line);
257 continue;
258 }
259
260 /* now read values from each column */
261 if (get_value(col[racol], &raval) != SH_SUCCESS) {
262 shError("read_data_file: can't read RA value from %s; skipping",
263 col[racol]);
264 continue;
265 }
266 if (get_value(col[deccol], &decval) != SH_SUCCESS) {
267 shError("read_data_file: can't read Dec value from %s; skipping",
268 col[deccol]);
269 continue;
270 }
271
272 /*
273 * check to see if the central RA and this star's RA are
274 * wrapped around zero from each other
275 */
276 if ((raval < 10) && (central_ra > 350)) {
277 delta_ra = (raval + 360) - central_ra;
278 }
279 else if ((raval > 350) && (central_ra < 10)) {
280 delta_ra = (raval - 360) - central_ra;
281 }
282 else {
283 delta_ra = raval - central_ra;
284 }
285 delta_ra *= DEGTORAD;
286
287
288 /*
289 * let's transform from (delta_RA, delta_Dec) to (xi, eta),
290 */
291 dec_rad = decval*DEGTORAD;
292 xx = cos(dec_rad)*sin(delta_ra);
293 yy = sin(cent_dec_rad)*sin(dec_rad) +
294 cos(cent_dec_rad)*cos(dec_rad)*cos(delta_ra);
295 xi = (xx/yy);
296
297 xx = cos(cent_dec_rad)*sin(dec_rad) -
298 sin(cent_dec_rad)*cos(dec_rad)*cos(delta_ra);
299 eta = (xx/yy);
300
301 #ifdef DEBUG
302 printf(" xi = %10.5f, eta = %10.5f\n", xi, eta);
303 #endif
304
305 /* if desired, convert xi and eta from radians to arcsec */
306 if (doASEC > 0) {
307 xi *= RADtoASEC;
308 eta *= RADtoASEC;
309 }
310 #ifdef DEBUG
311 printf(" now xi = %10.5f, eta = %10.5f\n", xi, eta);
312 #endif
313
314 /*
315 * now build up the output line. Note that we use slightly
316 * different output formats for (xi, eta) if they are expressed
317 * in radians or arcseconds.
318 */
319 line[0] = '\0';
320 for (i = 0; i < ncol; i++) {
321 if (i == racol) {
322 if (doASEC > 0) {
323 /* write the offsets in arcsec */
324 sprintf(line, "%s %12.5f", line, xi);
325 }
326 else {
327 sprintf(line, "%s %13.6e", line, xi);
328 }
329 }
330 else if (i == deccol) {
331 if (doASEC > 0) {
332 /* write the offsets in arcsec */
333 sprintf(line, "%s %12.5f", line, eta);
334 }
335 else {
336 sprintf(line, "%s %13.6e", line, eta);
337 }
338 }
339 else {
340 sprintf(line, "%s %s", line, col[i]);
341 }
342 }
343 fprintf(out_fp, "%s\n", line);
344
345 }
346
347
348 fclose(in_fp);
349 return(SH_SUCCESS);
350 }
351