# include "gastro2.h"

int getptolemy (CatStats *catstats, RefCatalog *Ref) {
  
  int i, j, k, Ns, Nptolemy, Nregions; 
  double FracArea;
  GSCdata *ptolemy;
  CatStats *regions;

  Ref[0].N = 0;
  Ref[0].Area = 0;
  ALLOCATE (Ref[0].stars, StarData, 1);

  regions = gregions (catstats, &Nregions);
  
  for (i = 0; i < Nregions; i++) {
    ptolemy = gptolemy (regions[i].name, &Nptolemy);
    if (Nptolemy == 0) continue;

    FracArea = ptolemy_area (ptolemy, Nptolemy, &regions[i]);

    Ns = Ref[0].N;
    Ref[0].N += Nptolemy;
    Ref[0].Area += regions[i].Area * FracArea;

    REALLOCATE (Ref[0].stars, StarData, MAX (1, Ref[0].N));
    for (k = Ns, j = 0; j < Nptolemy; k++, j++) {
      Ref[0].stars[k].R = ptolemy[j].R;
      Ref[0].stars[k].D = ptolemy[j].D;
      Ref[0].stars[k].M = ptolemy[j].M;
    }      
    free (ptolemy);
  }
  
  Ref[0].R0 = catstats[0].RA[0];
  Ref[0].R1 = catstats[0].RA[1];
  Ref[0].D0 = catstats[0].DEC[0];
  Ref[0].D1 = catstats[0].DEC[1];

  /* calculate luminosity function of stars */
  get_luminosity_func (Ref[0].stars, Ref[0].N, &Ref[0].lum);

  if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", Ref[0].N);
  return (TRUE);
}  

double ptolemy_area (GSCdata *ptolemy, int N, CatStats *patch) {

  int i, xb, yb, Nb;
  int bin[10][10];
  double frac, Rmin, Rmax, Dmin, Dmax, dR, dD;

  Rmin = patch[0].RA[0];
  Rmax = patch[0].RA[1];
  Dmin = patch[0].DEC[0];
  Dmax = patch[0].DEC[1];
  dR = Rmax - Rmin;
  dD = Dmax - Dmin;

  for (xb = 0; xb < 10; xb++) {
    for (yb = 0; yb < 10; yb++) {
      bin[xb][yb] = 0;
    }
  }

  for (i = 0; i < N; i++) {
    xb = MAX (MIN (0, 10 * (ptolemy[i].R - Rmin) / dR), 9);
    yb = MAX (MIN (0, 10 * (ptolemy[i].D - Dmin) / dD), 9);
    bin[xb][yb] ++;
  }

  Nb = 0;
  for (xb = 0; xb < 10; xb++) {
    for (yb = 0; yb < 10; yb++) {
      if (bin[xb][yb]) Nb ++;
    }
  }

  frac = Nb / 100.0;

  return (frac);
}
