/*****************************************************************************
******************************************************************************
**<AUTO>
**
** FILE:	catalogs.c
**
**<HTML>
**	This file contains C functions to match objects against other catalogs.
**</HTML>
**</AUTO>
**
**
** ENVIRONMENT:
**	ANSI C.
**
******************************************************************************
******************************************************************************
*/
#include <dervish.h>
#include <astrotools.h>
#include <slalib.h>
#include "taCalibObj.h"
#include "taCatalogs.h"
#include "arObjectStatus.h"


/* Binomial search by lambda on an array of TA_FIRST. Returns the index such
 * that the input lambda is in the range lambda(index) to lambda(index+1);
 * the input lambda is definitely greater than lambda(index), but could be
 * equal to lambda(index+1).
 */
int taFirstBinomialSearch(double lambda, TA_FIRST *array, int nObjects)
{
  int lo = -1;
  int hi = nObjects;
  while (hi - lo > 1)
    {
      int mid = (hi + lo) / 2;
      if (lambda > array[mid].lambda)
	lo = mid;
      else
	hi = mid;
    }
  return lo;
}

/* Binomial search by lambda on an array of TA_ROSAT. Returns the index such
 * that the input lambda is in the range lambda(index) to lambda(index+1);
 * the input lambda is definitely greater than lambda(index), but could be
 * equal to lambda(index+1).
 */
int taRosatBinomialSearch(double lambda, TA_ROSAT *array, int nObjects)
{
  int lo = -1;
  int hi = nObjects;
  while (hi - lo > 1)
    {
      int mid = (hi + lo) / 2;
      if (lambda > array[mid].lambda)
	lo = mid;
      else
	hi = mid;
    }
  return lo;
}

/* Binomial search by lambda on an array of TA_USNO. Returns the index such
 * that the input lambda is in the range lambda(index) to lambda(index+1);
 * the input lambda is definitely greater than lambda(index), but could be
 * equal to lambda(index+1).
 */
int taUSNOBinomialSearch(double lambda, TA_USNO *array, int nObjects)
{
  int lo = -1;
  int hi = nObjects;
  while (hi - lo > 1)
    {
      int mid = (hi + lo) / 2;
      if (lambda > array[mid].lambda)
	lo = mid;
      else
	hi = mid;
    }
  return lo;
}


