Changeset 32406 for trunk/ppTranslate/src/ppMopsRead.c
- Timestamp:
- Sep 15, 2011, 4:50:36 PM (15 years ago)
- File:
-
- 1 edited
-
trunk/ppTranslate/src/ppMopsRead.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/ppTranslate/src/ppMopsRead.c
r32396 r32406 20 20 psArray *detections = psArrayAlloc(num); // Array of detections, to return 21 21 for (int i = 0; i < num; i++) { 22 psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file 23 22 const char *name = inNames->data[i]; 23 24 psFits *fits = psFitsOpen(name, "r"); // FITS file 24 25 if (!fits) { 25 26 psError(PS_ERR_IO, false, "Unable to open input %d", i); 26 27 return false; 27 28 } 29 28 30 psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header 29 31 if (!header) { … … 74 76 } 75 77 ppMopsDetections *det = ppMopsDetectionsAlloc(size); 78 detections->data[i] = det; 79 det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension 80 det->num = size; 81 det->diffSkyfileId = diffSkyfileId; 76 82 77 83 psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); … … 90 96 psMetadataLookupF32(NULL, header, "FWHM_MIN")); 91 97 92 intnaxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns93 intnaxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows98 det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns 99 det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows 94 100 95 101 psFree(header); 96 102 97 #ifndef USE_OLD_READ_TABLE 98 // use new fits table functions 99 100 psFitsTable *table = psFitsReadTableNew(fits); // Table of interest 103 psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest 101 104 if (!table) { 102 105 psError(PS_ERR_IO, false, "Unable to read table %d", i); 103 return false; 104 } 106 return NULL; 107 } 108 det->table = table; 105 109 psFitsClose(fits); 110 if (args->version == 0) { 111 if (skyChipPsfVersion < 2) { 112 // XXX: TODO: Do we need to add dummy vectors for the missing columns? 113 } 114 } 115 116 psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); 117 psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); 118 119 det->raErr = psVectorAlloc(size, PS_TYPE_F64); 120 det->decErr = psVectorAlloc(size, PS_TYPE_F64); 121 det->mask = psVectorAlloc(size, PS_TYPE_U8); 122 123 // convert ra and dec to radians for use in the purge duplicates function 124 det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); 125 det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); 126 det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); 127 det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); 128 if (!det->ra || !det->dec || !det->x || !det->y) { 129 psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); 130 return NULL; 131 } 132 133 // Add our new vectors to the table so that duplicates and masked items may be purged 134 psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr); 135 psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr); 136 137 psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); 138 139 psVector *mag = psMetadataLookupVector(NULL, table, "PSF_INST_MAG"); 140 psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG"); 141 psVector *xErrV = psMetadataLookupVector(NULL, table, "X_PSF_SIG"); 142 psVector *yErrV = psMetadataLookupVector(NULL, table, "Y_PSF_SIG"); 143 psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE"); 144 psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE"); 145 psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS"); 106 146 107 147 double plateScale = 0.0; // Plate scale … … 109 149 for (long row = 0; row < size; row++) { 110 150 111 psU32 flags = psFitsTableGetU32(NULL, table, row, "FLAGS");151 psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS"); 112 152 if (flags & SOURCE_MASK) { 113 153 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags); 154 det->mask->data.U8[row] = 0xFF; 114 155 continue; 115 156 } 116 157 117 det->x->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "X_PSF");118 det->y->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "Y_PSF");119 det->ra->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "RA_PSF"));120 det->dec->data.F64[numGood] = DEG_TO_RAD(psFitsTableGetF64(NULL, table, row, "DEC_PSF"));121 det->mag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG");122 det->magErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_MAG_SIG");123 det->chi2->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_CHISQ");124 det->dof->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NDOF");125 det->cr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CR_NSIGMA");126 det->extended->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "EXT_NSIGMA");127 det->psfMajor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MAJOR");128 det->psfMinor->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_MINOR");129 det->psfTheta->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_THETA");130 det->quality->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF");131 det->numPix->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "PSF_NPIX");132 det->xxMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XX");133 det->xyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_XY");134 det->yyMoment->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_YY");135 det->flags->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "FLAGS");136 det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;137 det->naxis1->data.S32[numGood] = naxis1;138 det->naxis2->data.S32[numGood] = naxis2;139 det->nPos->data.S32[numGood] = psFitsTableGetS32(NULL, table, row, "DIFF_NPOS");140 det->fPos->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_FRATIO");141 det->nRatioBad->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_BAD");142 det->nRatioMask->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_MASK");143 det->nRatioAll->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_NRATIO_ALL");144 145 //Additions of 2010-10-25146 if (args->version == 2) {147 //Values are set only if the version is 2148 if (skyChipPsfVersion == 2) {149 det->psfInstFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX");150 det->psfInstFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_INST_FLUX_SIG");151 det->apMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG");152 det->apMagRaw->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RAW");153 det->apMagRadius->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_MAG_RADIUS");154 det->apFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX");155 det->apFluxSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "AP_FLUX_SIG");156 det->peakFluxAsMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PEAK_FLUX_AS_MAG");157 det->calPsfMag->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG");158 det->calPsfMagSig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "CAL_PSF_MAG_SIG");159 det->sky->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY");160 det->skySig->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "SKY_SIGMA");161 det->qualityPerfect->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "PSF_QF_PERFECT");162 det->momentsR1->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_R1");163 det->momentsRH->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "MOMENTS_RH");164 det->kronFlux->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX");165 det->kronFluxErr->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_ERR");166 det->kronFluxInner->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_INNER");167 det->kronFluxOuter->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "KRON_FLUX_OUTER");168 det->diffRP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_P");169 det->diffSnP->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_P");170 det->diffRM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_R_M");171 det->diffSnM->data.F32[numGood] = psFitsTableGetF32(NULL, table, row, "DIFF_SN_M");172 det->flags2->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "FLAGS2");173 det->ippIdet->data.U32[numGood] = psFitsTableGetU32(NULL, table, row, "IPP_IDET");174 det->nFrames->data.U16[numGood] = psFitsTableGetU16(NULL, table, row, "N_FRAMES");175 det->padding->data.S16[numGood] = psFitsTableGetS16(NULL, table, row, "PADDING");176 } else {177 det->psfInstFlux->data.F32[numGood] = NAN;178 det->psfInstFluxSig->data.F32[numGood] = NAN;179 det->apMag->data.F32[numGood] = NAN;180 det->apMagRaw->data.F32[numGood] = NAN;181 det->apMagRadius->data.F32[numGood] = NAN;182 det->apFlux->data.F32[numGood] = NAN;183 det->apFluxSig->data.F32[numGood] = NAN;184 det->peakFluxAsMag->data.F32[numGood] = NAN;185 det->calPsfMag->data.F32[numGood] = NAN;186 det->calPsfMagSig->data.F32[numGood] = NAN;187 det->sky->data.F32[numGood] = NAN;188 det->skySig->data.F32[numGood] = NAN;189 det->qualityPerfect->data.F32[numGood] = NAN;190 det->momentsR1->data.F32[numGood] = NAN;191 det->momentsRH->data.F32[numGood] = NAN;192 det->kronFlux->data.F32[numGood] = NAN;193 det->kronFluxErr->data.F32[numGood] = NAN;194 det->kronFluxInner->data.F32[numGood] = NAN;195 det->kronFluxOuter->data.F32[numGood] = NAN;196 det->diffRP->data.F32[numGood] = NAN;197 det->diffSnP->data.F32[numGood] = NAN;198 det->diffRM->data.F32[numGood] = NAN;199 det->diffSnM->data.F32[numGood] = NAN;200 det->flags2->data.U32[numGood] = 0;201 det->ippIdet->data.U32[numGood] = 0;202 det->nFrames->data.U16[numGood] = 0;203 det->padding->data.S16[numGood] = 0;204 }205 }206 207 158 // Calculate error in RA, Dec 208 double xErr = psFitsTableGetF64(NULL, table, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32 209 double yErr = psFitsTableGetF64(NULL, table, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32 210 double scale = psFitsTableGetF64(NULL, table, row, "PLTSCALE"); //SC: Warning! Promotion of F32 211 double angle = psFitsTableGetF64(NULL, table, row, "POSANGLE"); //SC: Warning! Promotion of F32 212 213 if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || 214 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || 215 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || 159 160 double xErr = xErrV->data.F32[row]; 161 double yErr = yErrV->data.F32[row]; 162 double scale = scaleV->data.F32[row]; 163 double angle = angleV->data.F32[row]; 164 165 if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) || 166 !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) || 167 !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) || 216 168 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { 217 169 psTrace("ppMops.read", 10, … … 219 171 "%f %f %lf %lf %f %f %f %f %f %f", 220 172 row, i, 221 det->x->data.F32[ numGood], det->y->data.F32[numGood],222 det->ra->data.F64[ numGood], det->dec->data.F64[numGood],223 det->mag->data.F32[numGood], det->magErr->data.F32[numGood],173 det->x->data.F32[row], det->y->data.F32[row], 174 det->ra->data.F64[row], det->dec->data.F64[row], 175 mag->data.F32[row], magErr->data.F32[row], 224 176 xErr, yErr, scale, angle); 177 det->mask->data.U8[row] = 0xFF; 225 178 continue; 226 179 } … … 231 184 double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); 232 185 double errScale = scale / 3600.0; 233 det->raErr->data.F64[ numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);234 det->decErr->data.F64[ numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);235 236 det->mask->data.U8[ numGood] = 0;186 det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); 187 det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); 188 189 det->mask->data.U8[row] = 0; 237 190 plateScale += scale; 238 191 numGood++; 239 192 } 240 #else241 // OLD way of reading fits tables read it into an array of metadata objects, one for each row. Uses lots and lots242 // and lots of memory.243 psArray *table = psFitsReadTable(fits); // Table of interest244 if (!table) {245 psError(PS_ERR_IO, false, "Unable to read table %d", i);246 return false;247 }248 psFitsClose(fits);249 250 double plateScale = 0.0; // Plate scale251 long numGood = 0; // Number of good rows252 for (long j = 0; j < size; j++) {253 psMetadata *row = table->data[j]; // Row of interest254 255 psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");256 if (flags & SOURCE_MASK) {257 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);258 continue;259 }260 261 det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");262 det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");263 det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));264 det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));265 det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");266 det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");267 det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ");268 det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF");269 det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA");270 det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");271 det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR");272 det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR");273 det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA");274 det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF");275 det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX");276 det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX");277 det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY");278 det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY");279 det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");280 det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;281 det->naxis1->data.S32[numGood] = naxis1;282 det->naxis2->data.S32[numGood] = naxis2;283 det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS");284 det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO");285 det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD");286 det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK");287 det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL");288 289 //Additions of 2010-10-25290 if (args->version == 2) {291 //Values are set only if the version is 2292 if (skyChipPsfVersion == 2) {293 det->psfInstFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX");294 det->psfInstFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_FLUX_SIG");295 det->apMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG");296 det->apMagRaw->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RAW");297 det->apMagRadius->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_MAG_RADIUS");298 det->apFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX");299 det->apFluxSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "AP_FLUX_SIG");300 det->peakFluxAsMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PEAK_FLUX_AS_MAG");301 det->calPsfMag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG");302 det->calPsfMagSig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CAL_PSF_MAG_SIG");303 det->sky->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY");304 det->skySig->data.F32[numGood] = psMetadataLookupF32(NULL, row, "SKY_SIGMA");305 det->qualityPerfect->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF_PERFECT");306 det->momentsR1->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_R1");307 det->momentsRH->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_RH");308 det->kronFlux->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX");309 det->kronFluxErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_ERR");310 det->kronFluxInner->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_INNER");311 det->kronFluxOuter->data.F32[numGood] = psMetadataLookupF32(NULL, row, "KRON_FLUX_OUTER");312 det->diffRP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_P");313 det->diffSnP->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_P");314 det->diffRM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_R_M");315 det->diffSnM->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_SN_M");316 det->flags2->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS2");317 det->ippIdet->data.U32[numGood] = psMetadataLookupU32(NULL, row, "IPP_IDET");318 det->nFrames->data.U16[numGood] = psMetadataLookupU16(NULL, row, "N_FRAMES");319 det->padding->data.S16[numGood] = psMetadataLookupS16(NULL, row, "PADDING");320 } else {321 det->psfInstFlux->data.F32[numGood] = NAN;322 det->psfInstFluxSig->data.F32[numGood] = NAN;323 det->apMag->data.F32[numGood] = NAN;324 det->apMagRaw->data.F32[numGood] = NAN;325 det->apMagRadius->data.F32[numGood] = NAN;326 det->apFlux->data.F32[numGood] = NAN;327 det->apFluxSig->data.F32[numGood] = NAN;328 det->peakFluxAsMag->data.F32[numGood] = NAN;329 det->calPsfMag->data.F32[numGood] = NAN;330 det->calPsfMagSig->data.F32[numGood] = NAN;331 det->sky->data.F32[numGood] = NAN;332 det->skySig->data.F32[numGood] = NAN;333 det->qualityPerfect->data.F32[numGood] = NAN;334 det->momentsR1->data.F32[numGood] = NAN;335 det->momentsRH->data.F32[numGood] = NAN;336 det->kronFlux->data.F32[numGood] = NAN;337 det->kronFluxErr->data.F32[numGood] = NAN;338 det->kronFluxInner->data.F32[numGood] = NAN;339 det->kronFluxOuter->data.F32[numGood] = NAN;340 det->diffRP->data.F32[numGood] = NAN;341 det->diffSnP->data.F32[numGood] = NAN;342 det->diffRM->data.F32[numGood] = NAN;343 det->diffSnM->data.F32[numGood] = NAN;344 det->flags2->data.U32[numGood] = 0;345 det->ippIdet->data.U32[numGood] = 0;346 det->nFrames->data.U16[numGood] = 0;347 det->padding->data.S16[numGood] = 0;348 }349 }350 351 // Calculate error in RA, Dec352 double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); //SC: Warning! Promotion of F32353 double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); //SC: Warning! Promotion of F32354 double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); //SC: Warning! Promotion of F32355 double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); //SC: Warning! Promotion of F32356 357 if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||358 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||359 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||360 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {361 psTrace("ppMops.read", 10,362 "Discarding row %ld from input %d because of non-finite values: "363 "%f %f %lf %lf %f %f %f %f %f %f",364 j, i,365 det->x->data.F32[numGood], det->y->data.F32[numGood],366 det->ra->data.F64[numGood], det->dec->data.F64[numGood],367 det->mag->data.F32[numGood], det->magErr->data.F32[numGood],368 xErr, yErr, scale, angle);369 continue;370 }371 372 // XXX Not at all sure I've got the angles around the right way here...373 double cosAngle = cos(angle), sinAngle = sin(angle);374 double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);375 double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);376 double errScale = scale / 3600.0;377 det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);378 det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);379 380 det->mask->data.U8[numGood] = 0;381 plateScale += scale;382 numGood++;383 }384 #endif385 193 det->seeing *= ((float) plateScale) / ((float) numGood); 386 194 387 det->x->n = numGood; 388 det->y->n = numGood; 389 det->ra->n = numGood; 390 det->dec->n = numGood; 391 det->raErr->n = numGood; 392 det->decErr->n = numGood; 393 det->mag->n = numGood; 394 det->magErr->n = numGood; 395 det->chi2->n = numGood; 396 det->dof->n = numGood; 397 det->cr->n = numGood; 398 det->extended->n = numGood; 399 det->psfMajor->n = numGood; 400 det->psfMinor->n = numGood; 401 det->psfTheta->n = numGood; 402 det->quality->n = numGood; 403 det->numPix->n = numGood; 404 det->xxMoment->n = numGood; 405 det->xyMoment->n = numGood; 406 det->yyMoment->n = numGood; 407 det->flags->n = numGood; 408 det->diffSkyfileId->n = numGood; 409 det->naxis1->n = numGood; 410 det->naxis2->n = numGood; 411 det->mask->n = numGood; 412 det->nPos->n = numGood; 413 det->fPos->n = numGood; 414 det->nRatioBad->n = numGood; 415 det->nRatioMask->n = numGood; 416 det->nRatioAll->n = numGood; 417 det->psfInstFlux->n = numGood; 418 det->psfInstFluxSig->n = numGood; 419 det->apMag->n = numGood; 420 det->apMagRaw->n = numGood; 421 det->apMagRadius->n = numGood; 422 det->apFlux->n = numGood; 423 det->apFluxSig->n = numGood; 424 det->peakFluxAsMag->n = numGood; 425 det->calPsfMag->n = numGood; 426 det->calPsfMagSig->n = numGood; 427 det->sky->n = numGood; 428 det->skySig->n = numGood; 429 det->qualityPerfect->n = numGood; 430 det->momentsR1->n = numGood; 431 det->momentsRH->n = numGood; 432 det->kronFlux->n = numGood; 433 det->kronFluxErr->n = numGood; 434 det->kronFluxInner->n = numGood; 435 det->kronFluxOuter->n = numGood; 436 det->diffRP->n = numGood; 437 det->diffSnP->n = numGood; 438 det->diffRM->n = numGood; 439 det->diffSnM->n = numGood; 440 det->flags2->n = numGood; 441 det->ippIdet->n = numGood; 442 det->nFrames->n = numGood; 443 det->padding->n = numGood; 444 445 det->num = numGood; 195 // Are we using numGood for anything outside of this function? 196 det->numGood = numGood; 446 197 447 198 if (isfinite(args->zp) && numGood > 0) { 448 psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); 449 } 450 451 psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); 452 453 psFree(table); 454 detections->data[i] = det; 199 psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); 200 } 201 202 psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name); 455 203 } 456 204
Note:
See TracChangeset
for help on using the changeset viewer.
