src/sky/projection.c

   1 /*
   2  * Copyright (C) 2008 Liam Girdwood
   3  *
   4  * This program is free software; you can redistribute it and/or modify
   5  * it under the terms of the GNU General Public License as published by
   6  * the Free Software Foundation; either version 2 of the License, or
   7  * (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software
  16  * Foundation, Inc., 59 Temple Place, Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include <strings.h>
  21 #include <math.h>
  22 #include "projection.h"
  23 #include "astro_object.h"
  24
  25 #define D2R  (1.7453292519943295769e-2)  /* deg->radian */
  26 #define R2D  (5.7295779513082320877e1)   /* radian->deg */
  27
  28 #define PROJ_MAX_RA             360.0
  29 #define PROJ_MIN_RA             0.0
  30
  31 /* Max/Min Declination != 90 or -90 to speed up rendering */
  32 #define PROJ_MAX_DEC            89.5
  33 #define PROJ_MIN_DEC            -89.5
  34
  35 #define PROJ_MIN_FOV            0.5
  36 #define PROJ_MAX_FOV            270.0
  37
  38 #define PROJ_MAG_BASE           4.2
  39 #define PROJ_MAG_FAINT          16
  40 #define PROJ_MAG_BRIGHT         -2
  41
  42 #define PROJ_MAG_CONST          0.005
  43 #define PROJ_MAG_FOV_DIFF       110
  44
  45 struct _mag_const {
  46         gfloat limit;
  47         gfloat coeff;
  48 };
  49
  50 static struct _mag_const mag_const[] = {
  51         {180.0, 2.5 * PROJ_MAG_CONST},
  52         {160.0, 2.6 * PROJ_MAG_CONST},
  53         {140.0, 2.7 * PROJ_MAG_CONST},
  54         {130.0, 2.7 * PROJ_MAG_CONST},
  55         {120.0, 3.0 * PROJ_MAG_CONST},
  56         {110.0, 3.3 * PROJ_MAG_CONST},
  57         {100.0, 3.7 * PROJ_MAG_CONST},
  58         {90.0, 3.9 * PROJ_MAG_CONST},
  59         {80.0, 4.1 * PROJ_MAG_CONST},
  60         {70.0, 4.3 * PROJ_MAG_CONST},
  61         {60.0, 4.5 * PROJ_MAG_CONST},
  62         {50.0, 4.7 * PROJ_MAG_CONST},
  63         {40.0, 5.0 * PROJ_MAG_CONST},
  64         {35.0, 5.2 * PROJ_MAG_CONST},
  65         {30.0, 5.4 * PROJ_MAG_CONST},
  66         {25.0, 5.6 * PROJ_MAG_CONST},
  67         {20.0, 5.8 * PROJ_MAG_CONST},
  68         {16.0, 6.0 * PROJ_MAG_CONST},
  69         {12.0, 6.5 * PROJ_MAG_CONST},
  70         {10.0, 7.0 * PROJ_MAG_CONST},
  71         {8.0, 7.5 * PROJ_MAG_CONST},
  72         {5.0, 7.8 * PROJ_MAG_CONST},
  73         {4.0, 8.0 * PROJ_MAG_CONST},
  74         {3.0, 8.5 * PROJ_MAG_CONST},
  75         {2.0, 9.0 * PROJ_MAG_CONST},
  76         {1.5, 10.0 * PROJ_MAG_CONST},
  77         {1.0, 11.0 * PROJ_MAG_CONST},
  78 };
  79
  80 gint projection_set_flip(struct projection *proj, enum projection_flip flip)
  81 {
  82         switch (flip) {
  83         case PF_NONE:
  84                 proj->trans = &transform_flip_none;
  85                 break;
  86         case PF_TB:
  87                 proj->trans = &transform_flip_tb;
  88                 break;
  89         case PF_LR:
  90                 proj->trans = &transform_flip_lr;
  91                 break;
  92         case PF_LR_TB:
  93                 proj->trans = &transform_flip_lr_tb;
  94                 break;
  95         }
  96
  97         proj->flip = flip;
  98         proj->clip_mag_faint = PROJ_MAG_FAINT;
  99         proj->clip_mag_bright = PROJ_MAG_BRIGHT;
 100         return 0;
 101 }
 102
 103 /*! \fn void void proj_calc_fov_bounds(struct projection* projection)
 104 *
 105 * Calculate the boundaries, magnitude level and field of view for the
 106 * virtual sky projection.
 107 */
 108 void projection_fov_bounds(struct projection* proj)
 109 {
 110         gfloat mag_coeff = PROJ_MAG_CONST;
 111         gint i = 0;
 112
 113         proj->sky_mid_width = proj->sky_width / 2.0;
 114         proj->sky_mid_height = proj->sky_height / 2.0;
 115
 116         if (proj->sky_width > proj->sky_height) {
 117                 proj->pixels_per_degree = proj->sky_width / proj->fov;
 118                 proj->pixels_per_radian = proj->pixels_per_degree / D2R;
 119         } else {
 120                 proj->pixels_per_degree = proj->sky_height / proj->fov;
 121                 proj->pixels_per_radian = proj->pixels_per_degree / D2R;
 122         }
 123         while (mag_const[i].limit > proj->fov && i < nsize(mag_const)) {
 124                 mag_coeff = mag_const[i].coeff;
 125                 i++;
 126         }
 127         proj->clip_mag_faint =
 128                 PROJ_MAG_BASE + (PROJ_MAG_FOV_DIFF - proj->fov) * mag_coeff;
 129
 130         proj->centre_ra_rad = proj->centre_ra * D2R;
 131         proj->centre_dec_rad = proj->centre_dec * D2R;
 132         proj->centre_ra_rad_sin = sin(proj->centre_ra_rad);
 133         proj->centre_dec_rad_sin = sin(proj->centre_dec_rad);
 134         proj->centre_ra_rad_cos = cos(proj->centre_ra_rad);
 135         proj->centre_dec_rad_cos = cos(proj->centre_dec_rad);
 136 }
 137
 138 void projection_check_bounds (struct projection* proj)
 139 {
 140         if (proj->centre_ra > PROJ_MAX_RA)
 141                 proj->centre_ra -= PROJ_MAX_RA;
 142         if (proj->centre_ra < PROJ_MIN_RA)
 143                 proj->centre_ra += PROJ_MAX_RA;
 144         if (proj->centre_dec > PROJ_MAX_DEC)
 145                 proj->centre_dec = PROJ_MAX_DEC;
 146         if (proj->centre_dec < PROJ_MIN_DEC)
 147                 proj->centre_dec = PROJ_MIN_DEC;
 148         if (proj->fov > PROJ_MAX_FOV)
 149                 proj->fov = PROJ_MAX_FOV;
 150         if (proj->fov < PROJ_MIN_FOV)
 151                 proj->fov = PROJ_MIN_FOV;
 152 }
 153
 154 static void sky_to_cairo(struct projection *proj, struct render_object *robject)
 155 {
 156         gdouble mid_x, mid_y, k;
 157         gdouble centre_ra = proj->centre_ra_rad;
 158         gdouble centre_dec_cos = proj->centre_dec_rad_cos;
 159         gdouble centre_dec_sin = proj->centre_dec_rad_sin;
 160         gint i;
 161
 162         mid_x = proj->sky_mid_width;
 163         mid_y = proj->sky_mid_height;
 164
 165         for (i = 0; i < robject->num_coords; i++) {
 166
 167                 gdouble ra = robject->coord[i].posn->ra * D2R;
 168                 gdouble dec = robject->coord[i].posn->dec * D2R;
 169                 gdouble sin_dec = sin(dec);
 170                 gdouble cos_dec = cos(dec);
 171                 gdouble cos_ra_delta = cos(ra - centre_ra);
 172
 173                 k = 2.0 / (1 + centre_dec_sin * sin_dec +
 174                         centre_dec_cos * cos_dec * cos_ra_delta);
 175
 176                 robject->coord[i].x =
 177                         k * (cos_dec * sin(ra - centre_ra));
 178                 robject->coord[i].y =
 179                         k * (centre_dec_cos * sin_dec - centre_dec_sin *
 180                                 cos_dec * cos_ra_delta);
 181
 182                 robject->coord[i].x =
 183                         mid_x - robject->coord[i].x * proj->pixels_per_radian;
 184                 robject->coord[i].y =
 185                         mid_y - robject->coord[i].y * proj->pixels_per_radian;
 186         }
 187 }
 188
 189 static void cairo_to_sky(struct projection *proj, struct render_object *robject)
 190 {
 191         gdouble mid_x, mid_y, p, c;
 192         gdouble centre_ra = proj->centre_ra_rad;
 193         gdouble centre_dec = proj->centre_dec_rad;
 194
 195         mid_x = proj->sky_mid_width;
 196         mid_y = proj->sky_mid_height;
 197
 198         robject->coord[0].x =
 199                 (mid_x - robject->coord[0].x) / proj->pixels_per_radian;
 200         robject->coord[0].y =
 201                 (mid_y - robject->coord[0].y) / proj->pixels_per_radian;
 202
 203         p = sqrt(robject->coord[0].x * robject->coord[0].x +
 204                         robject->coord[0].y * robject->coord[0].y);
 205         c = 2.0 * atan2(p, 2.0);
 206
 207         robject->coord[0].posn->dec = asin(cos(c) * sin(centre_dec) +
 208                 ((robject->coord[0].y * sin(c) * cos (centre_dec)) / p));
 209
 210         robject->coord[0].posn->ra = centre_ra +
 211                 atan2(robject->coord[0].x * sin(c),
 212                 p * cos(centre_dec) * cos(c) -
 213                 robject->coord[0].y * sin(centre_dec) * sin(c));
 214
 215         robject->coord[0].posn->dec = robject->coord[0].posn->dec * R2D;
 216         robject->coord[0].posn->ra = robject->coord[0].posn->ra * R2D;
 217 }
 218
 219 static void sky_to_cairo_tb(struct projection *proj, struct render_object *robject)
 220 {
 221         gdouble mid_x, mid_y, k;
 222         gdouble centre_ra = proj->centre_ra_rad;
 223         gdouble centre_dec = proj->centre_dec_rad;
 224         gint i;
 225
 226         mid_x = proj->sky_mid_width;
 227         mid_y = proj->sky_mid_height;
 228
 229         for (i = 0; i < robject->num_coords; i++) {
 230
 231                 gdouble ra = robject->coord[i].posn->ra * D2R;
 232                 gdouble dec = robject->coord[i].posn->dec * D2R;
 233
 234                 k = 2.0 / (1 + sin(centre_dec) * sin(dec) +
 235                         cos(centre_dec) * cos(dec) * cos(ra - centre_ra));
 236
 237                 robject->coord[i].x =
 238                         k * (cos(dec) * sin(ra - centre_ra));
 239                 robject->coord[i].y =
 240                         k * (cos(centre_dec) * sin(dec) - sin(centre_dec) *
 241                                 cos(dec) * cos(ra - centre_ra));
 242
 243                 robject->coord[i].x =
 244                         mid_x - robject->coord[i].x * proj->pixels_per_radian;
 245                 robject->coord[i].y =
 246                         mid_y + robject->coord[i].y * proj->pixels_per_radian;
 247         }
 248 }
 249
 250 static void cairo_to_sky_tb(struct projection *proj, struct render_object *robject)
 251 {
 252         gdouble mid_x, mid_y, p, c;
 253         gdouble centre_ra = proj->centre_ra_rad;
 254         gdouble centre_dec = proj->centre_dec_rad;
 255
 256         mid_x = proj->sky_mid_width;
 257         mid_y = proj->sky_mid_height;
 258
 259         robject->coord[0].x =
 260                 (mid_x + robject->coord[0].x) / proj->pixels_per_radian;
 261         robject->coord[0].y =
 262                 (mid_y - robject->coord[0].y) / proj->pixels_per_radian;
 263
 264         p = sqrt(robject->coord[0].x * robject->coord[0].x +
 265                         robject->coord[0].y * robject->coord[0].y);
 266         c = 2.0 * atan2(p, 2.0);
 267
 268         robject->coord[0].posn->dec = asin(cos(c) * sin(centre_dec) +
 269                 ((robject->coord[0].y * sin(c) * cos (centre_dec)) / p));
 270
 271         robject->coord[0].posn->ra = centre_ra +
 272                 atan2(robject->coord[0].x * sin(c),
 273                 p * cos(centre_dec) * cos(c) -
 274                 robject->coord[0].y * sin(centre_dec) * sin(c));
 275
 276         robject->coord[0].posn->dec = robject->coord[0].posn->dec * R2D;
 277         robject->coord[0].posn->ra = robject->coord[0].posn->ra * R2D;
 278 }
 279
 280 static void sky_to_cairo_lr(struct projection *proj, struct render_object *robject)
 281 {
 282         gdouble mid_x, mid_y, k;
 283         gdouble centre_ra = proj->centre_ra_rad;
 284         gdouble centre_dec = proj->centre_dec_rad;
 285         gint i;
 286
 287         mid_x = proj->sky_mid_width;
 288         mid_y = proj->sky_mid_height;
 289
 290         for (i = 0; i < robject->num_coords; i++) {
 291
 292                 gdouble ra = robject->coord[i].posn->ra * D2R;
 293                 gdouble dec = robject->coord[i].posn->dec * D2R;
 294
 295                 k = 2.0 / (1 + sin(centre_dec) * sin(dec) +
 296                         cos(centre_dec) * cos(dec) * cos(ra - centre_ra));
 297
 298                 robject->coord[i].x =
 299                         k * (cos(dec) * sin(ra - centre_ra));
 300                 robject->coord[i].y =
 301                         k * (cos(centre_dec) * sin(dec) - sin(centre_dec) *
 302                                 cos(dec) * cos(ra - centre_ra));
 303
 304                 robject->coord[i].x =
 305                         mid_x + robject->coord[i].x * proj->pixels_per_radian;
 306                 robject->coord[i].y =
 307                         mid_y - robject->coord[i].y * proj->pixels_per_radian;
 308         }
 309 }
 310
 311 static void cairo_to_sky_lr(struct projection *proj, struct render_object *robject)
 312 {
 313         gdouble mid_x, mid_y, p, c;
 314         gdouble centre_ra = proj->centre_ra_rad;
 315         gdouble centre_dec = proj->centre_dec_rad;
 316
 317         mid_x = proj->sky_mid_width;
 318         mid_y = proj->sky_mid_height;
 319
 320         robject->coord[0].x =
 321                 (mid_x - robject->coord[0].x) / proj->pixels_per_radian;
 322         robject->coord[0].y =
 323                 (mid_y + robject->coord[0].y) / proj->pixels_per_radian;
 324
 325         p = sqrt(robject->coord[0].x * robject->coord[0].x +
 326                         robject->coord[0].y * robject->coord[0].y);
 327         c = 2.0 * atan2(p, 2.0);
 328
 329         robject->coord[0].posn->dec = asin(cos(c) * sin(centre_dec) +
 330                 ((robject->coord[0].y * sin(c) * cos (centre_dec)) / p));
 331
 332         robject->coord[0].posn->ra = centre_ra +
 333                 atan2(robject->coord[0].x * sin(c),
 334                 p * cos(centre_dec) * cos(c) -
 335                 robject->coord[0].y * sin(centre_dec) * sin(c));
 336
 337         robject->coord[0].posn->dec = robject->coord[0].posn->dec * R2D;
 338         robject->coord[0].posn->ra = robject->coord[0].posn->ra * R2D;
 339 }
 340
 341 static void sky_to_cairo_lr_tb(struct projection *proj, struct render_object *robject)
 342 {
 343         gdouble mid_x, mid_y, k;
 344         gdouble centre_ra = proj->centre_ra_rad;
 345         gdouble centre_dec = proj->centre_dec_rad;
 346         gint i;
 347
 348         mid_x = proj->sky_mid_width;
 349         mid_y = proj->sky_mid_height;
 350
 351         for (i = 0; i < robject->num_coords; i++) {
 352
 353                 gdouble ra = robject->coord[i].posn->ra * D2R;
 354                 gdouble dec = robject->coord[i].posn->dec * D2R;
 355
 356                 k = 2.0 / (1 + sin(centre_dec) * sin(dec) +
 357                         cos(centre_dec) * cos(dec) * cos(ra - centre_ra));
 358
 359                 robject->coord[i].x =
 360                         k * (cos(dec) * sin(ra - centre_ra));
 361                 robject->coord[i].y =
 362                         k * (cos(centre_dec) * sin(dec) - sin(centre_dec) *
 363                                 cos(dec) * cos(ra - centre_ra));
 364
 365                 robject->coord[i].x =
 366                         mid_x + robject->coord[i].x * proj->pixels_per_radian;
 367                 robject->coord[i].y =
 368                         mid_y + robject->coord[i].y * proj->pixels_per_radian;
 369         }
 370 }
 371
 372 static void cairo_to_sky_lr_tb(struct projection *proj, struct render_object *robject)
 373 {
 374         gdouble mid_x, mid_y, p, c;
 375         gdouble centre_ra = proj->centre_ra_rad;
 376         gdouble centre_dec = proj->centre_dec_rad;
 377
 378         mid_x = proj->sky_mid_width;
 379         mid_y = proj->sky_mid_height;
 380
 381         robject->coord[0].x =
 382                 (mid_x + robject->coord[0].x) / proj->pixels_per_radian;
 383         robject->coord[0].y =
 384                 (mid_y + robject->coord[0].y) / proj->pixels_per_radian;
 385
 386         p = sqrt(robject->coord[0].x * robject->coord[0].x +
 387                         robject->coord[0].y * robject->coord[0].y);
 388         c = 2.0 * atan2(p, 2.0);
 389
 390         robject->coord[0].posn->dec = asin(cos(c) * sin(centre_dec) +
 391                 ((robject->coord[0].y * sin(c) * cos (centre_dec)) / p));
 392
 393         robject->coord[0].posn->ra = centre_ra +
 394                 atan2(robject->coord[0].x * sin(c),
 395                 p * cos(centre_dec) * cos(c) -
 396                 robject->coord[0].y * sin(centre_dec) * sin(c));
 397
 398         robject->coord[0].posn->dec = robject->coord[0].posn->dec * R2D;
 399         robject->coord[0].posn->ra = robject->coord[0].posn->ra * R2D;
 400 }
 401
 402
 403 struct transform transform_flip_none = {
 404         .sky_to_proj = sky_to_cairo,
 405         .proj_to_sky = cairo_to_sky,
 406 };
 407
 408 struct transform transform_flip_tb = {
 409         .sky_to_proj = sky_to_cairo_tb,
 410         .proj_to_sky = cairo_to_sky_tb,
 411 };
 412
 413 struct transform transform_flip_lr = {
 414         .sky_to_proj = sky_to_cairo_lr,
 415         .proj_to_sky = cairo_to_sky_lr,
 416 };
 417
 418 struct transform transform_flip_lr_tb = {
 419         .sky_to_proj = sky_to_cairo_lr_tb,
 420         .proj_to_sky = cairo_to_sky_lr_tb,
 421 };