/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	taObjMatchAgainstFirst
**
**<HTML>
**	Match a calibrated object (TA_CALIB_OBJ) against the FIRST catalog.
**	An object is a match if it isseparated from the optical source
**	by less than "matchRadius"
**
**	If a match is found, the fields "firstMatch", "firstDelta",
**	"firstPeak", "firstInt", "firstrms", "firstmajor", "firstminor",
**	and "firstpa" are filled in the TA_CALIB_OBJ object from the
**	matching TA_FIRST object; else those fields are all set to 0.
**
**	No FIRST-related target flags are modified (i.e., this routine
**	simply matches up, it does not select FIRST targets).
**</HTML>
**</AUTO>
******************************************************************************
****************************************************************************** */
void
taObjMatchAgainstFirst(TA_CALIB_OBJ *obj,
		       /* MODIFIED: calibrated object to match */
		       TA_FIRST_CATALOG *catalog
		       /* MOD: FIRST catalog */)
{
  /* Save position where we last started the search */
  static int previousStart;

  int index, closestIndex;
  double closestDistance, lambdaLo, lambdaHi, etaDiff, lambdaDiff;
  double cosLambda, distance, doubleSizeSq, matchRadiusSq;
  double ra,dec;
  int xpix, ypix;

  /* Set up useful numbers for matching */
  cosLambda = cos(obj->lambda * at_deg2Rad);

  /* Find the object at the smallest lambda which should be considered
     for the match algorithm --the max size of a composite radio source */

  lambdaLo = obj->lambda - catalog->doubleSize;
  if (catalog->newSearch)
    {
      previousStart = taFirstBinomialSearch(lambdaLo,
	catalog->obj, catalog->nObj) + 1;
      catalog->newSearch = 0;
    }
  else
    {
      while (previousStart > 0)
	{
	  if (catalog->obj[previousStart-1].lambda < lambdaLo)
	    break;
	  previousStart--;
	}
      while (previousStart < catalog->nObj)
	{
	  if (catalog->obj[previousStart].lambda >= lambdaLo)
	    break;
	  previousStart++;
	}
    }

  doubleSizeSq = catalog->doubleSize * catalog->doubleSize;
  matchRadiusSq = catalog->matchRadius * catalog->matchRadius;

  /* Now check each object until we're beyond the match radius or we run
   * out of the catalog */
  lambdaHi = obj->lambda + catalog->doubleSize;
  closestDistance = 1000. * doubleSizeSq;
  closestIndex = -1;
  index = previousStart;
  while (index < catalog->nObj)
    {
      /* Break out if we're beyond the radius we are interested in */
      if (catalog->obj[index].lambda > lambdaHi) break;

      etaDiff = obj->eta - catalog->obj[index].eta;
      if (etaDiff < -180) etaDiff += 360.;
      if (etaDiff > 180) etaDiff -= 360.;
      etaDiff *= cosLambda;
      if (fabs( etaDiff ) <= catalog->doubleSize)
	{
	  lambdaDiff = (obj->lambda - catalog->obj[index].lambda);
	  distance = lambdaDiff * lambdaDiff + etaDiff * etaDiff;
	  if (distance <= doubleSizeSq)
	    {
	      /* ----------------------------------------- */
	      /* We have a source within doubleSize */

	      /* wflag is a "warning" flag which indicates that this
               * object might be a sidelobe of a bright source nearby.
               * We do not want to include such objects. */
	      if (catalog->obj[index].wflag == 0 && distance <= matchRadiusSq)
		{
		  /* Simple case, we have a match to a single FIRST */
		  if (distance < closestDistance)
		    {
		      closestIndex = index;
		      closestDistance = distance;
		    }
		}
	      else
		{
		  /* Composite radio source -> save in chain? dist, flux, pos */
	        }
	    }
	}
      index++;
    }
      
  /* Do we have a simple match? */
  if (closestIndex != -1)
    {
      /* Yes. Copy the catalog info to the object. */
      obj->firstMatch = 1;
      obj->firstId    = catalog->obj[closestIndex].id;
      obj->firstLambda= catalog->obj[closestIndex].lambda;
      obj->firstEta   = catalog->obj[closestIndex].eta;
      obj->firstDelta = sqrt(closestDistance) * at_deg2Asec;
      obj->firstPeak  = catalog->obj[closestIndex].fpeak;
      obj->firstInt   = catalog->obj[closestIndex].fint;
      obj->firstrms   = catalog->obj[closestIndex].rms;
      obj->firstmajor = catalog->obj[closestIndex].maj;
      obj->firstminor = catalog->obj[closestIndex].min;
      obj->firstpa    = catalog->obj[closestIndex].pa;
       /*printf( "\nFIRST match: %8.4f %8.4f %5.3f %d\n",
	   obj->ra,obj->dec,obj->firstDelta ,obj->id);  */
    }
  else
    {
      /* No.  Set all first variables in object to 0. */
      obj->firstMatch = 0;
      obj->firstId    = 0;
      obj->firstLambda= 0;
      obj->firstEta   = 0;
      obj->firstDelta = 0;
      obj->firstPeak  = 0;
      obj->firstInt   = 0;
      obj->firstrms   = 0;
      obj->firstmajor = 0;
      obj->firstminor = 0;
      obj->firstpa    = 0;
		if(catalog->nrow >0 && catalog->ncol > 0) {
		/* check to see if we're within the catalog at all */
	atSurveyToEq(obj->lambda,obj->eta,&ra,&dec);
	xpix=(ra-catalog->crval1)/catalog->cdelt1+(catalog->crpix1-1);
	ypix=(dec-catalog->crval2)/catalog->cdelt2+(catalog->crpix2-1);
	if(xpix >= catalog->ncol || xpix < 0 || ypix >= catalog->nrow || ypix < 0) {
      		obj->firstMatch = -9999;
	} else {
	    if(catalog->covermap->rows_fl32[ypix][xpix] < 0.01) {
      		obj->firstMatch = -9999;
	     }
	}
	}
    }
  return;
}


