# include "astro.h"

# define GRID_SPACING 0.002

# define CHECKELEMENTS						\
  if (N == NELEMENTS) {						\
    NELEMENTS +=200;						\
    REALLOCATE (Xvec.elements.Flt, opihi_flt, NELEMENTS);	\
    REALLOCATE (Yvec.elements.Flt, opihi_flt, NELEMENTS);	\
  }

# define ADD_COORDINATE(RA,DEC)						\
  status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], (RA), (DEC), &graphmode.coords); \
  if ((Xvec.elements.Flt[N] >= graphmode.xmin) && (Xvec.elements.Flt[N] <= graphmode.xmax) && \
      (Yvec.elements.Flt[N] >= graphmode.ymin) && (Yvec.elements.Flt[N] <= graphmode.ymax) && status) { \
    N++;								\
    CHECKELEMENTS;							\
    OnPic = TRUE;							\
    if (!First) {							\
      Xvec.elements.Flt[N] = Xvec.elements.Flt[N-1];			\
      Yvec.elements.Flt[N] = Yvec.elements.Flt[N-1];			\
      N++;								\
      CHECKELEMENTS;							\
    } else {								\
      if (N > 1) {							\
	Xvec.elements.Flt[N-2] = Xvec.elements.Flt[N-1];		\
	Yvec.elements.Flt[N-2] = Yvec.elements.Flt[N-1];		\
	N--;								\
      }									\
      First = FALSE;							\
    }									\
  } else {								\
    LOnPic = FALSE;							\
    First = TRUE;							\
  }

# define ADD_DEC_LINE(RA)						\
  /* first, DEC increasing */						\
  LOnPic = TRUE;							\
  OnPic = FALSE;							\
  First = TRUE;								\
  for (d = firstDEC; (d < 90 + dD) && (LOnPic || NorthPole || SouthPole); d += dD) { \
    D = MAX (-90, MIN(90, d));						\
    ADD_COORDINATE((RA), D);						\
  }									\
  /* next, DEC decreasing */						\
  First = TRUE;								\
  LOnPic = TRUE;							\
  for (d = firstDEC; (d > -90 - dD) && (LOnPic || NorthPole || SouthPole); d -= dD) { \
    D = MAX (-90, MIN(90, d));						\
    ADD_COORDINATE((RA), D);						\
  } 

# define ADD_RA_LINE(DEC)						\
  D = MAX (-90, MIN(90, (DEC)));					\
  /* first, RA increasing */						\
  LOnPic = TRUE;							\
  OnPic = FALSE;							\
  First = TRUE;								\
  lastRA = graphmode.coords.crval1 + 180.0;				\
  for (r = firstRA; (r < lastRA + dR) && (LOnPic || NorthPole || SouthPole); r += dR) { \
    R = MIN (r, lastRA);					\
    ADD_COORDINATE(R, D);						\
  }									\
  /* next, RA decreasing */						\
  First = TRUE;								\
  LOnPic = TRUE;							\
  lastRA = graphmode.coords.crval1 - 180.0;				\
  for (r = firstRA; (r > lastRA - dR) && (LOnPic || NorthPole || SouthPole); r -= dR) { \
    R = MAX (r, lastRA);					\
    ADD_COORDINATE(R, D);						\
  }

