# include "setgalmodel.h"

int update_catalog_setgalmodel (Catalog *catalog) {

  off_t i, j, m;

  // choose this number more judiciously...
  time_t timeRef = ohana_date_to_sec("2000/01/01,00:00:00");

  Average *average = catalog[0].average;

  if ((catalog[0].Nstarpar_disk > 0) && (catalog[0].Nstarpar == 0)) {
    fprintf (stderr, "ERROR: failed to load starpar data\n");
    exit (2);
  }

  CoordTransform *transform = InitTransform (COORD_CELESTIAL, COORD_GALACTIC);

  // first set the proper motion based on Galactic rotation and solar motion
  for (i = 0; i < catalog[0].Naverage; i++) {

    { 
      m = average[i].measureOffset;
      Measure *measure = &catalog[0].measure[m];
      for (j = 0; j < average[i].Nmeasure; j++) {
	measure[j].RoffGAL = 0.0; // XXX divide by dcos? project and calculate?
	measure[j].DoffGAL = 0.0; // XXX divide by dcos? project and calculate?
      }
    }

    if (average[i].Nstarpar == 0) continue;

    m = average[i].starparOffset;
    StarPar *starpar = &catalog[0].starpar[m];

    // fake or real QSOs are marked with FeH = +/- 100.0
    if (fabs(starpar->FeH) > 99.0) continue;

    double distance = pow(10.0, 0.2*(starpar->DistMag + 5.0));

    double Lrad = starpar->galLon * RAD_DEG;
    double Brad = starpar->galLat * RAD_DEG;

    double uL_gal, uB_gal;
    GalaxyMotionModel_radians(&uL_gal, &uB_gal, Lrad, Brad);

    double uL_sol, uB_sol;
    SolarMotionModel_radians(&uL_sol, &uB_sol, Lrad, Brad, distance);

    double uL = uL_gal + uL_sol;
    double uB = uB_gal + uB_sol;

    // XXX: amplify motion to make tests easier:
    if (TESTING) {
      uL *= 100.0;
      uB *= 100.0;
    }
    
    double uR, uD;
    TransformProperMotion (&uR, &uD, uL, uB, average[i].R, average[i].D, transform);

    starpar->uRA  = uR;
    starpar->uDEC = uD;

    m = average[i].measureOffset;
    Measure *measure = &catalog[0].measure[m];

    // use the model uR,uD values to move the expected position based on the reference date;
    for (j = 0; j < average[i].Nmeasure; j++) {

      double dtime = (measure[j].t - timeRef) / 86400.0 / 365.25;

      // do not modify the original Xccd,Yccd values
      measure[j].RoffGAL = uR * dtime; // XXX divide by dcos? project and calculate?
      measure[j].DoffGAL = uD * dtime; // XXX divide by dcos? project and calculate?
    }
  }

  return (TRUE);
}

