/*****************************************************************************
******************************************************************************
**
**
** FILE: catalogs.c
**
**
** This file contains C functions to match objects against other catalogs.
**
**
**
**
** ENVIRONMENT:
** ANSI C.
**
******************************************************************************
******************************************************************************
*/
#include
#include
#include
#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;
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taObjMatchAgainstFirst
**
**
** 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).
**
**
******************************************************************************
****************************************************************************** */
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;
/* 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;
}
return;
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taObjMatchAgainstRosat
**
**
** 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).
**
**
******************************************************************************
******************************************************************************
*/
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;
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taObjMatchAgainstUSNO
**
**
** 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).
**
**
******************************************************************************
******************************************************************************
*/
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, closestAngle;
double negMatchRadius, cosDec, matchRadiusSq, distance, dEpoch;
double epj, catRa, catDec;
/* 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;
closestAngle = 90. * at_deg2Rad -
slaDbear(catRa * at_deg2Rad,
catDec * at_deg2Rad,
obj->ra * at_deg2Rad,
obj->dec * at_deg2Rad);
closestAngle = slaDranrm(closestAngle) * at_rad2Deg;
}
}
}
}
index++;
}
/* Do we have a match? */
if (closestIndex != -1)
{
/* Yes. */
field = catalog->obj[closestIndex].mags / 1000000;
epj = 2000. + (obj->field->mjd[2] - 51544.) / 365.25;
dEpoch = (epj - catalog->epoch[field]) / 100.;
obj->properMotionMatch = 1;
obj->properMotionDelta = sqrt(closestDistance) * at_deg2Asec;
obj->properMotion = obj->properMotionDelta / dEpoch;
obj->properMotionAngle = closestAngle;
obj->usnoBlue = ((catalog->obj[closestIndex].mags / 1000) % 1000) / 10.0;
obj->usnoRed = (catalog->obj[closestIndex].mags % 1000) / 10.0;
}
else
{
/* No. Set all first variables in object to 0. */
obj->properMotionMatch = 0;
obj->properMotionDelta = 0.;
obj->properMotion = 0.;
obj->properMotionAngle = 0.;
obj->usnoBlue = 0.;
obj->usnoRed = 0.;
}
return;
}