# include "relphot.h"
   
void plot_scatter (Catalog *catalog, int Ncatalog) {

  int i, j, k, m, N, Ntot;
  float Mrel, Mcal, Mmos, Mgrid;
  double *xlist, *ylist, *ilist;
  Graphdata graphdata;

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Ntot += catalog[i].average[j].Nmeasure;
    }
  }
  ALLOCATE (xlist, double, Ntot);
  ALLOCATE (ylist, double, Ntot);
  ALLOCATE (ilist, double, Ntot);

  N = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].average[j].flags & STAR_BAD) continue;  
      m = catalog[i].average[j].measureOffset;

      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
	Mcal = getMcal  (m, i);
	if (isnan(Mcal)) continue;
	Mmos = getMmos  (m, i);
	if (isnan(Mmos)) continue;
	Mgrid = getMgrid (m, i);
	if (isnan(Mgrid)) continue;

	Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
	xlist[N] = Mrel;
	ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;
	ilist[N] = PhotInst (&catalog[i].measure[m]);
	N++;
      }
    }
  }

  plot_defaults (&graphdata);
  graphdata.xmin = PlotMmin;
  graphdata.xmax = PlotMmax;
  graphdata.ymin = PlotdMmin;
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "%s.Mag.png", OUTROOT);

  plot_defaults (&graphdata);
  graphdata.ymin = PlotdMmin;
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "%s.iMag.png", OUTROOT);
  free (xlist);
  free (ylist);
  free (ilist);
}