/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	taObjMatchAgainstRosat
**
**<HTML>
**	Match a calibrated object (TA_CALIB_OBJ) against the ROSAT catalog.
**	Objects match if their separation is
**	less than or equal to the match radius.  If more than one match is
**	found, the closest match only is taken.  If a match is found, the
**	various "rosatXXX" fields are filled in the
**	TA_CALIB_OBJ object from the matching TA_ROSAT object; else those
**	fields are all set to 0.  No ROSAT-related target flags are modified
**	(i.e., this routine simply matches up, it does not select ROSAT
**	targets).
**</HTML>
**</AUTO>
******************************************************************************
******************************************************************************
*/
void
taObjMatchAgainstRosat(TA_CALIB_OBJ *obj,
		       /* MODIFIED: calibrated object to match */
		       TA_ROSAT_CATALOG *catalog
		       /* MOD: ROSAT catalog */)
{
  /* Save position where we last started the search */
  static int previousStart;

  int index, closestIndex;
  double closestDistance, lambdaLo, lambdaHi, etaDiff, lambdaDiff;
  double negMatchRadius, cosLambda, matchRadiusSq, distance;

  /* Set up useful numbers for matching */
  negMatchRadius = - catalog->matchRadius;
  cosLambda = cos(obj->lambda * at_deg2Rad);

  /* Find the object at the smallest lambda which is a possible match */
  lambdaLo = obj->lambda - catalog->matchRadius;
  if (catalog->newSearch)
    {
      previousStart = taRosatBinomialSearch(lambdaLo, catalog->obj,
					    catalog->nObj) + 1;
      catalog->newSearch = 0;
    }
  else
    {
      while (previousStart > 0)
	{
	  if (catalog->obj[previousStart-1].lambda < lambdaLo)
	    break;
	  previousStart--;
	}
      while (previousStart < catalog->nObj)
	{
	  if (catalog->obj[previousStart].lambda >= lambdaLo)
	    break;
	  previousStart++;
	}
    }

  /* Now check each object until we're beyond the match radius or we run
   * out of the catalog */
  lambdaHi = obj->lambda + catalog->matchRadius;
  matchRadiusSq = catalog->matchRadius * catalog->matchRadius;
  closestDistance = 1000. * matchRadiusSq;
  closestIndex = -1;
  index = previousStart;
  while (index < catalog->nObj)
    {
      /* Break out if we're beyond the match radius */
      if (catalog->obj[index].lambda > lambdaHi) break;

      etaDiff = obj->eta - catalog->obj[index].eta;
      if (etaDiff < -180) etaDiff += 360.;
      if (etaDiff > 180) etaDiff -= 360.;
      etaDiff *= cosLambda;
      if (etaDiff <= catalog->matchRadius && etaDiff >= negMatchRadius)
	{
	  lambdaDiff = (obj->lambda - catalog->obj[index].lambda);
	  distance = lambdaDiff * lambdaDiff + etaDiff * etaDiff;
	  if (distance <= matchRadiusSq)
	    {
	      /* This matches.  Save it if its the closest match thus far. */
	      if (distance < closestDistance)
		{
		  closestIndex = index;
		  closestDistance = distance;
		}
	    }
	}
      index++;
    }
      
  /* Do we have a match? */
  if (closestIndex != -1)
    {
      /* Yes.  Copy the catalog info to the object. */
      obj->rosatMatch   =catalog->obj[closestIndex].id;
      obj->rosatDelta   = sqrt(closestDistance) * at_deg2Asec;
      obj->rosatPosErr  = catalog->obj[closestIndex].posErr;
      obj->rosatCPS.val = catalog->obj[closestIndex].cps.val;
      obj->rosatCPS.err = catalog->obj[closestIndex].cps.err;
      obj->rosatHR1.val = catalog->obj[closestIndex].hr1.val;
      obj->rosatHR1.err = catalog->obj[closestIndex].hr1.err;
      obj->rosatHR2.val = catalog->obj[closestIndex].hr2.val;
      obj->rosatHR2.err = catalog->obj[closestIndex].hr2.err;
      obj->rosatExt     = catalog->obj[closestIndex].ext;
      obj->rosatExtLike = catalog->obj[closestIndex].extLike;
      obj->rosatDetectLike = catalog->obj[closestIndex].detectLike;
      obj->rosatExposure= catalog->obj[closestIndex].exposure;
    }
  else
    {
      /* No.  Set all first variables in object to 0. */
      obj->rosatMatch   = 0;
      obj->rosatDelta   = 0;
      obj->rosatPosErr  = 0;
      obj->rosatCPS.val = 0;
      obj->rosatCPS.err = 0;
      obj->rosatHR1.val = 0;
      obj->rosatHR1.err = 0;
      obj->rosatHR2.val = 0;
      obj->rosatHR2.err = 0;
      obj->rosatExt     = 0;
      obj->rosatExtLike = 0;
      obj->rosatDetectLike = 0;
      obj->rosatExposure= 0;
    }
  return;
}


