src/sky/grid.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  * Fast sky grid render engine. Mostly working, although probably needs
  20  * checking by a maths trig guru to further optimise (and clip RA line at 80).
  21  *
  22  * TODO: Clip RA grid lines at 80 and -80 on DEC. This needs cairo_arc_to() or
  23  * someone who knows another solution.
  24  */
  25
  26 #define _GNU_SOURCE /* for NAN and INF */
  27
  28 #include <math.h>
  29 #include <stdio.h>
  30 #include "astro_object.h"
  31
  32 #include "grid.h"
  33
  34 #define D2R  (1.7453292519943295769e-2)  /* deg->radian */
  35 #define R2D  (5.7295779513082320877e1)   /* radian->deg */
  36
  37 #define ARC_MAGIC 1000.0 /* need better method to determine line from arc */
  38
  39 #define GRID_NIGHT_ALPHA        0.35
  40 #define GRID_DAY_ALPHA          0.55
  41
  42
  43 /* RA grid in arcsecs */
  44 static const gdouble grid_hms_ra[] = {
  45         15.0 * 60.0 * 60.0,
  46         15.0 * 60.0 * 10.0,
  47         15.0 * 60.0 * 5.0,
  48         15.0 * 60.0,
  49         15.0 * 10.0,
  50         15.0 * 5.0,
  51         15.0,
  52         7.5,
  53         1.5,
  54 };
  55
  56 /* DEC drid in arcsecs */
  57 static const gdouble grid_dms_dec[] = {
  58         3600.0 * 10.0,
  59         3600.0 * 5.0,
  60         3600.0,
  61         1200.0,
  62         600.0,
  63         300.0,
  64         60.0,
  65         10.0,
  66         5.0,
  67         1.0,
  68 };
  69
  70 /* get grid RA step size for projection in arcsecs */
  71 static gdouble get_ra_step_delta(struct projection *proj)
  72 {
  73         gint i;
  74
  75         for (i = 0; i < nsize(grid_hms_ra); i++) {
  76                 if (proj->fov / 2.0 > grid_hms_ra[i] / 3600.0)
  77                         return grid_hms_ra[i] / 3600.0;
  78         }
  79         return grid_hms_ra[i] / 3600.0;
  80 }
  81
  82 /* get grid DEC step size for projection in arcsecs */
  83 static gdouble get_dec_step_delta(struct projection *proj)
  84 {
  85         gint i;
  86
  87         for (i = 0; i < nsize(grid_dms_dec); i++) {
  88                 if (proj->fov / 2.0 > grid_dms_dec[i] / 3600.0)
  89                         return grid_dms_dec[i] / 3600.0;
  90         }
  91         return grid_dms_dec[i] / 3600.0;
  92 }
  93
  94 /* clip DEC grid to current projection FOV */
  95 static void clip_dec_grid(struct projection *proj, gdouble step_size,
  96         gdouble *south, gdouble *north)
  97 {
  98         *south = proj->centre_dec - proj->fov / 2.0;
  99         if (*south <= -90.0)
 100                 *south = -90.0 + step_size;
 101         else
 102                 *south = floor(*south / step_size) * step_size;
 103
 104         *north = proj->centre_dec + proj->fov / 2.0;
 105         if (*north >= 90.0)
 106                 *north = 90.0 - step_size;
 107         else
 108                 *north = floor(*north / step_size) * step_size;
 109 }
 110
 111 /* clip RA grid to current projection FOV */
 112 static void clip_ra_grid(struct projection *proj, gdouble step_size,
 113         gdouble *start, gdouble *end)
 114 {
 115         struct render_object robject;
 116         struct ln_equ_posn pos[4];
 117         gdouble start_, end_;
 118         gint i;
 119
 120         /* get RA at each projection corner */
 121         robject.coord[0].x = 0;
 122         robject.coord[0].y = 0;
 123         robject.coord[0].posn = &pos[0];
 124         proj->trans->proj_to_sky(proj, &robject);
 125
 126         robject.coord[0].x = proj->sky_width;
 127         robject.coord[0].y = 0;
 128         robject.coord[0].posn = &pos[1];
 129         proj->trans->proj_to_sky(proj, &robject);
 130
 131         robject.coord[0].x = 0;
 132         robject.coord[0].y = proj->sky_height;
 133         robject.coord[0].posn = &pos[2];
 134         proj->trans->proj_to_sky(proj, &robject);
 135
 136         robject.coord[0].x = proj->sky_width;
 137         robject.coord[0].y = proj->sky_height;
 138         robject.coord[0].posn = &pos[3];
 139         proj->trans->proj_to_sky(proj, &robject);
 140
 141         /* get highest and lowest RA from corners */
 142         start_ = 360.0;
 143         end_ = 0.0;
 144         for (i = 0; i < 3; i++) {
 145                 if (start_ > pos[i].ra)
 146                         start_ = pos[i].ra;
 147                 if (end_ < pos[i].ra)
 148                         end_ = pos[i].ra;
 149         }
 150
 151         *start = start_;
 152         *start = floor(*start / step_size) * step_size;
 153
 154         *end = end_;
 155         *end = floor(*end / step_size) * step_size;
 156 }
 157
 158 /* should we render a line of this radius as an arc or a straight line */
 159 static inline gint is_arc_line(struct projection *proj, gdouble radius)
 160 {
 161         /* magic number - must refine */
 162         if (radius > proj->sky_width * ARC_MAGIC ||
 163                 radius < proj->sky_width * -ARC_MAGIC) {
 164                 //printf("%s rad %f\n", __func__, radius);
 165                 return 1;
 166         }
 167         return 0;
 168 }
 169
 170 /* is this arc visible in projection */
 171 static gint is_arc_visible(struct projection *proj, gdouble x, gdouble y,
 172         gdouble radius)
 173 {
 174         gdouble x1, y1, x2, y2;
 175
 176         /* sanity check */
 177         if (radius == FP_NAN
 178                 || radius == FP_INFINITE
 179                 || radius == -FP_INFINITE) {
 180                 return 0;
 181         }
 182
 183         /* create clip area sky size + radius */
 184         x1 = -radius;
 185         y1 = -radius;
 186         x2 = proj->sky_width + radius;
 187         y2 = proj->sky_height + radius;
 188
 189         /* are we inside ? */
 190         if (x >= x1 && x <= x2 && y >= y1 && y <= y2)
 191                 return 1;
 192
 193         return 0;
 194 }
 195
 196 static inline void do_ra_line(struct render_object *robject,
 197         struct projection *proj, gdouble x, gdouble y, gdouble my)
 198 {
 199         gdouble x1, y1, x2, y2;
 200
 201         x1 = x - my * y;
 202         y1 = 0;
 203         x2 = my * (proj->sky_height - y) + x;
 204         y2 = proj->sky_height;
 205
 206         cairo_move_to(robject->cr, x1, y1);
 207         cairo_line_to(robject->cr, x2, y2);
 208 }
 209
 210 /* draw a straight RA grid line */
 211 static inline void draw_ra_line(struct render_object *robject,
 212         struct projection *proj, gdouble step)
 213 {
 214         struct ln_equ_posn pos1, pos2;
 215         gdouble my;
 216
 217         /* calc based on small RA line */
 218         pos1.dec = proj->centre_dec + 1.0;
 219         pos2.dec = proj->centre_dec -1.0;
 220         robject->num_coords = 2;
 221
 222         pos1.ra = step;
 223         pos2.ra = step + 180.0;
 224
 225         robject->coord[0].posn = &pos1;
 226         robject->coord[1].posn = &pos2;
 227
 228         proj->trans->sky_to_proj(proj, robject);
 229
 230         /* calculate line gradient */
 231         my = (robject->coord[1].x - robject->coord[0].x) /
 232                 (robject->coord[1].y - robject->coord[0].y);
 233
 234         do_ra_line(robject, proj,
 235                 robject->coord[0].x, robject->coord[0].y, my);
 236 }
 237
 238 static inline void do_dec_line(struct render_object *robject,
 239         struct projection *proj, gdouble x, gdouble y, gdouble mx)
 240 {
 241         gdouble x1, y1, x2, y2;
 242
 243         x1 = 0;
 244         y1 = y -mx * y;
 245         x2 = proj->sky_width;
 246         y2 = mx * (proj->sky_width - x) + y;
 247
 248         cairo_move_to(robject->cr, x1, y1);
 249         cairo_line_to(robject->cr, x2, y2);
 250 }
 251
 252 /* draw a straight DEC grid line */
 253 static inline void draw_dec_line(struct render_object *robject,
 254         struct projection *proj, gdouble step)
 255 {
 256         struct ln_equ_posn pos1, pos2;
 257         gdouble mx;
 258
 259         /* calc based on small DEC line */
 260         pos1.ra = proj->centre_ra + 1.0;
 261         pos2.ra = proj->centre_ra -1.0;
 262         robject->num_coords = 2;
 263
 264         pos1.dec = step;
 265         pos2.dec = step;
 266
 267         robject->coord[0].posn = &pos1;
 268         robject->coord[1].posn = &pos2;
 269
 270         proj->trans->sky_to_proj(proj, robject);
 271
 272         /* calculate line gradient */
 273         mx = (robject->coord[1].y - robject->coord[0].y) /
 274                 (robject->coord[1].x - robject->coord[0].x);
 275
 276         do_dec_line(robject, proj,
 277                 robject->coord[0].x, robject->coord[0].y, mx);
 278 }
 279
 280 /* render all visible RA grid lines */
 281 static void grid_render_ra(struct render_object *robject,
 282         struct projection *proj, gint labels)
 283 {
 284         gdouble centre_x, centre_y, radius, x, y, step,
 285                 step_delta, ra_start, ra_end;
 286         struct ln_equ_posn pos1, pos2;
 287
 288         pos1.dec = -proj->centre_dec;
 289         pos2.dec = -proj->centre_dec;
 290         robject->num_coords = 2;
 291
 292         step_delta = get_ra_step_delta(proj);
 293         clip_ra_grid(proj, step_delta, &ra_start, &ra_end);
 294
 295         for (step = ra_start; step <= ra_end; step += step_delta) {
 296
 297                 pos1.ra = step;
 298                 pos2.ra = step + 180.0;
 299
 300                 robject->coord[0].posn = &pos1;
 301                 robject->coord[1].posn = &pos2;
 302
 303                 proj->trans->sky_to_proj(proj, robject);
 304
 305                 centre_x = robject->coord[0].x +
 306                         ((robject->coord[1].x - robject->coord[0].x) / 2.0);
 307                 centre_y = robject->coord[0].y +
 308                         ((robject->coord[1].y - robject->coord[0].y) / 2.0);
 309                 x = centre_x - robject->coord[0].x;
 310                 y = centre_y - robject->coord[0].y;
 311                 radius = sqrt(x * x + y * y);
 312
 313                 if (!is_arc_visible(proj, centre_x, centre_y, radius) &&
 314                         step != proj->centre_ra)
 315                         continue;
 316
 317 /* need cairo_arc_to() */
 318 #if 0
 319                 if (step == 0.0 || step == 90.0 ||
 320                         step == 180.0 || step == 270.0) {
 321                         struct ln_equ_posn pos1, pos2;
 322
 323                         pos1.dec = 80.0;
 324                         pos2.dec = -80.0;
 325                         pos1.ra = step;
 326                         pos2.ra = step;
 327                         robject->num_coords = 2;
 328
 329                         robject->coord[0].posn = &pos1;
 330                         robject->coord[1].posn = &pos2;
 331                         proj->trans->sky_to_proj(proj, robject);
 332
 333                         cairo_new_sub_path(robject->cr);
 334                         cairo_arc_to(robject->cr,
 335                                 robject->coord[0].x, robject->coord[0].y,
 336                                 robject->coord[1].x, robject->coord[1].y,
 337                                 radius);
 338                 } else {
 339                         if (is_arc_line(proj, centre_x))
 340                                 draw_ra_line(robject, proj, step);
 341                         else {
 342                                 cairo_new_sub_path(robject->cr);
 343                                 cairo_arc(robject->cr, centre_x, centre_y,
 344                                         radius, circle_start, circle_end);
 345                         }
 346                 }
 347 #else
 348                 if (is_arc_line(proj, centre_x))
 349                         draw_ra_line(robject, proj, step);
 350                 else {
 351                         cairo_new_sub_path(robject->cr);
 352                         cairo_arc(robject->cr, centre_x, centre_y,
 353                                         radius, 0, 2 * M_PI);
 354                 }
 355 #endif
 356         }
 357 }
 358
 359 /* render all visible DEC grid lines */
 360 static void grid_render_dec(struct render_object *robject,
 361         struct projection *proj, gint labels)
 362 {
 363         gdouble centre_x, centre_y, radius, x, y, step,
 364                 step_delta, dec_start, dec_end;
 365         struct ln_equ_posn pos1, pos2;
 366
 367         pos1.ra = proj->centre_ra;
 368         pos2.ra = proj->centre_ra + 180.0;
 369         robject->num_coords = 2;
 370
 371         step_delta = get_dec_step_delta(proj);
 372         clip_dec_grid(proj, step_delta, &dec_start, &dec_end);
 373
 374         for (step = dec_start; step <= dec_end; step += step_delta) {
 375
 376                 pos1.dec = step;
 377                 pos2.dec = step;
 378
 379                 robject->coord[0].posn = &pos1;
 380                 robject->coord[1].posn = &pos2;
 381
 382                 proj->trans->sky_to_proj(proj, robject);
 383
 384                 centre_x = robject->coord[0].x +
 385                         ((robject->coord[1].x - robject->coord[0].x) / 2.0);
 386                 centre_y = robject->coord[0].y +
 387                         ((robject->coord[1].y - robject->coord[0].y) / 2.0);
 388                 x = centre_x - robject->coord[0].x;
 389                 y = centre_y - robject->coord[0].y;
 390                 radius = sqrt(x * x + y * y);
 391
 392                 if (!is_arc_visible(proj, centre_x, centre_y, radius) &&
 393                         step != proj->centre_dec)
 394                         continue;
 395
 396                 cairo_new_sub_path(robject->cr);
 397
 398                 if (is_arc_line(proj, centre_x))
 399                         draw_dec_line(robject, proj, step);
 400                 else
 401                         cairo_arc(robject->cr, centre_x, centre_y,
 402                                         radius, 0, 2 * M_PI);
 403         }
 404 }
 405
 406 static void render_dec_labels_at(struct projection *proj,
 407         struct render_object *robject, gdouble dec)
 408 {
 409         gchar text[32];
 410         gdouble step_delta, delta_div = 2.0;
 411         struct ln_equ_posn pos_start, pos_end;
 412         struct ln_dms dms;
 413
 414         step_delta = get_dec_step_delta(proj);
 415
 416         pos_start.ra = proj->centre_ra;
 417         pos_end.ra = proj->centre_ra;
 418         pos_start.dec = pos_end.dec = dec;
 419         robject->num_coords = 2;
 420         robject->coord[0].posn = &pos_start;
 421         robject->coord[1].posn = &pos_end;
 422
 423         /* find start position */
 424         do {
 425                 pos_start.ra -= step_delta;
 426                 proj->trans->sky_to_proj(proj, robject);
 427         } while (projection_is_visible(proj, 0, robject) &&
 428                 proj->centre_ra - pos_start.ra < 90.0);
 429         pos_start.ra += step_delta;
 430
 431         /* binary chop start */
 432         do {
 433                 pos_start.ra -= step_delta / delta_div;
 434                 proj->trans->sky_to_proj(proj, robject);
 435                 if (!projection_is_visible(proj, 0, robject))
 436                         pos_start.ra += step_delta / delta_div;
 437                 delta_div *= 2.0;
 438         } while (delta_div < 128.0);
 439
 440         /* find end position */
 441         do {
 442                 pos_end.ra += step_delta;
 443                 proj->trans->sky_to_proj(proj, robject);
 444         } while (projection_is_visible(proj, 1, robject) &&
 445                 pos_end.ra - proj->centre_ra < 90.0);
 446         pos_end.ra -= step_delta;
 447
 448         /* binary chop end */
 449         delta_div = 2.0;
 450         do {
 451                 pos_end.ra += step_delta / delta_div;
 452                 proj->trans->sky_to_proj(proj, robject);
 453                 if (!projection_is_visible(proj, 1, robject))
 454                         pos_end.ra -= step_delta / delta_div;
 455                 delta_div *= 2.0;
 456         } while (delta_div < 128.0);
 457
 458         ln_deg_to_dms(dec, &dms);
 459         if (dms.neg)
 460                 text[0] = '-';
 461         else
 462                 text[0] = '+';
 463         if (dms.minutes == 0 && dms.seconds == 0.0)
 464                 sprintf(&text[1], "%2.2dº", dms.degrees);
 465         else if (dms.seconds < 0.1)
 466                 sprintf(&text[1], "%2.2dº%2.2dm", dms.degrees, dms.minutes);
 467         else
 468                 sprintf(&text[1], "%2.2dº%2.2dm%2.0f", dms.degrees, dms.minutes,
 469                         dms.seconds);
 470
 471         cairo_move_to(robject->cr,
 472                         robject->coord[0].x - 30.0,
 473                         robject->coord[0].y - 1.0);
 474         cairo_show_text(robject->cr, text);
 475         cairo_move_to(robject->cr,
 476                         robject->coord[1].x + 1.0,
 477                         robject->coord[1].y - 1.0);
 478         cairo_show_text(robject->cr, text);
 479 }
 480
 481 static void render_ra_labels_at(struct projection *proj,
 482         struct render_object *robject, gdouble ra)
 483 {
 484         gchar text[32];
 485         gdouble step_delta, delta_div = 2.0;
 486         struct ln_equ_posn pos_start, pos_end;
 487         struct ln_hms hms;
 488
 489         step_delta = get_ra_step_delta(proj);
 490
 491         pos_start.dec = proj->centre_dec;
 492         pos_end.dec = proj->centre_dec;
 493         pos_start.ra = pos_end.ra = ra;
 494         robject->num_coords = 2;
 495         robject->coord[0].posn = &pos_start;
 496         robject->coord[1].posn = &pos_end;
 497
 498         /* find start position */
 499         do {
 500                 pos_start.dec -= step_delta;
 501                 proj->trans->sky_to_proj(proj, robject);
 502         } while (projection_is_visible(proj, 0, robject) &&
 503                 proj->centre_dec - pos_start.dec < 90.0);
 504         pos_start.dec += step_delta;
 505
 506         /* binary chop start */
 507         do {
 508                 pos_start.dec -= step_delta / delta_div;
 509                 proj->trans->sky_to_proj(proj, robject);
 510                 if (!projection_is_visible(proj, 0, robject))
 511                         pos_start.dec += step_delta / delta_div;
 512                 delta_div *= 2.0;
 513         } while (delta_div < 128.0);
 514
 515         /* find end position */
 516         do {
 517                 pos_end.dec += step_delta;
 518                 proj->trans->sky_to_proj(proj, robject);
 519         } while (projection_is_visible(proj, 1, robject) &&
 520                 pos_end.dec - proj->centre_dec < 90.0);
 521         pos_end.dec -= step_delta;
 522
 523         /* binary chop end */
 524         delta_div = 2.0;
 525         do {
 526                 pos_end.dec += step_delta / delta_div;
 527                 proj->trans->sky_to_proj(proj, robject);
 528                 if (!projection_is_visible(proj, 1, robject))
 529                         pos_end.dec -= step_delta / delta_div;
 530                 delta_div *= 2.0;
 531         } while (delta_div < 128.0);
 532
 533         ln_deg_to_hms(ra, &hms);
 534         if (hms.minutes == 0 && hms.seconds == 0.0)
 535                 sprintf(text, "%2.0dh", hms.hours);
 536         else if (hms.seconds < 0.1)
 537                 sprintf(text, "%2.0dh%2.2dm", hms.hours, hms.minutes);
 538         else
 539                 sprintf(text, "%2.0dh%2.2dm%2.0f", hms.hours, hms.minutes,
 540                         hms.seconds);
 541
 542         cairo_move_to(robject->cr,
 543                         robject->coord[0].x,
 544                         robject->coord[0].y);
 545         cairo_show_text(robject->cr, text);
 546         cairo_move_to(robject->cr,
 547                         robject->coord[1].x,
 548                         robject->coord[1].y + 15.0); // font size
 549         cairo_show_text(robject->cr, text);
 550 }
 551
 552
 553 void grid_render_dec_labels(struct render_object *robject, struct projection *proj)
 554 {
 555         gdouble step, step_delta, dec_start, dec_end;
 556         struct ln_equ_posn pos1;
 557
 558         pos1.ra = proj->centre_ra;
 559         robject->num_coords = 1;
 560
 561         step_delta = get_dec_step_delta(proj);
 562         clip_dec_grid(proj, step_delta, &dec_start, &dec_end);
 563         robject->coord[0].posn = &pos1;
 564
 565         for (step = dec_start; step <= dec_end; step += step_delta) {
 566                 pos1.dec = step;
 567                 proj->trans->sky_to_proj(proj, robject);
 568                 render_dec_labels_at(proj, robject, step);
 569         }
 570 }
 571
 572 void grid_render_ra_labels(struct render_object *robject, struct projection *proj)
 573 {
 574         gdouble step, step_delta, ra_start, ra_end;
 575         struct ln_equ_posn pos1;
 576
 577         pos1.ra = proj->centre_dec;
 578         robject->num_coords = 1;
 579
 580         step_delta = get_ra_step_delta(proj);
 581         clip_ra_grid(proj, step_delta, &ra_start, &ra_end);
 582         robject->coord[0].posn = &pos1;
 583
 584         for (step = ra_start; step <= ra_end; step += step_delta) {
 585                 pos1.ra = step;
 586                 proj->trans->sky_to_proj(proj, robject);
 587                 render_ra_labels_at(proj, robject, step);
 588         }
 589 }
 590
 591 void grid_render(struct render_object *robject, struct projection *proj,
 592                         gint labels)
 593 {
 594         gdouble alpha;
 595
 596         if (robject->flags.night_mode)
 597                 alpha = GRID_NIGHT_ALPHA;
 598         else
 599                 alpha = GRID_DAY_ALPHA;
 600
 601         cairo_save(robject->cr);
 602         cairo_set_source_rgba(robject->cr, 0, 0.35, 0.55, alpha);
 603
 604         if (robject->type == RT_FAST)
 605                 cairo_set_tolerance (robject->cr, 1.0);
 606
 607         grid_render_ra(robject, proj, labels);
 608         grid_render_dec(robject, proj, labels);
 609
 610         cairo_stroke(robject->cr);
 611
 612         cairo_set_font_size (robject->cr, 13.0);
 613         grid_render_dec_labels(robject, proj);
 614         grid_render_ra_labels(robject, proj);
 615
 616         if (robject->type == RT_FAST)
 617                 cairo_set_tolerance (robject->cr, 0.1); /* do we need this */
 618         cairo_restore(robject->cr);
 619 }
 620
 621