text: /*
   # This file is part of the Astrometry.net suite.
   # Licensed under a 3-clause BSD style license - see LICENSE
   */
  #include <string.h>
  #include <math.h>
  
  #include "os-features.h"
  #include "plothealpix.h"
  #include "sip-utils.h"
  #include "cairoutils.h"
  #include "ioutils.h"
  #include "log.h"
  #include "errors.h"
  #include "healpix-utils.h"
  #include "healpix.h"
  
  DEFINE_PLOTTER(healpix);
  
  plothealpix_t* plot_healpix_get(plot_args_t* pargs) {
      return plotstuff_get_config(pargs, "healpix");
  }
  
  void* plot_healpix_init(plot_args_t* plotargs) {
      plothealpix_t* args = calloc(1, sizeof(plothealpix_t));
      args->nside = 1;
      args->stepsize = 50;
      return args;
  }
  
  int plot_healpix_plot(const char* command,
                        cairo_t* cairo, plot_args_t* pargs, void* baton) {
      plothealpix_t* args = (plothealpix_t*)baton;
      double ra,dec,rad;
      il* hps;
      int i;
      double hpstep;
      int minx[12], maxx[12], miny[12], maxy[12];
  
      plotstuff_builtin_apply(cairo, pargs);
  
      if (plotstuff_get_radec_center_and_radius(pargs, &ra, &dec, &rad)) {
          ERROR("Failed to get RA,Dec center and radius");
          return -1;
      }
      hps = healpix_rangesearch_radec(ra, dec, rad, args->nside, NULL);
      logmsg("Found %zu healpixes in range.\n", il_size(hps));
      hpstep = args->nside * args->stepsize * plotstuff_pixel_scale(pargs) / 60.0 / healpix_side_length_arcmin(args->nside);
      hpstep = MIN(1, hpstep);
      logmsg("Taking steps of %g in healpix space\n", hpstep);
  
      // For each of the 12 top-level healpixes, find the range of healpixes covered by this image.
      for (i=0; i<12; i++) {
          maxx[i] = maxy[i] = -1;
          minx[i] = miny[i] = args->nside+1;
      }
      for (i=0; i<il_size(hps); i++) {
          int hp = il_get(hps, i);
          int hpx, hpy;
          int bighp;
          healpix_decompose_xy(hp, &bighp, &hpx, &hpy, args->nside);
          logverb("  hp %i: bighp %i, x,y (%i,%i)\n", i, bighp, hpx, hpy);
          minx[bighp] = MIN(minx[bighp], hpx);
          maxx[bighp] = MAX(maxx[bighp], hpx);
          miny[bighp] = MIN(miny[bighp], hpy);
          maxy[bighp] = MAX(maxy[bighp], hpy);
      }
      il_free(hps);
  
      for (i=0; i<12; i++) {
          int hx,hy;
          int hp;
          double d, frac;
          double x,y;
  
          if (maxx[i] == -1)
              continue;
          logverb("Big healpix %i: x range [%i, %i], y range [%i, %i]\n",
                  i, minx[i], maxx[i], miny[i], maxy[i]);
  
          for (hy = miny[i]; hy <= maxy[i]; hy++) {
              logverb("  y=%i\n", hy);
              for (d=minx[i]; d<=maxx[i]; d+=hpstep) {
                  hx = floor(d);
                  frac = d - hx;
                  hp = healpix_compose_xy(i, hx, hy, args->nside);
                  healpix_to_radecdeg(hp, args->nside, frac, 0.0, &ra, &dec);
                  if (!plotstuff_radec2xy(pargs, ra, dec, &x, &y))
                      continue;
                  if (d == minx[i])
                      cairo_move_to(pargs->cairo, x, y);
                  else
                      cairo_line_to(pargs->cairo, x, y);
              }
              cairo_stroke(pargs->cairo);
          }
          for (hx = minx[i]; hx <= maxx[i]; hx++) {
              for (d=miny[i]; d<=maxy[i]; d+=hpstep) {
                  hy = floor(d);
                  frac = d - hy;
                  hp = healpix_compose_xy(i, hx, hy, args->nside);
                  healpix_to_radecdeg(hp, args->nside, 0.0, frac, &ra, &dec);
                  if (!plotstuff_radec2xy(pargs, ra, dec, &x, &y))
                      continue;
                  if (d == miny[i])
                      cairo_move_to(pargs->cairo, x, y);
                  else
                      cairo_line_to(pargs->cairo, x, y);
              }
              cairo_stroke(pargs->cairo);
          }
      }
      return 0;
  }
  
  int plot_healpix_command(const char* cmd, const char* cmdargs,
                           plot_args_t* pargs, void* baton) {
      plothealpix_t* args = (plothealpix_t*)baton;
      if (streq(cmd, "healpix_nside")) {
          args->nside = atoi(cmdargs);
      } else if (streq(cmd, "healpix_stepsize")) {
          args->stepsize = atoi(cmdargs);
      } else {
          ERROR("Did not understand command \"%s\"", cmd);
          return -1;
      }
      return 0;
  }
  
  void plot_healpix_free(plot_args_t* plotargs, void* baton) {
      plothealpix_t* args = (plothealpix_t*)baton;
      free(args);
  }
  
  