# include "gastro2.h"

int Nv, NV;
float *xv, *yv;
Graphdata gv;

void plot_init_gridplot () {
  
  NV = 1000;
  Nv = 0;

  ALLOCATE (xv, float, NV);
  ALLOCATE (yv, float, NV);

  PlotReset (2);

  gv.xmin = -400;
  gv.xmax = +400;
  gv.ymin = -400;
  gv.ymax = +400;

  gv.style = 2;
  gv.ltype = 0;
  gv.etype = 0;
  gv.ebar  = 0;
  gv.size = 0.3;
  gv.lweight = 0;
  gv.color = 0;
  gv.ptype = 2;

}

int plot_addpt_gridplot (double x, double y) {

  if (x < gv.xmin) return (0);
  if (x > gv.xmax) return (0);
  if (y < gv.ymin) return (0);
  if (y > gv.ymax) return (0);

  xv[Nv] = x;
  yv[Nv] = y;
  Nv ++;
  
  if (Nv == NV) {
    NV += 1000;
    REALLOCATE (xv, float, NV);
    REALLOCATE (yv, float, NV);
  }
  return (1);
}

void plot_done_gridplot () {

  char c;

  /* send to Kapa */
  PrepPlotting (Nv, &gv, 2);
  PlotVector (Nv, xv, 0, 2);
  PlotVector (Nv, yv, 1, 2);

  DonePlotting (&gv, 2);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);
  free (xv);
  free (yv);
}
  
static double Xm[2];
static int Xg[2] = {0, 2};

void plot_resid_init (int version, double xmax) {

  if (version > 1) return;
  if (version < 0) return;
  
  Xm[version] = xmax;
}


void plot_resid_plot (int version, float *xvect, float *yvect, int Nvect) {
 
  int i;
  char c;
  float xmin, xmax, ymin, ymax;
  Graphdata graphdata;
  
  xmin = xmax = xvect[0];
  ymin = ymax = yvect[0];
  for (i = 0; i < Nvect; i++) {
    xmax = MAX (xvect[i], xmax);
    xmin = MIN (xvect[i], xmin);
    ymax = MAX (yvect[i], ymax);
    ymin = MIN (yvect[i], ymin);
  }
    
  PlotReset (Xg[version]);

  graphdata.xmin = MIN (xmin,   0);
  graphdata.xmax = MAX (xmax, Xm[version]);
  graphdata.ymin = MIN (ymin, -10);
  graphdata.ymax = MAX (ymax, +10);

  graphdata.style = 2;
  graphdata.ltype = 0;
  graphdata.etype = 0;
  graphdata.ebar  = 0;
  graphdata.size = 0.5;
  graphdata.lweight = 0;
  graphdata.color = 0;
  graphdata.ptype = 2;

  /* send to Kapa */
  PrepPlotting (Nvect, &graphdata, Xg[version]);
  PlotVector (Nvect, xvect, 0, Xg[version]);
  PlotVector (Nvect, yvect, 1, Xg[version]);

  DonePlotting (&graphdata, Xg[version]);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);
}

void plot_gridpts (double *pts, int Npts) {

  char c;
  int i;
  float *xvect, *yvect, ymax;
  Graphdata graphdata;
  
  ymax = 0;
  ALLOCATE (xvect, float, Npts);
  ALLOCATE (yvect, float, Npts);
  for (i = 0; i < Npts; i++) {
    xvect[i] = i;
    yvect[i] = pts[i];
    ymax = MAX (yvect[i], ymax);
  }
    
  PlotReset (0);

  graphdata.xmin = -10;
  graphdata.xmax = Npts + 10;
  graphdata.ymin = -1;
  graphdata.ymax = ymax + 1;

  graphdata.style = 1;
  graphdata.ltype = 0;
  graphdata.etype = 0;
  graphdata.ebar  = 0;
  graphdata.size = 0.5;
  graphdata.lweight = 0;
  graphdata.color = 0;
  graphdata.ptype = 0;

  /* send to Kapa */
  PrepPlotting (Npts, &graphdata, 0);
  PlotVector (Npts, xvect, 0, 0);
  PlotVector (Npts, yvect, 1, 0);
  free (xvect);
  free (yvect);

  DonePlotting (&graphdata, 0);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);
}

static Graphdata gf;