/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	taRecalculateAstrometry
**
**<HTML>
** Recalculate position and proper motion in USNO-B catalog, including
** SDSS position.
**</HTML>
**</AUTO>
******************************************************************************
******************************************************************************
*/
void
taRecalculateAstrometry(TA_CALIB_OBJ *obj,        /* MOD: SDSS object */
			TA_USNO *usno,            /* IN: matching USNO object*/
			TA_USNO_CATALOG *usnoCat, /* IN: USNO catalog */
			double *pmRa,             /* IN: RA proper motion
						   * (mas/yr, cos(dec) incl)*/
			double *pmDec             /* IN: Dec proper motion
						   * (mas/yr) */)
{
  int i;
  double dt, X, XX, raY, raXY, decY, decXY, delta, weight, ra, dec;
  /* Error per measurement squared for SDSS and USNO-B (arcsec). */
  double sdssErr2 = 0.045 * 0.045, usnoErr2 = 0.12 * 0.12;

  /* Fitting ra = ra(t0) + dra * (t - t0), dec = dec(t0) + ddec * (t - t0)
   * where t0 = epoch of SDSS observation (usno->catEpoch).
   */

  /* Verify inputs */
  shAssert(obj != NULL);
  shAssert(usno != NULL);
  shAssert(usnoCat != NULL);
  shAssert(pmRa != NULL);
  shAssert(pmDec != NULL);

  /* Initialize sums with SDSS position */
  dt = slaEpj(obj->field->mjd[2]) - usnoCat->catEpoch;
  X = (dt / sdssErr2);
  XX = (dt * dt / sdssErr2);
  ra = obj->ra;
  dec = obj->dec;
  raY = (ra / sdssErr2);
  raXY = (dt * ra / sdssErr2);
  decY = (dec / sdssErr2);
  decXY = (dt * dec / sdssErr2);
  weight = (1. / sdssErr2);

  /* Add in USNO-B detections */
  for (i = 0; i < 5; i++)
    {
      /* Skip if not a detection */
      if (usno->field[i] == 0) continue;
      shAssert(usnoCat->epoch[usno->field[i]] > 1930. &&
	       usnoCat->epoch[usno->field[i]] < 2010.);

      /* Update sums */
      dt = usnoCat->epoch[usno->field[i]] - usnoCat->catEpoch;
      X += (dt / usnoErr2);
      XX += (dt * dt / usnoErr2);
      ra = usno->plateRa[i] / 360000.;
      dec = usno->plateSpd[i] / 360000. - 90.;
      raY += (ra / usnoErr2);
      raXY += (dt * ra / usnoErr2);
      decY += (dec / usnoErr2);
      decXY += (dt * dec / usnoErr2);
      weight += (1. / usnoErr2);
    }

  /* Now calculate the new values */
  delta = weight * XX - X * X;
  dec = (XX * decY - X * decXY) / delta;
  *pmRa = (weight * raXY - X * raY) * 3600000. * cos(dec * at_deg2Rad) / delta;
  *pmDec = (weight * decXY - X * decY) * 3600000. / delta;
}


