/***************************************************************************** ****************************************************************************** ** ** ** FILE: calibrate.c ** ** ** This file contains C functions to apply photometric and astrometric ** calibrations. ** ** ** ** ** ENVIRONMENT: ** ANSI C. ** ****************************************************************************** ****************************************************************************** */ #include #include #include #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. ** ** ****************************************************************************** ****************************************************************************** */ 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) */) { 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; 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]; /* 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; } } } /* 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] = obj->detection[m].psfCounts.val; transMagErr[m] = obj->detection[m].psfCounts.err; } else { transMag[m] = obj->detection[m].counts_model.val; transMagErr[m] = obj->detection[m].counts_model.err; } 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. */ 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; }