void plot_fullfield (CmpCatalog *Target, RefCatalog *Ref) {

  int i, Nvect;
  float *xvect, *yvect, *zvect, M, dM;
  char c;
  
  dM = Target[0].lum.Mmax - Target[0].lum.Mmin;

  PlotReset (1);

  gf.xmin = -100;
  gf.xmax = Target[0].header.Naxis[0] + 100;
  gf.ymin = -100;
  gf.ymax = Target[0].header.Naxis[1] + 100;

  gf.style = 2;
  gf.ltype = 0;
  gf.etype = 0;
  gf.ebar  = 0;
  gf.size = -1;
  gf.lweight = 0;

  /* fill in vectors */
  Nvect = Target[0].N;
  ALLOCATE (xvect, float, Nvect);
  ALLOCATE (yvect, float, Nvect);
  ALLOCATE (zvect, float, Nvect);
  for (i = 0; i < Nvect; i++) {
    xvect[i] = Target[0].stars[i].X;
    yvect[i] = Target[0].stars[i].Y;
    M = (Target[0].lum.Mmax - Target[0].stars[i].M) / dM;
    zvect[i] = MIN (1.0, MAX (0.01, M));
  }

  /* send to Kapa */
  gf.ptype = 1;
  gf.color = 0;
  PrepPlotting (Nvect, &gf, 1);
  PlotVector (Nvect, xvect, 0, 1);
  PlotVector (Nvect, yvect, 1, 1);
  PlotVector (Nvect, zvect, 1, 1);
  free (xvect);
  free (yvect);
  free (zvect);

  /* fill in vectors */
  Nvect = Ref[0].N;
  ALLOCATE (xvect, float, Nvect);
  ALLOCATE (yvect, float, Nvect);
  ALLOCATE (zvect, float, Nvect);
  for (i = 0; i < Nvect; i++) {
    xvect[i] = Ref[0].stars[i].X;
    yvect[i] = Ref[0].stars[i].Y;
    M = (Target[0].lum.Mmax - Ref[0].stars[i].M) / dM;
    zvect[i] = MIN (1.0, MAX (0.01, M));
  }

  /* send to Kapa */
  gf.color = 2;
  gf.ptype = 2;
  PrepPlotting (Nvect, &gf, 1);
  PlotVector (Nvect, xvect, 0, 1);
  PlotVector (Nvect, yvect, 1, 1);
  PlotVector (Nvect, zvect, 2, 1);
  free (xvect);
  free (yvect);
  free (zvect);

  DonePlotting (&gf, 1);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);
}

void plot_fullfield_pairs (float *x, float *y, int n) {

  char c;

  /* send to Kapa */
  gf.color = 6;
  gf.ptype = 100;
  gf.size =  1.0;
  
  PrepPlotting (n, &gf, 1);
  PlotVector (n, x, 0, 1);
  PlotVector (n, y, 1, 1);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);
}  

/* mag range -5 - 35, dmag = 0.25, Nbin = 160 */
# define MMIN -5
# define MMAX 35
# define dM 0.5
# define NMBIN 90

void plot_lumfunc (CmpCatalog *Target, RefCatalog *Ref) {

  int i, Nr, Nt;
  float tbin[NMBIN], tval[NMBIN], rbin[NMBIN], rval[NMBIN];
  double ymin, ymax;
  char c;
  Graphdata graphdata;

  fill_lumfunc (Target[0].stars, Target[0].N, tval, tbin, &Nt);
  fill_lumfunc (Ref[0].stars, Ref[0].N, rval, rbin, &Nr);

  ymin = 5; ymax = -5;
  for (i = 0; i < Nt; i++) {
    tval[i] = tval[i] - log (Target[0].Area) / log (10.0);
    ymin = MIN (tval[i], ymin);
    ymax = MAX (tval[i], ymax);
  }

  for (i = 0; i < Nr; i++) {
    rval[i] = rval[i] - log (Ref[0].Area) / log (10.0);
    rbin[i] = rbin[i] + Ref[0].Moff;
    ymin = MIN (rval[i], ymin);
    ymax = MAX (rval[i], ymax);
  }

  graphdata.xmin = MMIN;
  graphdata.xmax = MMAX;
  graphdata.ymin = ymin - 0.1;
  graphdata.ymax = ymax + 0.1;

  graphdata.style = 1;
  graphdata.ptype = 2;
  graphdata.ltype = 0;
  graphdata.etype = 0;
  graphdata.ebar  = 0;
  graphdata.color = 0;

  graphdata.lweight = 0;
  graphdata.size = 0.5;

  PlotReset (0);

  PrepPlotting (Nt, &graphdata, 0);
  PlotVector (Nt, tbin, 0, 0);
  PlotVector (Nt, tval, 1, 0);

  graphdata.color = 2;
  graphdata.style = 1;
  graphdata.lweight = 0;
  PrepPlotting (Nr, &graphdata, 0);
  PlotVector (Nr, rbin, 0, 0);
  PlotVector (Nr, rval, 1, 0);

  /* plot truncated lum func */
  for (i = 0; i < Nr; i++) {
    if (rbin[i] < Target[0].lum.Mmin - 0.5) rval[i] = -1;
    if (rbin[i] > Target[0].lum.Mmax + 0.5) rval[i] = -1;
  }
  graphdata.ltype = 1;
  PrepPlotting (Nr, &graphdata, 0);
  PlotVector (Nr, rbin, 0, 0);
  PlotVector (Nr, rval, 1, 0);

  DonePlotting (&graphdata, 0);
  fprintf (stderr, "type return to continue");
  fscanf (stdin, "%c", &c);

}

void fill_lumfunc (StarData *stars, int N, float *lbin, float *bin, int *nb) {

  int i, j, Nb;
  double mbin[NMBIN];

  bzero (mbin, NMBIN * sizeof (double));

  /* sum histogram */
  for (i = 0; i < N; i++) {
    if (stars[i].M < MMIN) continue;
    if (stars[i].M > MMAX) continue;

    j = (stars[i].M - MMIN) / dM;
    j = MIN (MAX (j, 0), (NMBIN - 1));
    mbin[j] ++;
  }

  /* select filled bins */
  for (Nb = i = 0; i < NMBIN; i++) {
    if (mbin[i] > 0) {
      bin[Nb]  = i * dM + MMIN;
      lbin[Nb] = log (mbin[i]) / log (10.0);
      Nb++;
    }
  }

  *nb = Nb;
}

