1 /***************************************************************************
2 * This File is a part of the UVS Zodiacal light model software *
3 * Copyright (C) 2007 by The Tauvex Software Team, *
4 * Indian Institute of Astrophysics, *
5 * Bangalore 560 034 *
6 * tauvex_AT_iiap.res.in *
7 * *
8 * This program is free software; you can redistribute it and/or modify *
9 * it under the terms of the GNU General Public License as published by *
10 * the Free Software Foundation; either version 2 of the License, or *
11 * (at your option) any later version. *
12 * *
13 * This program is distributed in the hope that it will be useful, *
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16 * GNU General Public License for more details. *
17 * *
18 * You should have received a copy of the GNU General Public License *
19 * along with this program; if not, write to the *
20 * Free Software Foundation, Inc., *
21 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
22 ***************************************************************************/
23 /*
24 Coordinate conversion routines cannibalized from the SkyView routines of IPAC
25 Copyright (C) 1995, California Institute of Technology.
26 U.S. Government Sponsorship under NASA Contract NAS7-918 is acknowledged.
27 */
29 #define M_PI 3.141592653589793
30 static double dtr = M_PI/180;
31 static double rtd = 180/M_PI;
32 #define debug 0 /*Do not debug*/
33 #define TRUE 1
34 #define FALSE 0
35 #define debug_file "COORDINATES.DEBUG"
36 #include "sky_model.h"
38 /*This program converts J2000 into galactic coordinates using a rotation matrix*/
39 int CONV_EQ_TO_GAL(double ra, double dec, double *gl, double *gb)
40 {
41 int i, j;
42 double x[3], xp[3], r, d, lat, lon;
43 /* The following matrix comes from _The Hipparcos & Tycho */
44 /* Catalogues: Introduction & Guide to the Data_, p 92: */
45 double A[3][3] = {
46 -.0548755604, -.8734370902, -.4838350155,
47 .4941094279, -.4448296300, .7469822445,
48 -.8676661490, -.1980763734, .4559837762 };
50 r = DEGRAD(ra);
51 d = DEGRAD(dec);
52 x[0] = cos(r)*cos(d);
53 x[1] = sin(r)*cos(d);
54 x[2] = sin(d);
56 for (i = 0; i<3; ++i) {
57 xp[i] = 0;
58 for (j = 0; j < 3; ++j) {
59 xp[i] += A[i][j]*x[j];
60 }
61 }
63 if (xp[2] > 1.)
64 xp[2] = 1;
65 if (xp[2] < -1.)
66 xp[2] = -1;
67 lat = asin(xp[2]);
68 lon = atan2(xp[1], xp[0]);
69 if (lon < 0)
70 lon += 2.*PI;
71 *gl = RADDEG(lon);
72 *gb = RADDEG(lat);
74 return(EXIT_SUCCESS);
75 }
77 /*********************************************************************/
78 /* */
79 /* EC_to_EQ */
80 /* converts ecliptic coordinates to equatorial (both in degrees) */
81 /* */
82 /*********************************************************************/
84 void
85 EC_to_EQ(coord)
86 double coord[2];
87 {
88 int i, j;
89 double ra, dec, lat, lon, x[3], xp[3];
90 double dt1, dt2, dt0;
92 static double A[3][3] = {1.000000000, 0.000000000, 0.000000000,
93 0.000000000, 0.917436945, -0.397881203,
94 0.000000000, 0.397881203, 0.917436945};
96 if (debug == TRUE)
97 printf("EC_to_EQ: in %-g %-g\n", coord[0], coord[1]);
99 lon = coord[0] * dtr;
100 lat = coord[1] * dtr;
102 x[0] = cos(lon) * cos(lat);
103 x[1] = sin(lon) * cos(lat);
104 x[2] = sin(lat);
106 for (i = 0; i < 3; ++i)
107 {
108 xp[i] = 0.;
109 for (j = 0; j < 3; ++j)
110 {
111 xp[i] += A[i][j] * x[j];
112 }
113 }
115 if (xp[2] > 1.)
116 xp[2] = 1.;
118 if (xp[2] < -1.)
119 xp[2] = -1.;
120 dt2 = xp[2];
121 dt1 = xp[1];
122 dt0 = xp[0];
123 dec = asin(dt2);
124 ra = atan2(dt1, dt0);
126 if (ra < 0.)
127 ra += 2.* M_PI;
129 coord[0] = ra / dtr;
130 coord[1] = dec / dtr;
132 if (debug == TRUE)
133 printf("EC_to_EQ: out %-g %-g\n", coord[0], coord[1]);
135 }
139 /*********************************************************************/
140 /* */
141 /* GA_to_EQ */
142 /* converts galactic coordinates to equatorial (both in degrees) */
143 /* */
144 /*********************************************************************/
146 void
147 GA_to_EQ(coord)
148 double coord[2];
149 {
150 int i, j;
151 double ra, dec, lat, lon, x[3], xp[3];
153 static double A[3][3] = {-0.066988740, 0.492728470, -0.867600820,
154 -0.872755770, -0.450346960, -0.188374600,
155 -0.483538920, 0.744584640, 0.460199790};
157 if (debug == TRUE)
158 printf("GA_to_EQ: in %-g %-g\n", coord[0], coord[1]);
160 lon = coord[0] * dtr;
161 lat = coord[1] * dtr;
163 x[0] = cos(lon) * cos(lat);
164 x[1] = sin(lon) * cos(lat);
165 x[2] = sin(lat);
167 for (i = 0; i < 3; ++i)
168 {
169 xp[i] = 0.;
170 for (j = 0; j < 3; ++j)
171 {
172 xp[i] += A[i][j] * x[j];
173 }
174 }
176 if (xp[2] > 1.)
177 xp[2] = 1.;
179 if (xp[2] < -1.)
180 xp[2] = -1.;
182 dec = asin(xp[2]);
183 ra = atan2(xp[1], xp[0]);
185 if (ra < 0.)
186 ra += 2.* M_PI;
188 coord[0] = ra * rtd;
189 coord[1] = dec * rtd;
191 if (debug == TRUE)
192 printf("GA_to_EQ: out %-g %-g\n", coord[0], coord[1]);
193 }
196 /*********************************************************************/
197 /* */
198 /* EQ_to_EC */
199 /* converts equatorial coordinates to ecliptic (both in degrees) */
200 /* */
201 /*********************************************************************/
203 void
204 EQ_to_EC(coord)
205 double coord[2];
206 {
207 int i, j;
208 double ra, dec, lat, lon, x[3], xp[3];
210 static double A[3][3] = {1.000000000, 0.000000000, 0.000000000,
211 0.000000000, 0.917436945, 0.397881203,
212 0.000000000, -0.397881203, 0.917436945};
214 if (debug == TRUE)
215 printf("EQ_to_EC: in %-g %-g\n", coord[0], coord[1]);
217 ra = coord[0] * dtr;
218 dec = coord[1] * dtr;
220 x[0] = cos(ra) * cos(dec);
221 x[1] = sin(ra) * cos(dec);
222 x[2] = sin(dec);
224 for (i = 0; i < 3; ++i)
225 {
226 xp[i] = 0.;
227 for (j = 0; j < 3; ++j)
228 {
229 xp[i] += A[i][j] * x[j];
230 }
231 }
233 if (xp[2] > 1.)
234 xp[2] = 1.;
236 if (xp[2] < -1.)
237 xp[2] = -1.;
239 lat = asin(xp[2]);
240 lon = atan2(xp[1], xp[0]);
242 if (lon < 0.)
243 lon += 2.* M_PI;
245 coord[0] = lon * rtd;
246 coord[1] = lat * rtd;
248 if (debug == TRUE)
249 printf("EQ_to_EC: out %-g %-g\n", coord[0], coord[1]);
250 }
253 /*********************************************************************/
254 /* */
255 /* GA_to_EC */
256 /* converts galactic coordinates to ecliptic (both in degrees) */
257 /* */
258 /*********************************************************************/
260 void
261 GA_to_EC(coord)
262 double coord[2];
263 {
264 if (debug == TRUE)
265 printf("GA_to_EC: in %-g %-g\n", coord[0], coord[1]);
267 GA_to_EQ(coord);
268 EQ_to_EC(coord);
270 if (debug == TRUE)
271 printf("GA_to_EC: out %-g %-g\n", coord[0], coord[1]);
272 }
276 /*********************************************************************/
277 /* */
278 /* EQ_to_GA */
279 /* converts equatorial coordinates to galactic (both in degrees) */
280 /* */
281 /*********************************************************************/
283 void
284 EQ_to_GA(coord)
285 double coord[2];
286 {
287 int i, j;
288 double ra, dec, lat, lon, x[3], xp[3];
290 static double A[3][3] = {-0.066988740, -0.872755770, -0.483538920,
291 0.492728470, -0.450346960, 0.744584640,
292 -0.867600820, -0.188374600, 0.460199790};
294 ra = coord[0] * dtr;
295 dec = coord[1] * dtr;
297 x[0] = cos(ra) * cos(dec);
298 x[1] = sin(ra) * cos(dec);
299 x[2] = sin(dec);
301 for (i = 0; i < 3; ++i)
302 {
303 xp[i] = 0.;
304 for (j = 0; j < 3; ++j)
305 {
306 xp[i] += A[i][j] * x[j];
307 }
308 }
310 if (xp[2] > 1.)
311 xp[2] = 1.;
313 if (xp[2] < -1.)
314 xp[2] = -1.;
316 lat = asin(xp[2]);
317 lon = atan2(xp[1], xp[0]);
319 if (lon < 0.)
320 lon += 2.* M_PI;
322 coord[0] = lon * rtd;
323 coord[1] = lat * rtd;
325 }
329 /*********************************************************************/
330 /* */
331 /* EC_to_GA */
332 /* converts ecliptic coordinates to galactic (both in degrees) */
333 /* */
334 /*********************************************************************/
336 void
337 EC_to_GA(coord)
338 double coord[2];
339 {
340 if (debug == TRUE)
341 printf("EC_to_GA: in %-g %-g\n", coord[0], coord[1]);
343 EC_to_EQ(coord);
344 EQ_to_GA(coord);
346 if (debug == TRUE)
347 printf("EC_to_GA: in %-g %-g\n", coord[0], coord[1]);
348 }
353 /**********************************************************************/
354 /* */
355 /* COORD_TO_STRING */
356 /* converts the coordinates to strings to be printed out. The */
357 /* variable 'sex' is used to decide whether to output in decimal */
358 /* degrees or sexigesimal. */
359 /* */
360 /**********************************************************************/
363 void
364 coord_to_string(coord, s_coord)
365 double coord[2];
366 char s_coord[2][40];
367 {
368 char sign;
369 int hours, min, degs;
370 double dsec, c_hours, c_degs;
371 int sex = TRUE;
372 int s1_server_mode = FALSE;
373 int output_coord_sys = 1;
374 int EQ = 1;
376 if (debug == TRUE)
377 printf("coord_to_string: lon,lat %-.8g %-.8g \n", coord[0], coord[1]);
379 if ((sex == TRUE) && (s1_server_mode == FALSE)) /* server gets decimal */
380 {
381 if (coord[0] < 0.)
382 sign = '-';
383 else
384 sign = ' ';
386 /* Convert longitude */
387 if (output_coord_sys == EQ)
388 {
389 c_hours = coord[0] / 15; /* use hours, min, sec */
390 if (c_hours < 0)
391 c_hours += 24; /* no negative hours */
392 c_hours += .005 / 3600; /* round up one-half of .01 sec */
393 }
394 else
395 {
396 c_hours = fabs(coord[0]); /* use deg, min, sec */
397 c_hours += .05 / 3600; /* round up one-half .1 sec */
398 }
400 hours = c_hours;
401 min = (c_hours - hours) * 60.;
402 dsec = (c_hours - hours - min / 60.) * 3600.;
404 if (output_coord_sys == EQ)
405 {
406 /* remove the rounding (f format does it) */
407 dsec -= .005;
408 if (dsec < 0.0)
409 dsec = 0.0;
410 sprintf(s_coord[0], "%dh%02dm%05.2fs", hours, min, dsec);
411 }
412 else
413 {
414 /* remove the rounding (f format does it) */
415 dsec -= .05;
416 if (dsec < 0.0)
417 dsec = 0.0;
418 sprintf(s_coord[0], "%c%dd%02dm%04.1fs", sign, hours, min, dsec);
419 }
423 /* Convert latitude */
424 if (coord[1] < 0.)
425 sign = '-';
426 else
427 sign = ' ';
429 c_degs = fabs(coord[1]);/* Convert latitude */
430 c_degs += .05 / 3600; /* round up one-half of .1 sec */
432 degs = c_degs;
433 min = (c_degs - degs) * 60.;
434 dsec = (c_degs - degs - min / 60.) * 3600.;
436 /* remove the rounding (f format does it) */
437 dsec -= .05;
438 if (dsec < 0.0)
439 dsec = 0.0;
440 sprintf(s_coord[1], "%c%dd%02dm%04.1fs", sign, degs, min, dsec);
441 }
443 else /* if sex == FALSE */
444 {
445 sprintf(s_coord[0], "%f", coord[0]);
446 sprintf(s_coord[1], "%f", coord[1]);
447 }
448 }