/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	taObjMatchAgainstUSNO
**
**<HTML>
**	Match a calibrated object (TA_CALIB_OBJ) against the USNO catalog.
**	Objects match if their separation is
**	less than or equal to the match radius.  If more than one match is
**	found, the closest match only is taken.  If a match is found, the
**	various proper motions fields are filled in the
**	TA_CALIB_OBJ object from the matching TA_USNO object; else those
**	fields are all set to 0.  No proper-motion-related target flags are
**	modified (i.e., this routine simply matches up, it does not select
**	targets).
**</HTML>
**</AUTO>
******************************************************************************
******************************************************************************
*/
void
taObjMatchAgainstUSNO(TA_CALIB_OBJ *obj,
		      /* MODIFIED: calibrated object to match */
		      TA_USNO_CATALOG *catalog
		      /* MOD: USNO catalog */)
{
  /* Save position where we last started the search */
  static int previousStart;

  int index, closestIndex, field;
  double closestDistance, lambdaLo, lambdaHi, raDiff, decDiff;
  double negMatchRadius, cosDec, matchRadiusSq, distance;
  double catRa, catDec, pmRa, pmDec;
  TA_USNO *usno;

  /* Set up useful numbers for matching */
  negMatchRadius = - catalog->matchRadius;
  cosDec = cos(obj->dec * at_deg2Rad);

  /* Find the object at the smallest lambda which is a possible match */
  lambdaLo = (obj->lambda - catalog->matchRadius) * 360000.;
  if (catalog->newSearch)
    {
      previousStart = taUSNOBinomialSearch(lambdaLo, catalog->obj,
					   catalog->nObj) + 1;
      catalog->newSearch = 0;
    }
  else
    {
      while (previousStart > 0)
	{
	  if (catalog->obj[previousStart-1].lambda < lambdaLo)
	    break;
	  previousStart--;
	}
      while (previousStart < catalog->nObj)
	{
	  if (catalog->obj[previousStart].lambda >= lambdaLo)
	    break;
	  previousStart++;
	}
    }

  /* Now check each object until we're beyond the match radius or we run
   * out of the catalog */
  lambdaHi = (obj->lambda + catalog->matchRadius) * 360000.;
  matchRadiusSq = catalog->matchRadius * catalog->matchRadius;
  closestDistance = 1000. * matchRadiusSq;
  closestIndex = -1;
  index = previousStart;
  while (index < catalog->nObj)
    {
      /* Break out if we're beyond the match radius */
      if (catalog->obj[index].lambda > lambdaHi) break;

      catDec = catalog->obj[index].spd / 360000. - 90.;
      decDiff = obj->dec - catDec;
      if (decDiff <= catalog->matchRadius && decDiff >= negMatchRadius)
	{
	  catRa = catalog->obj[index].ra / 360000.;
	  raDiff = obj->ra - catRa;
	  if (raDiff < -180) raDiff += 360.;
	  if (raDiff > 180) raDiff -= 360.;
	  raDiff *= cosDec;
	  if (raDiff <= catalog->matchRadius && raDiff >= negMatchRadius)
	    {
	      distance = raDiff * raDiff + decDiff * decDiff;
	      if (distance <= matchRadiusSq)
		{
		  /* This matches.  Save if its the closest match thus far. */
		  if (distance < closestDistance)
		    {
		      closestIndex = index;
		      closestDistance = distance;
		    }
		}
	    }
	}
      index++;
    }
      
  /* Do we have a match? */
  if (closestIndex != -1)
    {
      /* Yes. */
      usno = &(catalog->obj[closestIndex]);
      obj->properMotionMatch = 1;
      obj->properMotionDelta = sqrt(closestDistance) * at_deg2Asec;
      taRecalculateAstrometry(obj, usno, catalog, &pmRa, &pmDec);
      obj->properMotion = sqrt(pmRa * pmRa + pmDec * pmDec) / 10.;
      obj->properMotionAngle = slaDranrm(atan2(pmRa, pmDec)) * at_rad2Deg;
      obj->usnoBlue = ((usno->mag[2] > 0 && usno->mag[2] < 2500) ?
		       usno->mag[2] : usno->mag[0]) / 100.;
      obj->usnoRed = ((usno->mag[3] > 0 && usno->mag[3] < 2500) ?
		       usno->mag[3] : usno->mag[1]) / 100.;

    }
  else
    {
      /* No.  Set all USNO variables in object to 0. */
      obj->properMotionMatch = 0;
      obj->properMotionDelta = 0.;
      obj->properMotion = 0.;
      obj->properMotionAngle = 0.;
      obj->usnoBlue = 0.;
      obj->usnoRed = 0.;
    }
  return;
}