- Timestamp:
- Sep 21, 2007, 4:20:17 PM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branch_20070921/psModules/src/objects/pmPSF_IO.c
r14969 r14980 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.22.2. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-09-2 1 18:55:12$8 * @version $Revision: 1.22.2.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-22 02:20:11 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 148 148 149 149 // for a pmPSF supplied on the analysis metadata, we write out 150 // if needed: 151 // - a PHU blank header 150 152 // - image header : FITS Image NAXIS = 0 151 // - psf table (+header) : FITS Table152 // XXX if we are using psImageMap, this is also a FITS Image, with header data153 // - psf resid (+header) : FITS Image154 // if needed, we also write out a PHU blank header153 // if (trendMode == MAP) 154 // - psf resid (+header) : FITS Image 155 // else 156 // - psf table (+header) : FITS Table 155 157 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) 156 158 { … … 216 218 217 219 // write out the IMAGE header segment 218 // this header block is new, write it to disk220 // if this header block is new, write it to disk 219 221 if (hdu->header != file->header) { 220 222 // add EXTNAME, EXTHEAD, EXTTYPE to header … … 244 246 245 247 // write the PSF model parameters in a FITS table 248 assert (psf->psfTrendMode != PM_TREND_NONE); 246 249 { 247 250 // we need to write a header for the table, … … 267 270 } 268 271 272 char *modeName = pmTrend2DModeToString(trend->mode); 273 269 274 // other required information describing the PSF 270 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); 271 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); 272 psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); 273 psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); 275 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); 276 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); 277 psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); 278 psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); 279 psMetadataAddStr (header, PS_LIST_TAIL, "MODE", 0, "", modeName); 280 psFree (modeName); 274 281 275 282 // XXX can we drop this now? … … 279 286 psArray *psfTable = psArrayAllocEmpty (100); 280 287 for (int i = 0; i < nPar; i++) { 281 psPolynomial2D *poly = psf->params->data[i]; 282 if (poly == NULL) continue; // skip unset parameters (eg, XPOS) 283 for (int ix = 0; ix <= poly->nX; ix++) { 284 for (int iy = 0; iy <= poly->nY; iy++) { 285 286 psMetadata *row = psMetadataAlloc (); 287 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 288 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 289 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 290 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 291 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 292 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 293 294 psArrayAdd (psfTable, 100, row); 295 psFree (row); 288 pmTrend2D *trend = psf->params->data[i]; 289 if (trend == NULL) continue; // skip unset parameters (eg, XPOS) 290 291 if (trend->mode == PM_TREND_MAP) { 292 // write the image components into a table: this is needed because they may each be a different size 293 psImageMap *map = trend->map; 294 for (int ix = 0; ix <= map->map->numCols; ix++) { 295 for (int iy = 0; iy <= map->map->numRows; iy++) { 296 psMetadata *row = psMetadataAlloc (); 297 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 298 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 299 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 300 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", map->map->data.F32[iy][ix]); 301 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", map->error->data.F32[ix][iy]); 302 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", 0); // no cells are masked 303 304 psArrayAdd (psfTable, 100, row); 305 psFree (row); 306 } 307 } 308 } else { 309 // write the polynomial components into a table 310 psPolynomial2D *poly = trend->poly; 311 for (int ix = 0; ix <= poly->nX; ix++) { 312 for (int iy = 0; iy <= poly->nY; iy++) { 313 psMetadata *row = psMetadataAlloc (); 314 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 315 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 316 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 317 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 318 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 319 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 320 321 psArrayAdd (psfTable, 100, row); 322 psFree (row); 323 } 296 324 } 297 325 } … … 535 563 int nYorder = psMetadataLookupS32 (&status, header, name); 536 564 if (!status) { 537 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i);565 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i); 538 566 return false; 539 567 } 568 snprintf (name, 9, "PAR%02d_MD", i); 569 char *modeName = psMetadataLookupStr (&status, header, name); 570 if (!status) { 571 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i); 572 return false; 573 } 574 // XXX allocate pmTrend2d based on the values here 540 575 psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); 541 576 } … … 549 584 // XXX can we drop this now? 550 585 psf->skyBias = psMetadataLookupF32 (&status, header, "SKY_BIAS"); 586 587 char *modeName = psMetadataLookupStr (&status, header, "MODE"); 588 trend->mode = pmTrend2DModeToString(modeName); 551 589 552 590 // read the raw table data … … 561 599 // XXX sanity check here 562 600 563 p sPolynomial2D *poly= psf->params->data[iPar];564 if ( poly== NULL) {565 psError(PS_ERR_UNKNOWN, true, " values for parameter %d, but missing NX", iPar);601 pmTrend2D *trend = psf->params->data[iPar]; 602 if (trend == NULL) { 603 psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar); 566 604 return false; 567 605 } 568 606 569 poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); 570 poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); 571 poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); 607 if (trend->mode == PM_TREND_MAP) { 608 psImageMap *map = trend->map; 609 map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); 610 map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); 611 // poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); 612 } else { 613 psPolynomial2D *poly = trend->poly; 614 poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); 615 poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); 616 poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); 617 } 572 618 } 573 619 psFree (header);
Note:
See TracChangeset
for help on using the changeset viewer.
