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

int cgrid (int argc, char **argv) {
  
  double range, lrange, factor, mantis, power, fmantis;
  double firstRA, firstDEC, minorRA, minorDEC;
  double r, d, dR, dD, D;
  double x, y;
  Vector Xvec, Yvec;
  int kapa, NorthPole, SouthPole, N, OnPic, LOnPic, status, NELEMENTS, First;
  Graphdata graphmode;

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

  if (argc != 1) {
    gprint (GP_ERR, "USAGE: cgrid [style]\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;
  lrange = log10(MAX(fabs(range), 1e-30));
  factor = (int) (lrange);
  if (lrange < 0) { factor -= 1; }
  mantis = lrange - factor;
  power = pow(10.0, factor);
  fmantis = pow(10.0, mantis);
  if ((fmantis >= 1.0) && (fmantis <=  2.0)) {
    minorRA = 0.1 * power;
  }
  if ((fmantis > 2.0) && (fmantis <=  4.0)) {
    minorRA = 0.2 * power;
  }
  if ((fmantis > 4.0) && (fmantis <=  6.0)) {
    minorRA = 0.5 * power;
  }
  if ((fmantis > 6.0) && (fmantis <=  10.0)) {
    minorRA = 0.5 * power;
  }
  dR = range / 100.0;

  /* set spacings for DEC */
  range = MIN (fabs(graphmode.coords.cdelt1*(graphmode.xmax-graphmode.xmin)), fabs(graphmode.coords.cdelt2*(graphmode.ymax-graphmode.ymin)));
  lrange = log10(MAX(fabs(range), 1e-30));
  factor = (int) (lrange);
  if (lrange < 0) { factor -= 1; }
  mantis = lrange - factor;
  power = pow(10.0, factor);
  fmantis = pow(10.0, mantis);
  if ((fmantis >= 1.0) && (fmantis <=  2.0)) {
    minorDEC = 0.1 * power;
  }
  if ((fmantis > 2.0) && (fmantis <=  4.0)) {
    minorDEC = 0.2 * power;
  }
  if ((fmantis > 4.0) && (fmantis <=  6.0)) {
    minorDEC = 0.5 * power;
  }
  if ((fmantis > 6.0) && (fmantis <=  10.0)) {
    minorDEC = 0.5 * power;
  }
  dD = range / 100.0;

  /* 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);
  
  /***  do consecutive RA lines, first increasing **/
  OnPic = TRUE;
  for (r = firstRA; (r <= firstRA + 180) && (OnPic); r += minorRA) {
    /* 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));
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
    /* next, DEC decreasing */
    First = TRUE;
    LOnPic = TRUE;
    for (d = firstDEC; (d > -90 - dD) && (LOnPic || NorthPole || SouthPole); d -= dD) {
      D = MAX (-90, MIN(90, d));
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      } 
    }
  }

  /***  do consecutive RA lines, decreasing **/
  OnPic = TRUE;
  for (r = firstRA; (r >=  firstRA - 180) && (OnPic); r -= minorRA) {
    /* first, DEC increasing */
    First = TRUE;
    LOnPic = TRUE;
    OnPic = FALSE;
    for (d = firstDEC; (d < 90 + dD) && (LOnPic || NorthPole || SouthPole); d += dD) {
      D = MAX (-90, MIN(90, d));
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
    /* next, DEC decreasing */
    First = TRUE;
    LOnPic = TRUE;
    for (d = firstDEC; (d > -90 - dD) && (LOnPic || NorthPole || SouthPole); d -= dD) {
      D = MAX (-90, MIN(90, d));
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
  }

  /***  do consecutive DEC lines, first increasing **/
  OnPic = TRUE;
  for (d = firstDEC; (d < 90 + dD) && (OnPic); d += minorDEC) {
    D = MAX (-90, MIN(90, d));
    /* first, RA increasing */
    LOnPic = TRUE;
    OnPic = FALSE;
    First = TRUE;
    for (r = firstRA; (r < firstRA + 180) && (LOnPic || NorthPole || SouthPole); r += dR) {
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
    /* next, RA decreasing */
    First = TRUE;
    LOnPic = TRUE;
    for (r = firstRA; (r > firstRA - 180) && (LOnPic || NorthPole || SouthPole); r -= dR) {
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
  }

  /***  do consecutive DEC lines, decreasing **/
  OnPic = TRUE;
  for (d = firstDEC; (d > -90 - dD) && (OnPic); d -= minorDEC) {
    D = MAX (-90, MIN(90, d));
    /* first, RA increasing */
    LOnPic = TRUE;
    OnPic = FALSE;
    First = TRUE;
    for (r = firstRA; (r < firstRA + 180) && (LOnPic || NorthPole || SouthPole); r += dR) {
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
    /* next, RA decreasing */
    First = TRUE;
    LOnPic = TRUE;
    for (r = firstRA; (r > firstRA - 180) && (LOnPic || NorthPole || SouthPole); r -= dR) {
      status = RD_to_XY (&Xvec.elements.Flt[N], &Yvec.elements.Flt[N], r, D, &graphmode.coords);
      /*
      if ((fabs(Xvec.elements.Flt[N] - Xvec.elements.Flt[N-1]) > 10) && (fabs(Yvec.elements.Flt[N] - Yvec.elements.Flt[N-1]) > 10))
	First = TRUE;
	*/
      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;
      }
    }
  }
  
  /* 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, &graphmode);

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

}
