#include "taCalibObj.h"
#include "arDetectionFlag.h"
#include "taCosmic.h"
/* Error = -9999 represents an undefined value and error, error = -1000
* represents undefined error (but value is defined), and for values that
* lack an associated error, -9999 represents an undefined value. */
#define taUndefined(x) (fabs((x) + 9999.) < 0.01 ? 1 : 0)
#define taErrUndefined(x) (fabs((x) + 1000.) < 0.01 ? 1 : 0)
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taGood
**
**
** Determine whether an object is "GOOD" (i.e., it is considered to be
** a good measured object, and is thus eligible as either a primary or
** secondary detection in the survey) based on the values of its
** "objc_flags". The "objc_flags" are passed in as an integer. The
** routine returns 1 if the object is good, 0 if it is bad.
**
** An object is considered good if it is not deblended
** and not a BRIGHT object. The elimination of parents with children
** and bright objects cuts out multiple versions of the same object
** before matching/resolving the edges of scan lines.
**
**
******************************************************************************
******************************************************************************
*/
int
taGood(int flags /* IN: object's "objc_flags" */)
{
return ! (flags & AR_DFLAG_BRIGHT) &&
(! (flags & AR_DFLAG_BLENDED) ||
((flags & AR_DFLAG_BLENDED) && (flags & AR_DFLAG_NODEBLEND)));
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taTransCosmicSet
**
**
** Set the cosmic magnitudes and scatters used by the astrotools TRANS
** routines, based on the passed in cosmic colors and scatters.
** This is unfortunately convoluted, as the TRANS functions want cosmic
** magnitudes, but we think in terms of cosmic colors. This routine
** makes the necessary transformations.
**
** NOTE: this code works only if the scatter decreases along the
** cosmicScatter array (that is, cosmicScatter[i] > cosmicScatter[i+1]).
**
**
******************************************************************************
******************************************************************************
*/
void
taTransCosmicSet(const double cosmicColor[4], /* IN: cosmic colors, in order
* u-g, g-r, r-i, i-z */
const double cosmicScatter[4] /* IN: cosmic scatter, in order
* u-g, g-r, r-i, i-z */)
{
float cosmicMag[5], cosmicMagScatter[5];
int i;
/* Start by setting z = 0 and zErr = 0, and work backwards.
* NOTE: this code works only if the scatter decreases along the
* cosmicScatter array (that is, cosmicScatter[i] > cosmicScatter[i+1]). */
cosmicMag[4] = 0.;
cosmicMagScatter[4] = 0.;
for (i = 3; i >= 0; i--)
{
cosmicMag[i] = cosmicColor[i] + cosmicMag[i+1]; /* e.g., i=(i-z)+z */
cosmicMagScatter[i] = sqrt(cosmicScatter[i]*cosmicScatter[i] -
cosmicMagScatter[i+1]*cosmicMagScatter[i+1]);
}
atCosmicMagSet(cosmicMag);
atCosmicMagScatterSet(cosmicMagScatter);
return;
}
/*****************************************************************************
******************************************************************************
**
** taLCaliApply --- Handle bad magnitude measurements.
**
******************************************************************************
******************************************************************************
*/
int taLCaliApply
(
double counts, /* IN: counts in the reference filter */
double countsErr, /* IN: error in counts */
double cCounts, /* IN: count in the adjacent filter */
double cCountsErr, /* IN: error in cCounts */
double intTime, /* IN: exposure time in the reference filter */
double cIntTime, /* IN: exposure time in the adjacent filter */
double airmass, /* IN: airmass (sec(Z)) */
double zeropt, /* IN: zeropoint (a') */
double extinct, /* IN: extinction prime (k') */
double colorTerm, /* IN: color term (b') */
double secColorTerm, /* IN: second order color term (c') */
int sign, /* IN: usually +1 for u g r i; -1 for z */
double cosmicColor, /* IN: cosmic reference-adjacent color */
double cosmicError, /* IN: scatter in cosmicColor */
double *calMag, /* OUT: calibrated luptitude */
double *calMagErr, /* OUT: error in calMag */
int *status /* OUT: status */
)
{
int returnValue = 0;
/* If counts undefined, return -9999 indicating undefined results. */
if (taUndefined(countsErr))
{
*calMag = -9999.;
*calMagErr = -9999.;
return returnValue;
}
/* If adjacent filter counts and/or error are undefined, set adjacent
* counts negative which will force atLCaliApply to use the cosmic color. */
if (taUndefined(cCountsErr) || taErrUndefined(cCountsErr)) cCounts = -1000.;
/* Calibrate */
returnValue = atLCaliApply(counts, countsErr, cCounts, cCountsErr, intTime,
cIntTime, airmass, zeropt, extinct, colorTerm,
secColorTerm, sign, cosmicColor, cosmicError,
calMag, calMagErr, status);
/* If counts error undefined, return -1000 indicating the error only is
* undefined. */
if (taErrUndefined(countsErr)) *calMagErr = -1000.;
return returnValue;
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taObjCalibrate
**
**
** Calibrate a single TA_CALIB_OBJ. The calibration is done
** in place. Thus, the TA_CALIB_OBJ should contain instrumental
** quantities on input, and will contain calibrated quantities on output.
** The argument "trans" is a pointer to a 5-element array of pointers to
** TRANS structures, in the order u, g, r, i, z. These are the current
** 'best' astrometric transformations. The argument "transPhoto" is
** the TRANS structures used when photo was run, corrected using
** photo's "rowOffset" and "colOffset", and are required to
** reproduce the position offsets in each filter.
** For DCR corrections in astrometric transformations, PSF magnitudes are
** used for stars and model magnitudes for galaxies.
**
** If "raw" is set to 1, then lengths, magnitudes, surface brightnesses
** and angles are not calibrated. Celestial coordinates and extinction
** values are still calculated, however.
**
**
******************************************************************************
******************************************************************************
*/
RET_CODE
taObjCalibrate(TA_CALIB_OBJ *obj, /* MOD: Object to calibrate */
const TA_PTRANS *ptrans,/* IN: Photometric calibration */
const TRANS *trans[5], /* IN: Current 'best' astrometric
* calibration (ugriz) */
const TRANS *transPhoto[5], /* IN: astrometric calibration
* used by photo, corrected
* for photo's determined
* rowOffset and colOffset
* (ugriz)*/
double node, /* IN: GC ascending node (degrees) */
double incl /* IN: GC inclination (degrees) */,
char refFilter, /* IN: Ref. filter (primary position)*/
double expTime /* IN: Exposure time (seconds) */,
int raw /* IN: 1=don't calibrate lengths,
* magnitudes, surface brightesses,
* and angles; 0=do */)
{
TA_DETECTION *det = NULL, *adjDet = NULL;
int m, p;
double scale, scale2, mag, magErr, mu, nu, ra, dec, cosDec, dust, ra2, dec2;
double pa, refRa, refDec;
double psf[TA_NFILTERS], psfErr[TA_NFILTERS];
double model[TA_NFILTERS], modelErr[TA_NFILTERS];
float transMag[TA_NFILTERS], transMagErr[TA_NFILTERS];
/* Filter order in array fields */
const char *filters = "ugriz";
/*adjacent filter to use in colorterms*/
const int adj[TA_NFILTERS] = {1, 2, 3, 4, 3};
/* sign to use for color terms */
const int sign[TA_NFILTERS] = {1, 1, 1, 1, -1};
int status, refIndex = -1, nProf, nOtherProf, c;
/* Error = -9999 represents an undefined value and error, error = -1000
* represents undefined error (but value is defined), and for values that
* lack an associated error, -9999 represents an undefined value. Don't
* calibrate undefined values. */
/* For each filter ... */
for (m = 0; m < TA_NFILTERS; m++)
{
/* Get pointers to the detection, as well as the detection in the
* adjacent filter used to create the color used in the photometric
* calibrations and the cosmic color index. */
det = &(obj->detection[m]);
adjDet = &(obj->detection[adj[m]]);
c = cosmicIdx[m];
if (raw == 0)
{
/* Scale calibrations (pixels -> arcsecs) */
scale = at_deg2Asec * sqrt((double)(trans[m]->b * trans[m]->f -
trans[m]->e * trans[m]->c));
if (! taUndefined(det->petroRad.err))
{
det->petroRad.val *= scale;
if (! taErrUndefined(det->petroRad.err))
det->petroRad.err *= scale;
}
if (! taUndefined(det->petroR50.err))
{
det->petroR50.val *= scale;
if (! taErrUndefined(det->petroR50.err))
det->petroR50.err *= scale;
}
if (! taUndefined(det->petroR90.err))
{
det->petroR90.val *= scale;
if (! taErrUndefined(det->petroR90.err))
det->petroR90.err *= scale;
}
if (! taUndefined(det->r_deV.err))
{
det->r_deV.val *= scale;
if (! taErrUndefined(det->r_deV.err)) det->r_deV.err *= scale;
}
if (! taUndefined(det->r_exp.err))
{
det->r_exp.val *= scale;
if (! taErrUndefined(det->r_exp.err)) det->r_exp.err *= scale;
}
/* Isophotal values */
if (! taUndefined(det->iso_rowc.err))
det->iso_rowc.val -= det->iso_rowc.grad * ptrans->dZero[m];
if (! taUndefined(det->iso_colc.err))
det->iso_colc.val -= det->iso_colc.grad * ptrans->dZero[m];
if (! taUndefined(det->iso_a.err))
{
det->iso_a.grad *= scale;
det->iso_a.val = det->iso_a.val * scale
- det->iso_a.grad * ptrans->dZero[m];
if (! taErrUndefined(det->iso_a.err)) det->iso_a.err *= scale;
}
if (! taUndefined(det->iso_b.err))
{
det->iso_b.grad *= scale;
det->iso_b.val = det->iso_b.val * scale
- det->iso_b.grad * ptrans->dZero[m];
if (! taErrUndefined(det->iso_b.err)) det->iso_b.err *= scale;
}
if (! taUndefined(det->iso_phi.err))
det->iso_phi.val -= det->iso_phi.grad * ptrans->dZero[m];
/* Flux calibrations (counts -> luptitudes) */
taLCaliApply(det->psfCounts.val, det->psfCounts.err,
adjDet->psfCounts.val, adjDet->psfCounts.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->psfCounts.val = mag;
det->psfCounts.err = magErr;
taLCaliApply(det->fiberCounts.val, det->fiberCounts.err,
adjDet->fiberCounts.val, adjDet->fiberCounts.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->fiberCounts.val = mag;
det->fiberCounts.err = magErr;
taLCaliApply(det->petroCounts.val, det->petroCounts.err,
adjDet->petroCounts.val, adjDet->petroCounts.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->petroCounts.val = mag;
det->petroCounts.err = magErr;
taLCaliApply(det->counts_deV.val, det->counts_deV.err,
adjDet->counts_deV.val, adjDet->counts_deV.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->counts_deV.val = mag;
det->counts_deV.err = magErr;
taLCaliApply(det->counts_exp.val, det->counts_exp.err,
adjDet->counts_exp.val, adjDet->counts_exp.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->counts_exp.val = mag;
det->counts_exp.err = magErr;
taLCaliApply(det->counts_model.val, det->counts_model.err,
adjDet->counts_model.val, adjDet->counts_model.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
det->counts_model.val = mag;
det->counts_model.err = magErr;
/* Surface brightnesses (counts/pixel^2 -> luptitudes/arcsecs^2) */
scale2 = scale * scale;
taLCaliApply(det->sky.val / scale2, det->sky.err / scale2,
adjDet->sky.val / scale2, adjDet->sky.err / scale2,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicSkyColor[c], cosmicSkyScatter[c],
&mag, &magErr, &status);
det->sky.val = mag;
det->sky.err = magErr;
nProf = det->nProf;
nOtherProf = adjDet->nProf;
for (p = 0; p < nProf; p++)
{
if (p < nOtherProf)
{
taLCaliApply(det->profile[p].val / scale2,
det->profile[p].err / scale2,
adjDet->profile[p].val / scale2,
adjDet->profile[p].err / scale2,
expTime, expTime, trans[m]->airmass,
ptrans->a[m].val, ptrans->k[m].val,
ptrans->b[m].val, ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c], &mag,
&magErr, &status);
det->profile[p].val = mag;
det->profile[p].err = magErr;
}
else
{
/* Radial bin doesn't exist for other filter. Force use of
* cosmic color by setting other filter's mag = 0,
* err = -1000.*/
taLCaliApply(det->profile[p].val / scale2,
det->profile[p].err / scale2,
0., -1000.,
expTime, expTime, trans[m]->airmass,
ptrans->a[m].val, ptrans->k[m].val,
ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c], &mag,
&magErr, &status);
det->profile[p].val = mag;
det->profile[p].err = magErr;
}
}
/* Save following magnitudes for use in astrometric calibrations */
psf[m] = det->psfCounts.val;
psfErr[m] = det->psfCounts.err;
model[m] = det->counts_model.val;
modelErr[m] = det->counts_model.err;
}
else
{
/* Must calculate magnitudes for use in astrometric calibrations */
taLCaliApply(det->psfCounts.val, det->psfCounts.err,
adjDet->psfCounts.val, adjDet->psfCounts.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
psf[m] = mag;
psfErr[m] = magErr;
taLCaliApply(det->counts_model.val, det->counts_model.err,
adjDet->counts_model.val, adjDet->counts_model.err,
expTime, expTime, trans[m]->airmass, ptrans->a[m].val,
ptrans->k[m].val, ptrans->b[m].val,ptrans->c[m].val,
sign[m], cosmicColor[c], cosmicScatter[c],
&mag, &magErr, &status);
model[m] = mag;
modelErr[m] = magErr;
}
}
/* Now positions. Use PSF magnitudes for DCR corrections for stars, model
* magnitudes for galaxies.. */
for (m = 0; m < TA_NFILTERS; m++)
{
if (obj->objc_type == AR_OBJECT_TYPE_STAR)
{
transMag[m] = psf[m];
transMagErr[m] = psfErr[m];
}
else
{
transMag[m] = model[m];
transMagErr[m] = modelErr[m];
}
if (refFilter == filters[m]) refIndex = m;
}
if (refIndex == -1)
{
shErrStackPush("taObjCalibrate: illegal refFilter");
return SH_GENERIC_ERROR;
}
atTransApply(trans[refIndex], refFilter, obj->objc_rowc.val, 0.,
obj->objc_colc.val, 0., transMag, transMagErr, &mu, NULL,
&nu, NULL);
atGCToEq(mu, nu, &(obj->ra), &(obj->dec), node, incl);
atEqToSurvey(obj->ra, obj->dec, &(obj->lambda), &(obj->eta));
atEqToGal(obj->ra, obj->dec, &(obj->l), &(obj->b));
/* Photo's position in reference filter */
atTransApply(transPhoto[refIndex], refFilter,
obj->detection[refIndex].rowc.val, 0.,
obj->detection[refIndex].colc.val, 0., transMag, transMagErr,
&mu, NULL, &nu, NULL);
atGCToEq(mu, nu, &refRa, &refDec, node, incl);
cosDec = cos(refDec * at_deg2Rad);
for (m = 0; m < TA_NFILTERS; m++)
{
/* Photo measured position offsets in each filter */
if (m == refIndex)
{
obj->offsetRa[m] = 0.;
obj->offsetDec[m] = 0.;
}
else
{
atTransApply(transPhoto[m], filters[m], obj->detection[m].rowc.val,
0., obj->detection[m].colc.val, 0., transMag,
transMagErr, &mu, NULL, &nu, NULL);
atGCToEq(mu, nu, &ra, &dec, node, incl);
obj->offsetRa[m] = (ra - refRa) * cosDec * at_deg2Asec;
obj->offsetDec[m] = (dec - refDec) * at_deg2Asec;
}
/* Position angles. */
if (raw == 0)
{
atTransApply(trans[m], filters[m], obj->detection[m].rowc.val, 0.,
obj->detection[m].colc.val, 0., transMag, transMagErr,
&mu, NULL, &nu, NULL);
atGCToEq(mu, nu, &ra, &dec, node, incl);
atTransApply(trans[m], filters[m],
obj->detection[m].rowc.val + 100., 0.,
obj->detection[m].colc.val, 0., transMag, transMagErr,
&mu, NULL, &nu, NULL);
atGCToEq(mu, nu, &ra2, &dec2, node, incl);
pa = 90. - slaDbear(ra * at_deg2Rad, dec * at_deg2Rad,
ra2 * at_deg2Rad, dec2 * at_deg2Rad) *at_rad2Deg;
if (! taUndefined(obj->detection[m].iso_phi.err))
obj->detection[m].iso_phi.val =
slaDranrm((obj->detection[m].iso_phi.val+pa)*at_deg2Rad) *
at_rad2Deg;
if (! taUndefined(obj->detection[m].phi_deV.err))
obj->detection[m].phi_deV.val =
slaDranrm((obj->detection[m].phi_deV.val+pa)*at_deg2Rad) *
at_rad2Deg;
if (! taUndefined(obj->detection[m].phi_exp.err))
obj->detection[m].phi_exp.val =
slaDranrm((obj->detection[m].phi_exp.val+pa)*at_deg2Rad) *
at_rad2Deg;
}
}
/* Get the extinction */
if (taExtinction(obj->l, obj->b, &(obj->reddening[0]), &(obj->reddening[1]),
&(obj->reddening[2]), &(obj->reddening[3]),
&(obj->reddening[4])) != SH_SUCCESS)
return SH_GENERIC_ERROR;
/* Succesful return */
return SH_SUCCESS;
}
/*****************************************************************************
******************************************************************************
**
**
** ROUTINE: taExtinction
**
**
** Return the extinction in each filter at the specified galactic
** coordinates.
**
**
******************************************************************************
******************************************************************************
*/
RET_CODE
taExtinction(double l, /* IN: Galactic longitude (degrees) */
double b, /* IN: Galactic latitude (degrees) */
float *u, /* OUT: extinction in u band */
float *g, /* OUT: extinction in g band */
float *r, /* OUT: extinction in r band */
float *i, /* OUT: extinction in i band */
float *z /* OUT: extinction in z band */
)
{
double ebv;
/* HARDWIRE: hardwiring the following conversions from E(B-V) to
* extinction (Rv factors). */
double uRv = 5.155, gRv = 3.793, rRv = 2.751, iRv = 2.086, zRv = 1.479;
/* Get the reddening from dust maps. */
if (atDustGetvalFast(l, b, "Ebv", 1, 0, &ebv) != SH_SUCCESS)
return SH_GENERIC_ERROR;
/* Return the extinctions */
*u = uRv * ebv;
*g = gRv * ebv;
*r = rRv * ebv;
*i = iRv * ebv;
*z = zRv * ebv;
return SH_SUCCESS;
}