# include "astro.h"

int czplot (int argc, char **argv) {
  
  int i, N, Npts, Ngraph;
  double min, range, Rmin, Rmax;
  float *in, *out, *r, *d, *x, *y;
  Vector Xvec, Yvec, Zvec, *xvec, *yvec, *zvec;
  Graphdata graphmode;

  Ngraph = -1;
  if (N = get_argument (argc, argv, "-n")) {
    remove_argument (N, &argc, argv);
    Ngraph = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }
  if (!GetGraph (&graphmode, NULL, &Ngraph)) return (FALSE);

  if (argc != 6) {
    fprintf (stderr, "USAGE: czplot <x> <y> <z> min max\n");
    return (FALSE);
  }

  min = atof(argv[4]);
  range = atof(argv[5]) - min;
  Rmin = graphmode.coords.crval1 - 180.0;
  Rmax = graphmode.coords.crval1 + 180.0;

  /* find vectors */
  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
  if ((zvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);

  if (xvec[0].Nelements != yvec[0].Nelements) {
    fprintf (stderr, "vectors %s and %s not the same length\n", argv[1], argv[2]);
    return (FALSE);
  }
  if (xvec[0].Nelements != zvec[0].Nelements) {
    fprintf (stderr, "vectors %s and %s not the same length\n", argv[1], argv[3]);
    return (FALSE);
  }
  Xvec.Nelements = xvec[0].Nelements;
  Yvec.Nelements = xvec[0].Nelements;
  Zvec.Nelements = zvec[0].Nelements;
  ALLOCATE (Xvec.elements, float, Xvec.Nelements);
  ALLOCATE (Yvec.elements, float, Yvec.Nelements);
  ALLOCATE (Zvec.elements, float, Zvec.Nelements);
  
  r   = xvec[0].elements;
  d   = yvec[0].elements;
  in  = zvec[0].elements;
  x   = Xvec.elements;
  y   = Yvec.elements;
  out = Zvec.elements;
  for (i = 0; i < Zvec.Nelements; i++, in++, out++, r++, d++, x++, y++) {
    *out = MIN (1.0, MAX (0.01, (*in - min) / range));
    while (*r < Rmin) *r += 360.0;
    while (*r > Rmax) *r -= 360.0;
    fRD_to_XY (x, y, *r, *d, &graphmode.coords);
  }

  graphmode.style = 2;
  graphmode.size = -1; /* point size determined by Zvec */
  graphmode.etype = 0;
  Npts = Xvec.Nelements;
  PrepPlotting (Npts, &graphmode);

  PlotVector (Npts, Xvec.elements);
  PlotVector (Npts, Yvec.elements);
  PlotVector (Npts, Zvec.elements);

  free (Xvec.elements);
  free (Yvec.elements);
  free (Zvec.elements);

  return (TRUE);

}
