Changeset 28484 for branches/pap/ppTranslate/src/ppMopsRead.c
- Timestamp:
- Jun 24, 2010, 2:59:09 PM (16 years ago)
- Location:
- branches/pap
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
ppTranslate/src/ppMopsRead.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/pap
- Property svn:mergeinfo changed
-
branches/pap/ppTranslate/src/ppMopsRead.c
r28027 r28484 4 4 5 5 #include <stdio.h> 6 #include <string.h> 6 7 #include <pslib.h> 7 8 … … 16 17 psArray *detections = psArrayAlloc(num); // Array of detections, to return 17 18 for (int i = 0; i < num; i++) { 18 psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file 19 const char *name = inNames->data[i]; // File name 20 psFits *fits = psFitsOpen(name, "r"); // FITS file 19 21 if (!fits) { 20 22 psError(PS_ERR_IO, false, "Unable to open input %d", i); … … 24 26 if (!header) { 25 27 psError(PS_ERR_IO, false, "Unable to read header %d", i); 26 return false;27 }28 29 psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image30 if (diffSkyfileId == 0) {31 psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);32 28 return false; 33 29 } … … 47 43 continue; 48 44 } 49 ppMopsDetections *det = ppMopsDetectionsAlloc(size); 50 51 psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]); 45 ppMopsDetections *det = detections->data[i] = ppMopsDetectionsAlloc(); 46 det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension 47 det->header = header; 48 det->num = size; 52 49 53 50 det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA")); … … 64 61 psMetadataLookupF32(NULL, header, "FWHM_MIN")); 65 62 66 int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns67 int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows63 det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); 64 det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); 68 65 69 psFree(header); 70 71 psArray *table = psFitsReadTable(fits); // Table of interest 72 if (!table) { 73 psError(PS_ERR_IO, false, "Unable to read table %d", i); 66 det->psfHeader = psFitsReadHeader(NULL, fits); 67 if (!det->psfHeader) { 68 psError(psErrorCodeLast(), false, "Unable to read header %d", i); 74 69 return false; 75 70 } 76 psFitsClose(fits); 77 78 double plateScale = 0.0; // Plate scale 79 long numGood = 0; // Number of good rows 80 for (long j = 0; j < size; j++) { 81 psMetadata *row = table->data[j]; // Row of interest 82 83 psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS"); 84 if (flags & SOURCE_MASK) { 85 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags); 86 continue; 87 } 88 89 det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF"); 90 det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF"); 91 det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF")); 92 det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF")); 93 det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG"); 94 det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG"); 95 det->chi2->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_CHISQ"); 96 det->dof->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NDOF"); 97 det->cr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "CR_NSIGMA"); 98 det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA"); 99 det->psfMajor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MAJOR"); 100 det->psfMinor->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_MINOR"); 101 det->psfTheta->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_THETA"); 102 det->quality->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_QF"); 103 det->numPix->data.S32[numGood] = psMetadataLookupS32(NULL, row, "PSF_NPIX"); 104 det->xxMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XX"); 105 det->xyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_XY"); 106 det->yyMoment->data.F32[numGood] = psMetadataLookupF32(NULL, row, "MOMENTS_YY"); 107 det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS"); 108 det->diffSkyfileId->data.S64[numGood] = diffSkyfileId; 109 det->naxis1->data.S32[numGood] = naxis1; 110 det->naxis2->data.S32[numGood] = naxis2; 111 112 det->nPos->data.S32[numGood] = psMetadataLookupS32(NULL, row, "DIFF_NPOS"); 113 det->fPos->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_FRATIO"); 114 det->nRatioBad->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_BAD"); 115 det->nRatioMask->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_MASK"); 116 det->nRatioAll->data.F32[numGood] = psMetadataLookupF32(NULL, row, "DIFF_NRATIO_ALL"); 117 118 // Calculate error in RA, Dec 119 double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG"); 120 double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG"); 121 double scale = psMetadataLookupF64(NULL, row, "PLTSCALE"); 122 double angle = psMetadataLookupF64(NULL, row, "POSANGLE"); 123 124 if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) || 125 !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) || 126 !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) || 127 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) { 128 psTrace("ppMops.read", 10, 129 "Discarding row %ld from input %d because of non-finite values: " 130 "%f %f %lf %lf %f %f %f %f %f %f", 131 j, i, 132 det->x->data.F32[numGood], det->y->data.F32[numGood], 133 det->ra->data.F64[numGood], det->dec->data.F64[numGood], 134 det->mag->data.F32[numGood], det->magErr->data.F32[numGood], 135 xErr, yErr, scale, angle); 136 continue; 137 } 138 139 // XXX Not at all sure I've got the angles around the right way here... 140 double cosAngle = cos(angle), sinAngle = sin(angle); 141 double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle); 142 double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr); 143 double errScale = scale / 3600.0; 144 det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2); 145 det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2); 146 147 det->mask->data.U8[numGood] = 0; 148 plateScale += scale; 149 numGood++; 150 } 151 det->seeing *= plateScale / numGood; 152 153 det->x->n = numGood; 154 det->y->n = numGood; 155 det->ra->n = numGood; 156 det->dec->n = numGood; 157 det->raErr->n = numGood; 158 det->decErr->n = numGood; 159 det->mag->n = numGood; 160 det->magErr->n = numGood; 161 det->chi2->n = numGood; 162 det->dof->n = numGood; 163 det->cr->n = numGood; 164 det->extended->n = numGood; 165 det->psfMajor->n = numGood; 166 det->psfMinor->n = numGood; 167 det->psfTheta->n = numGood; 168 det->quality->n = numGood; 169 det->numPix->n = numGood; 170 det->xxMoment->n = numGood; 171 det->xyMoment->n = numGood; 172 det->yyMoment->n = numGood; 173 det->flags->n = numGood; 174 det->diffSkyfileId->n = numGood; 175 det->naxis1->n = numGood; 176 det->naxis2->n = numGood; 177 det->mask->n = numGood; 178 det->nPos->n = numGood; 179 det->fPos->n = numGood; 180 det->nRatioBad->n = numGood; 181 det->nRatioMask->n = numGood; 182 det->nRatioAll->n = numGood; 183 184 det->num = numGood; 185 186 if (isfinite(args->zp) && numGood > 0) { 187 psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32)); 71 psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest 72 if (!table) { 73 psError(psErrorCodeLast(), false, "Unable to read table %d", i); 74 return false; 188 75 } 189 76 190 psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]); 77 psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF"); 78 psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF"); 191 79 192 psFree(table); 193 detections->data[i] = det; 80 det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); 81 det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64)); 82 det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF")); 83 det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF")); 84 if (!det->ra || !det->dec || !det->x || !det->y) { 85 psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns"); 86 return false; 87 } 88 89 psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component); 90 91 if (!psFitsMoveExtName(fits, "SkyChip.deteff")) { 92 psWarning("No detection efficiencies included in %s", det->component); 93 psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component); 94 psErrorClear(); 95 continue; 96 } 97 98 det->deteffHeader = psFitsReadHeader(NULL, fits); 99 if (!det->deteffHeader) { 100 psWarning("Unable to read detection efficiency header in %s", det->component); 101 psErrorClear(); 102 continue; 103 } 104 det->deteffTable = psFitsReadTableAllColumns(fits); 105 if (!det->deteffTable) { 106 psFree(det->deteffHeader); 107 det->deteffHeader = NULL; 108 psWarning("Unable to read detection efficiency table in %s", det->component); 109 psErrorClear(); 110 continue; 111 } 112 psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component); 113 psFitsClose(fits); 194 114 } 195 115
Note:
See TracChangeset
for help on using the changeset viewer.
