# 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].Nm;
    }
  }
  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].code & STAR_BAD) continue;  
      m = catalog[i].average[j].offset;

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

	Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
	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", "Mag.png");

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