int cgrid (int argc, char **argv) {
  
  double range, minor, major;
  double firstRA, lastRA, firstDEC, minorRA, minorDEC;
  double r, d, dR, dD, R, D;
  double x, y;
  Vector Xvec, Yvec;
  int kapa, NorthPole, SouthPole, N, OnPic, LOnPic, status, NELEMENTS;
  int First, RAbyHour, Labels;
  Graphdata graphmode;

  RAbyHour = FALSE;
  if ((N = get_argument (argc, argv, "-ra-by-hour"))) {
    remove_argument (N, &argc, argv);
    RAbyHour = TRUE;
  }

  Labels = FALSE;
  if ((N = get_argument (argc, argv, "-labels"))) {
    remove_argument (N, &argc, argv);
    Labels = TRUE;
  }

  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;

  if (argc != 1) {
    gprint (GP_ERR, "USAGE: cgrid [style] [-ra-by-hour] [-labels]\n");
    return (FALSE);
  }

  /* are we plotting one of the poles? */
  NorthPole = SouthPole = FALSE;
  status = RD_to_XY (&x, &y, 0.0, 90.0, &graphmode.coords);
  if ((x >= graphmode.xmin) && (x <= graphmode.xmax) && 
      (y >= graphmode.ymin) && (y <= graphmode.ymax) && status)
    NorthPole = TRUE;
  status = RD_to_XY (&x, &y, 0.0, -90.0, &graphmode.coords);
  if ((x >= graphmode.xmin) && (x <= graphmode.xmax) && 
      (y >= graphmode.ymin) && (y <= graphmode.ymax) && status)
    SouthPole = TRUE;

  /* set spacings for RA */
  minorRA = minorDEC = 0.1;
  range = MIN (fabs(graphmode.coords.cdelt1*(graphmode.xmax-graphmode.xmin)), fabs(graphmode.coords.cdelt2*(graphmode.ymax-graphmode.ymin)));
  if (NorthPole || SouthPole) range = 360;

  if (RAbyHour) {
    SetGridScales (&major, &minor, range / 15.0);
    minorRA = minor * 15.0;
  } else {
    SetGridScales (&major, &minor, range);
    minorRA = minor;
  }
  dR = range * GRID_SPACING;

  /* set spacings for DEC */
  range = MIN (fabs(graphmode.coords.cdelt1*(graphmode.xmax-graphmode.xmin)), fabs(graphmode.coords.cdelt2*(graphmode.ymax-graphmode.ymin)));
  SetGridScales (&major, &minor, range);
  minorDEC = minor;
  dD = range * GRID_SPACING;

  /* choose a starting point */
  if ((int)(graphmode.coords.crval1/minorRA) == (graphmode.coords.crval1/minorRA)) {
    firstRA = graphmode.coords.crval1;
  } else {
    firstRA = minorRA + minorRA*((int)(graphmode.coords.crval1/minorRA));
  }
  if ((int)(graphmode.coords.crval2/minorDEC) == (graphmode.coords.crval2/minorDEC)) {
    firstDEC = graphmode.coords.crval2;
  } else {
    firstDEC = minorDEC + minorDEC*((int)(graphmode.coords.crval2/minorDEC));
  }
  if (SouthPole) firstDEC = -90;
  if (NorthPole) firstDEC = 90;
  
  /* prepare vectors to hold data */
  N = 0;
  NELEMENTS = 200;
  SetVector (&Xvec, OPIHI_FLT, NELEMENTS);
  SetVector (&Yvec, OPIHI_FLT, NELEMENTS);
  
  { // sanity check
    float Nelem;
    Nelem = 180.0 / minorRA;
    if (!isfinite(Nelem) || (fabs(Nelem) > 10000)) { fprintf (stderr, "absurd cgrid spacing %f\n", minorRA); return FALSE; }
    Nelem = 180.0 / minorDEC;
    if (!isfinite(Nelem) || (fabs(Nelem) > 10000)) { fprintf (stderr, "absurd cgrid spacing %f\n", minorDEC); return FALSE; }
  }

  /***  do consecutive RA lines, first increasing **/
  OnPic = TRUE;
  lastRA = graphmode.coords.crval1 + 180.0;
  for (r = firstRA; (r <= lastRA) && (OnPic); r += minorRA) {
    ADD_DEC_LINE (r);
  }
  if (r != lastRA) {
    ADD_DEC_LINE (lastRA);
  }

  /***  do consecutive RA lines, decreasing **/
  OnPic = TRUE;
  lastRA = graphmode.coords.crval1 - 180.0;
  for (r = firstRA; (r >=  lastRA) && (OnPic); r -= minorRA) {
    ADD_DEC_LINE (r);
  }
  if (r != lastRA) {
    ADD_DEC_LINE (lastRA);
  }

  /***  do consecutive DEC lines, first increasing **/
  OnPic = TRUE;
  for (d = firstDEC; (d < 90 + dD) && (OnPic); d += minorDEC) {
    ADD_RA_LINE (d);
  }

  /***  do consecutive DEC lines, decreasing **/
  OnPic = TRUE;
  for (d = firstDEC; (d > -90 - dD) && (OnPic); d -= minorDEC) {
    ADD_RA_LINE (d);
  }
  
  // add labels for center lines:
  if (Labels) { 
    char line[16], format[8];
    double xt, yt, dx, dy, frac;
    dx = +0.01 * (graphmode.xmax - graphmode.xmin);
    dy = -0.02 * (graphmode.ymax - graphmode.ymin);

    for (r = firstRA; r <= graphmode.coords.crval1 + 180.0; r += minorRA) {
      status = RD_to_XY (&xt, &yt, r, graphmode.coords.crval2, &graphmode.coords);
      if (!status) continue;
      if (xt < graphmode.xmin) continue;
      if (xt > graphmode.xmax) continue;
      if (yt < graphmode.ymin) continue;
      if (yt > graphmode.ymax) continue;
      frac = -1.0 * log10(minorRA);
      if (frac != (int)frac) {
	frac += 1.0;
      }
      if (frac <= 0.0) frac = 0.0;
      snprintf (format, 8, "%%.%df", (int) frac);
      if (RAbyHour) {
	snprintf (line, 16, format, r / 15.0);
      } else {
	snprintf (line, 16, format, r);
      }
      KapaSendTextline (kapa, line, xt + dx, yt + dy, 0.0);
    }
    for (r = firstRA; r >= graphmode.coords.crval1 - 180.0; r -= minorRA) {
      status = RD_to_XY (&xt, &yt, r, graphmode.coords.crval2, &graphmode.coords);
      if (!status) continue;
      if (xt < graphmode.xmin) continue;
      if (xt > graphmode.xmax) continue;
      if (yt < graphmode.ymin) continue;
      if (yt > graphmode.ymax) continue;
      frac = -1.0 * log10(minorRA);
      if (frac != (int)frac) {
	frac += 1.0;
      }
      if (frac <= 0.0) frac = 0.0;
      snprintf (format, 8, "%%.%df", (int) frac);
      if (RAbyHour) {
	snprintf (line, 16, format, r / 15.0);
      } else {
	snprintf (line, 16, format, r);
      }
      KapaSendTextline (kapa, line, xt + dx, yt + dy, 0.0);
    }
    for (d = firstDEC; d <= graphmode.coords.crval2 + 90.0; d += minorDEC) {
      status = RD_to_XY (&xt, &yt, graphmode.coords.crval1, d, &graphmode.coords);
      if (!status) continue;
      if (xt < graphmode.xmin) continue;
      if (xt > graphmode.xmax) continue;
      if (yt < graphmode.ymin) continue;
      if (yt > graphmode.ymax) continue;
      frac = -1.0 * log10(minorDEC);
      if (frac != (int)frac) {
	frac += 1.0;
      }
      if (frac <= 0.0) frac = 0.0;
      snprintf (format, 8, "%%.%df", (int) frac);
      snprintf (line, 16, format, d);
      KapaSendTextline (kapa, line, xt + dx, yt + dy, 0.0);
    }
    for (d = firstDEC; d >= graphmode.coords.crval2 - 90.0; d -= minorDEC) {
      status = RD_to_XY (&xt, &yt, graphmode.coords.crval1, d, &graphmode.coords);
      if (!status) continue;
      if (xt < graphmode.xmin) continue;
      if (xt > graphmode.xmax) continue;
      if (yt < graphmode.ymin) continue;
      if (yt > graphmode.ymax) continue;
      frac = -1.0 * log10(minorDEC);
      if (frac != (int)frac) {
	frac += 1.0;
      }
      if (frac <= 0.0) frac = 0.0;
      snprintf (format, 8, "%%.%df", (int) frac);
      snprintf (line, 16, format, d);
      KapaSendTextline (kapa, line, xt + dx, yt + dy, 0.0);
    }
  }

  /* send the line segments as connect-points */
  Xvec.Nelements = Yvec.Nelements = N;
  graphmode.style = 2; /* points */
  graphmode.ptype = 100; /* connect a pair */
  graphmode.etype = 0;
  PlotVectorPair (kapa, &Xvec, &Yvec, NULL, &graphmode);

  free (Xvec.elements.Ptr);
  free (Yvec.elements.Ptr);
  return (TRUE);

}
