src/coordinates.c

   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 */
  28
  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"
  37
  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 };
  49
  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);
  55
  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     }
  62
  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);
  73
  74     return(EXIT_SUCCESS);
  75 }
  76
  77 /*********************************************************************/
  78 /*                                                                   */
  79 /* EC_to_EQ                                                          */
  80 /* converts ecliptic coordinates to equatorial (both in degrees)     */
  81 /*                                                                   */
  82 /*********************************************************************/
  83
  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;
  91
  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};
  95
  96     if (debug == TRUE)
  97         printf("EC_to_EQ: in %-g %-g\n", coord[0], coord[1]);
  98
  99     lon = coord[0] * dtr;
 100     lat = coord[1] * dtr;
 101
 102     x[0] = cos(lon) * cos(lat);
 103     x[1] = sin(lon) * cos(lat);
 104     x[2] = sin(lat);
 105
 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     }
 114
 115     if (xp[2] > 1.)
 116         xp[2] = 1.;
 117
 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);
 125
 126     if (ra < 0.)
 127         ra += 2.* M_PI;
 128
 129     coord[0] = ra / dtr;
 130     coord[1] = dec / dtr;
 131
 132     if (debug == TRUE)
 133         printf("EC_to_EQ: out %-g %-g\n", coord[0], coord[1]);
 134
 135 }
 136
 137
 138
 139 /*********************************************************************/
 140 /*                                                                   */
 141 /* GA_to_EQ                                                          */
 142 /* converts galactic coordinates to equatorial (both in degrees)     */
 143 /*                                                                   */
 144 /*********************************************************************/
 145
 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];
 152
 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};
 156
 157     if (debug == TRUE)
 158         printf("GA_to_EQ: in %-g %-g\n", coord[0], coord[1]);
 159
 160     lon = coord[0] * dtr;
 161     lat = coord[1] * dtr;
 162
 163     x[0] = cos(lon) * cos(lat);
 164     x[1] = sin(lon) * cos(lat);
 165     x[2] = sin(lat);
 166
 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     }
 175
 176     if (xp[2] > 1.)
 177         xp[2] = 1.;
 178
 179     if (xp[2] < -1.)
 180         xp[2] = -1.;
 181
 182     dec = asin(xp[2]);
 183     ra = atan2(xp[1], xp[0]);
 184
 185     if (ra < 0.)
 186         ra += 2.* M_PI;
 187
 188     coord[0] = ra * rtd;
 189     coord[1] = dec * rtd;
 190
 191     if (debug == TRUE)
 192         printf("GA_to_EQ: out %-g %-g\n", coord[0], coord[1]);
 193 }
 194
 195
 196 /*********************************************************************/
 197 /*                                                                   */
 198 /* EQ_to_EC                                                          */
 199 /* converts equatorial coordinates to ecliptic (both in degrees)     */
 200 /*                                                                   */
 201 /*********************************************************************/
 202
 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];
 209
 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};
 213
 214     if (debug == TRUE)
 215         printf("EQ_to_EC: in %-g %-g\n", coord[0], coord[1]);
 216
 217     ra = coord[0] * dtr;
 218     dec = coord[1] * dtr;
 219
 220     x[0] = cos(ra) * cos(dec);
 221     x[1] = sin(ra) * cos(dec);
 222     x[2] = sin(dec);
 223
 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     }
 232
 233     if (xp[2] > 1.)
 234         xp[2] = 1.;
 235
 236     if (xp[2] < -1.)
 237         xp[2] = -1.;
 238
 239     lat = asin(xp[2]);
 240     lon = atan2(xp[1], xp[0]);
 241
 242     if (lon < 0.)
 243         lon += 2.* M_PI;
 244
 245     coord[0] = lon * rtd;
 246     coord[1] = lat * rtd;
 247
 248     if (debug == TRUE)
 249         printf("EQ_to_EC: out %-g %-g\n", coord[0], coord[1]);
 250 }
 251
 252
 253 /*********************************************************************/
 254 /*                                                                   */
 255 /* GA_to_EC                                                          */
 256 /* converts galactic coordinates to ecliptic (both in degrees)       */
 257 /*                                                                   */
 258 /*********************************************************************/
 259
 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]);
 266
 267     GA_to_EQ(coord);
 268     EQ_to_EC(coord);
 269
 270     if (debug == TRUE)
 271         printf("GA_to_EC: out %-g %-g\n", coord[0], coord[1]);
 272 }
 273
 274
 275
 276 /*********************************************************************/
 277 /*                                                                   */
 278 /* EQ_to_GA                                                          */
 279 /* converts equatorial coordinates to galactic (both in degrees)     */
 280 /*                                                                   */
 281 /*********************************************************************/
 282
 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];
 289
 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};
 293
 294     ra = coord[0] * dtr;
 295     dec = coord[1] * dtr;
 296
 297     x[0] = cos(ra) * cos(dec);
 298     x[1] = sin(ra) * cos(dec);
 299     x[2] = sin(dec);
 300
 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     }
 309
 310     if (xp[2] > 1.)
 311         xp[2] = 1.;
 312
 313     if (xp[2] < -1.)
 314         xp[2] = -1.;
 315
 316     lat = asin(xp[2]);
 317     lon = atan2(xp[1], xp[0]);
 318
 319     if (lon < 0.)
 320         lon += 2.* M_PI;
 321
 322     coord[0] = lon * rtd;
 323     coord[1] = lat * rtd;
 324
 325 }
 326
 327
 328
 329 /*********************************************************************/
 330 /*                                                                   */
 331 /* EC_to_GA                                                          */
 332 /* converts ecliptic coordinates to galactic (both in degrees)       */
 333 /*                                                                   */
 334 /*********************************************************************/
 335
 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]);
 342
 343     EC_to_EQ(coord);
 344     EQ_to_GA(coord);
 345
 346     if (debug == TRUE)
 347         printf("EC_to_GA: in %-g %-g\n", coord[0], coord[1]);
 348 }
 349
 350
 351
 352
 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 /**********************************************************************/
 361
 362
 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;
 375
 376     if (debug == TRUE)
 377         printf("coord_to_string: lon,lat %-.8g %-.8g \n", coord[0], coord[1]);
 378
 379     if ((sex == TRUE) && (s1_server_mode == FALSE))  /* server gets decimal */
 380     {
 381         if (coord[0] < 0.)
 382             sign = '-';
 383         else
 384             sign = ' ';
 385
 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         }
 399
 400         hours = c_hours;
 401         min = (c_hours - hours) * 60.;
 402         dsec = (c_hours - hours - min / 60.) * 3600.;
 403
 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         }
 420
 421
 422
 423         /* Convert latitude */
 424         if (coord[1] < 0.)
 425             sign = '-';
 426         else
 427             sign = ' ';
 428
 429         c_degs = fabs(coord[1]);/* Convert latitude     */
 430         c_degs += .05 / 3600;   /* round up one-half of .1 sec */
 431
 432         degs = c_degs;
 433         min = (c_degs - degs) * 60.;
 434         dsec = (c_degs - degs - min / 60.) * 3600.;
 435
 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     }
 442
 443     else                        /* if sex == FALSE  */
 444     {
 445         sprintf(s_coord[0], "%f", coord[0]);
 446         sprintf(s_coord[1], "%f", coord[1]);
 447     }
 448 }