IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 5, 2007, 12:47:04 PM (19 years ago)
Author:
Paul Price
Message:

Need to pass down a corrected view so that the particular chip can be identified.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmPSF_IO.c

    r15053 r15228  
    66 *  @author EAM, IfA
    77 *
    8  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-09-28 00:38:42 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-10-05 22:47:04 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    126126bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
    127127{
     128    pmFPAview *thisView = pmFPAviewAlloc (view->nRows);
     129    *thisView = *view;
    128130
    129131    for (int i = 0; i < fpa->chips->n; i++) {
    130132        pmChip *chip = fpa->chips->data[i];
    131         if (!pmPSFmodelWriteChip (chip, view, file, config)) {
     133        thisView->chip = i;
     134        if (!pmPSFmodelWriteChip (chip, thisView, file, config)) {
    132135            psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
    133             return false;
    134         }
    135     }
     136            psFree(thisView);
     137            return false;
     138        }
     139    }
     140    psFree(thisView);
    136141    return true;
    137142}
    … …  
    141146{
    142147    if (!pmPSFmodelWrite (chip->analysis, view, file, config)) {
    143         psError(PS_ERR_IO, false, "Failed to write PSF for chip");
    144         return false;
     148        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
     149        return false;
    145150    }
    146151    return true;
    … …  
    151156//   - a PHU blank header
    152157// - image header        : FITS Image NAXIS = 0
    153 // if (trendMode == MAP) 
     158// if (trendMode == MAP)
    154159//   - psf resid (+header) : FITS Image
    155160// else
    … …  
    186191    // define the EXTNAME values used for image header, table data, and residual image segments
    187192    {
    188         // lookup the EXTNAME values used for table data and image header segments
    189         char *rule = NULL;
    190 
    191         // Menu of EXTNAME rules
    192         psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
    193         if (!menu) {
    194             psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
    195             return false;
    196         }
    197 
    198         // EXTNAME for image header
    199         rule = psMetadataLookupStr(&status, menu, "PSF.HEAD");
    200         if (!rule) {
    201             psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
    202             return false;
    203         }
    204         headName = pmFPAfileNameFromRule (rule, file, view);
    205 
    206         // EXTNAME for table data
    207         rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
    208         if (!rule) {
    209             psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
    210             psFree (headName);
    211             return false;
    212         }
    213         tableName = pmFPAfileNameFromRule (rule, file, view);
    214 
    215         // EXTNAME for resid data
    216         rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
    217         if (!rule) {
    218             psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
    219             psFree (headName);
    220             psFree (tableName);
    221             return false;
    222         }
    223         residName = pmFPAfileNameFromRule (rule, file, view);
     193        // lookup the EXTNAME values used for table data and image header segments
     194        char *rule = NULL;
     195
     196        // Menu of EXTNAME rules
     197        psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
     198        if (!menu) {
     199            psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
     200            return false;
     201        }
     202
     203        // EXTNAME for image header
     204        rule = psMetadataLookupStr(&status, menu, "PSF.HEAD");
     205        if (!rule) {
     206            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
     207            return false;
     208        }
     209        headName = pmFPAfileNameFromRule (rule, file, view);
     210
     211        // EXTNAME for table data
     212        rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
     213        if (!rule) {
     214            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
     215            psFree (headName);
     216            return false;
     217        }
     218        tableName = pmFPAfileNameFromRule (rule, file, view);
     219
     220        // EXTNAME for resid data
     221        rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
     222        if (!rule) {
     223            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
     224            psFree (headName);
     225            psFree (tableName);
     226            return false;
     227        }
     228        residName = pmFPAfileNameFromRule (rule, file, view);
    224229    }
    225230
    … …  
    227232    // if this header block is new, write it to disk
    228233    if (hdu->header != file->header) {
    229         // add EXTNAME, EXTHEAD, EXTTYPE to header
    230         psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName);
    231         psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName);
    232         psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE");
    233         if (!file->wrote_phu) {
    234             // this hdu->header acts as the PHU: set EXTEND to be true
    235             psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true);
    236             file->wrote_phu = true;
    237         }
    238 
    239         psFitsWriteBlank (file->fits, hdu->header, headName);
    240         psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type);
    241         file->header = hdu->header;
    242         psFree (headName);
     234        // add EXTNAME, EXTHEAD, EXTTYPE to header
     235        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName);
     236        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName);
     237        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE");
     238        if (!file->wrote_phu) {
     239            // this hdu->header acts as the PHU: set EXTEND to be true
     240            psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true);
     241            file->wrote_phu = true;
     242        }
     243
     244        psFitsWriteBlank (file->fits, hdu->header, headName);
     245        psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type);
     246        file->header = hdu->header;
     247        psFree (headName);
    243248    }
    244249
    … …  
    246251    pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
    247252    if (!psf) {
    248         psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
    249         psFree (tableName);
    250         psFree (residName);
    251         return false;
     253        psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
     254        psFree (tableName);
     255        psFree (residName);
     256        return false;
    252257    }
    253258
    254259    // write the PSF model parameters in a FITS table
    255     { 
    256         // we need to write a header for the table,
    257         psMetadata *header = psMetadataAlloc();
    258 
    259         char *modelName = pmModelClassGetName (psf->type);
    260         psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
    261 
    262         psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
    263         psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
    264         psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
    265 
    266         int nPar = pmModelClassParameterCount (psf->type)    ;
    267         psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
    268 
    269         psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
    270         psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
    271         psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
    272         psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
    273 
    274         // extract PSF Clump info
    275         pmPSFClump psfClump;
     260    {
     261        // we need to write a header for the table,
     262        psMetadata *header = psMetadataAlloc();
     263
     264        char *modelName = pmModelClassGetName (psf->type);
     265        psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
     266
     267        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
     268        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
     269        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
     270
     271        int nPar = pmModelClassParameterCount (psf->type)    ;
     272        psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
     273
     274        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
     275        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
     276        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
     277        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
     278
     279        // extract PSF Clump info
     280        pmPSFClump psfClump;
    276281
    277282        psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");   assert (status);
    … …  
    280285        psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  assert (status);
    281286
    282         psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
    283         psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
    284         psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
    285         psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
    286 
    287         // save the dimensions of each parameter
    288         for (int i = 0; i < nPar; i++) {
    289             char name[9];
    290             int nX, nY;
    291 
    292             pmTrend2D *trend = psf->params->data[i];
    293             if (trend == NULL) continue;
    294 
    295             if (trend->mode == PM_TREND_MAP) {
    296               nX = trend->map->map->numCols;
    297               nY = trend->map->map->numRows;
    298             } else {
    299               nX = trend->poly->nX;
    300               nY = trend->poly->nY;
    301             }
    302             snprintf (name, 9, "PAR%02d_NX", i);
    303             psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
    304             snprintf (name, 9, "PAR%02d_NY", i);
    305             psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
    306             snprintf (name, 9, "PAR%02d_MD", i);
    307             char *modeName = pmTrend2DModeToString (trend->mode);
    308             psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
    309             psFree (modeName);
    310         }
    311 
    312         // other required information describing the PSF
    313         psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
    314         psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
    315         psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
    316         psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
    317 
    318         // XXX can we drop this now?
    319         psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
    320 
    321         // build a FITS table of the PSF parameters
    322         psArray *psfTable = psArrayAllocEmpty (100);
    323         for (int i = 0; i < nPar; i++) {
    324             pmTrend2D *trend = psf->params->data[i];
    325             if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
    326 
    327             if (trend->mode == PM_TREND_MAP) {
    328                 // write the image components into a table: this is needed because they may each be a different size
    329                 psImageMap *map = trend->map;
    330                 for (int ix = 0; ix < map->map->numCols; ix++) {
    331                     for (int iy = 0; iy < map->map->numRows; iy++) {
    332                         psMetadata *row = psMetadataAlloc ();
    333                         psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
    334                         psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
    335                         psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
    336                         psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
    337                         psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
    338                         psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
    339 
    340                         psArrayAdd (psfTable, 100, row);
    341                         psFree (row);
    342                     }
    343                 }
    344             } else {
    345                 // write the polynomial components into a table
    346                 psPolynomial2D *poly = trend->poly;
    347                 for (int ix = 0; ix <= poly->nX; ix++) {
    348                     for (int iy = 0; iy <= poly->nY; iy++) {
    349                         psMetadata *row = psMetadataAlloc ();
    350                         psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
    351                         psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
    352                         psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
    353                         psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
    354                         psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
    355                         psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
    356 
    357                         psArrayAdd (psfTable, 100, row);
    358                         psFree (row);
    359                     }
    360                 }
    361             }
    362         }
    363 
    364         // write an empty FITS segment if we have no PSF information
    365         if (psfTable->n == 0) {
    366             // XXX this is probably an error (if we have a PSF, how do we have no data?)
    367             psFitsWriteBlank (file->fits, header, tableName);
    368         } else {
    369             psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
    370             if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
    371                 psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
    372                 psFree (tableName);
    373                 psFree (residName);
    374                 psFree (psfTable);
    375                 psFree (header);
    376                 return false;
    377             }
    378         }
    379         psFree (tableName);
    380         psFree (psfTable);
    381         psFree (header);
     287        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
     288        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
     289        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
     290        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
     291
     292        // save the dimensions of each parameter
     293        for (int i = 0; i < nPar; i++) {
     294            char name[9];
     295            int nX, nY;
     296
     297            pmTrend2D *trend = psf->params->data[i];
     298            if (trend == NULL) continue;
     299
     300            if (trend->mode == PM_TREND_MAP) {
     301              nX = trend->map->map->numCols;
     302              nY = trend->map->map->numRows;
     303            } else {
     304              nX = trend->poly->nX;
     305              nY = trend->poly->nY;
     306            }
     307            snprintf (name, 9, "PAR%02d_NX", i);
     308            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
     309            snprintf (name, 9, "PAR%02d_NY", i);
     310            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
     311            snprintf (name, 9, "PAR%02d_MD", i);
     312            char *modeName = pmTrend2DModeToString (trend->mode);
     313            psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
     314            psFree (modeName);
     315        }
     316
     317        // other required information describing the PSF
     318        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
     319        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
     320        psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
     321        psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
     322
     323        // XXX can we drop this now?
     324        psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
     325
     326        // build a FITS table of the PSF parameters
     327        psArray *psfTable = psArrayAllocEmpty (100);
     328        for (int i = 0; i < nPar; i++) {
     329            pmTrend2D *trend = psf->params->data[i];
     330            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
     331
     332            if (trend->mode == PM_TREND_MAP) {
     333                // write the image components into a table: this is needed because they may each be a different size
     334                psImageMap *map = trend->map;
     335                for (int ix = 0; ix < map->map->numCols; ix++) {
     336                    for (int iy = 0; iy < map->map->numRows; iy++) {
     337                        psMetadata *row = psMetadataAlloc ();
     338                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     339                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     340                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     341                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
     342                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
     343                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
     344
     345                        psArrayAdd (psfTable, 100, row);
     346                        psFree (row);
     347                    }
     348                }
     349            } else {
     350                // write the polynomial components into a table
     351                psPolynomial2D *poly = trend->poly;
     352                for (int ix = 0; ix <= poly->nX; ix++) {
     353                    for (int iy = 0; iy <= poly->nY; iy++) {
     354                        psMetadata *row = psMetadataAlloc ();
     355                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
     356                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
     357                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
     358                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
     359                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
     360                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
     361
     362                        psArrayAdd (psfTable, 100, row);
     363                        psFree (row);
     364                    }
     365                }
     366            }
     367        }
     368
     369        // write an empty FITS segment if we have no PSF information
     370        if (psfTable->n == 0) {
     371            // XXX this is probably an error (if we have a PSF, how do we have no data?)
     372            psFitsWriteBlank (file->fits, header, tableName);
     373        } else {
     374            psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
     375            if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
     376                psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
     377                psFree (tableName);
     378                psFree (residName);
     379                psFree (psfTable);
     380                psFree (header);
     381                return false;
     382            }
     383        }
     384        psFree (tableName);
     385        psFree (psfTable);
     386        psFree (header);
    382387    }
    383388
    384389    // write the residual images (3D)
    385390    {
    386         psMetadata *header = psMetadataAlloc ();
    387         if (psf->residuals == NULL) {
    388             // set some header keywords to make it clear there are no residuals?
    389             psFitsWriteBlank (file->fits, header, residName);
    390             psFree (residName);
    391             psFree (header);
    392             return true;
    393         }
    394 
    395         psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
    396         psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
    397         psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
    398         psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
    399 
    400         // write the residuals as three planes of the image
    401         // this call creates an extension with NAXIS3 = 3
    402         if (psf->residuals->Rx) {
    403             // this call creates an extension with NAXIS3 = 3
    404             psArray *images = psArrayAllocEmpty (3);
    405             psArrayAdd (images, 1, psf->residuals->Ro);
    406             psArrayAdd (images, 1, psf->residuals->Rx);
    407             psArrayAdd (images, 1, psf->residuals->Ry);
    408 
    409             psFitsWriteImageCube (file->fits, header, images, residName);
    410             psFree (images);
    411         } else {
    412             // this call creates an extension with NAXIS3 = 1
    413             psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName);
    414         }
    415         psFree (residName);
    416         psFree (header);
     391        psMetadata *header = psMetadataAlloc ();
     392        if (psf->residuals == NULL) {
     393            // set some header keywords to make it clear there are no residuals?
     394            psFitsWriteBlank (file->fits, header, residName);
     395            psFree (residName);
     396            psFree (header);
     397            return true;
     398        }
     399
     400        psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
     401        psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
     402        psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
     403        psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
     404
     405        // write the residuals as three planes of the image
     406        // this call creates an extension with NAXIS3 = 3
     407        if (psf->residuals->Rx) {
     408            // this call creates an extension with NAXIS3 = 3
     409            psArray *images = psArrayAllocEmpty (3);
     410            psArrayAdd (images, 1, psf->residuals->Ro);
     411            psArrayAdd (images, 1, psf->residuals->Rx);
     412            psArrayAdd (images, 1, psf->residuals->Ry);
     413
     414            psFitsWriteImageCube (file->fits, header, images, residName);
     415            psFree (images);
     416        } else {
     417            // this call creates an extension with NAXIS3 = 1
     418            psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName);
     419        }
     420        psFree (residName);
     421        psFree (header);
    417422    }
    418423
    … …  
    465470    psMetadata *outhead = psMetadataAlloc();
    466471    if (phu) {
    467         psMetadataCopy (outhead, phu->header);
     472        psMetadataCopy (outhead, phu->header);
    468473    } else {
    469         pmConfigConformHeader (outhead, file->format);
    470 
    471         psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information
    472         const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME");
    473         if (fpaNameHdr && strlen(fpaNameHdr) > 0) {
    474             const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME");
    475             psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName);
    476         }
     474        pmConfigConformHeader (outhead, file->format);
     475
     476        psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information
     477        const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME");
     478        if (fpaNameHdr && strlen(fpaNameHdr) > 0) {
     479            const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME");
     480            psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName);
     481        }
    477482    }
    478483
    … …  
    524529{
    525530    if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
    526         psError(PS_ERR_IO, false, "Failed to write PSF for chip");
    527         return false;
     531        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
     532        return false;
    528533    }
    529534    return true;
    … …  
    628633        snprintf (name, 9, "PAR%02d_NX", i);
    629634        binning->nXruff = psMetadataLookupS32 (&status, header, name);
    630         if (!status) continue;          // not all parameters are defined
     635        if (!status) continue;          // not all parameters are defined
    631636
    632637        snprintf (name, 9, "PAR%02d_NY", i);
    633638        binning->nYruff = psMetadataLookupS32 (&status, header, name);
    634         if (!status) {
    635             psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
    636             return false;
    637         }
     639        if (!status) {
     640            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
     641            return false;
     642        }
    638643
    639644        snprintf (name, 9, "PAR%02d_MD", i);
    640645        char *modeName = psMetadataLookupStr (&status, header, name);
    641         if (!status) {
    642             psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
    643             return false;
    644         }
    645         pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
    646         if (psfTrendMode == PM_TREND_NONE) {
    647             psfTrendMode = PM_TREND_POLY_ORD;
    648         }
    649 
    650         psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
    651         psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
    652         psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
     646        if (!status) {
     647            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
     648            return false;
     649        }
     650        pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
     651        if (psfTrendMode == PM_TREND_NONE) {
     652            psfTrendMode = PM_TREND_POLY_ORD;
     653        }
     654
     655        psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
     656        psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
     657        psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
    653658    }
    654659    psFree (binning);
    … …  
    674679
    675680        pmTrend2D *trend = psf->params->data[iPar];
    676         if (trend == NULL) {
    677             psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
    678             return false;
    679         }
    680 
    681         if (trend->mode == PM_TREND_MAP) {
    682             psImageMap *map = trend->map;
    683             assert (map);
    684             assert (map->map);
    685             assert (map->error);
    686             assert (xPow >= 0);
    687             assert (yPow >= 0);
    688             assert (xPow < map->map->numCols);
    689             assert (yPow < map->map->numRows);
    690             map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    691             map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
    692         } else {
    693             psPolynomial2D *poly = trend->poly;
    694             assert (poly);
    695             assert (xPow >= 0);
    696             assert (yPow >= 0);
    697             assert (xPow <= poly->nX);
    698             assert (yPow <= poly->nY);
    699             poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
    700             poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
    701             poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
    702         }
     681        if (trend == NULL) {
     682            psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
     683            return false;
     684        }
     685
     686        if (trend->mode == PM_TREND_MAP) {
     687            psImageMap *map = trend->map;
     688            assert (map);
     689            assert (map->map);
     690            assert (map->error);
     691            assert (xPow >= 0);
     692            assert (yPow >= 0);
     693            assert (xPow < map->map->numCols);
     694            assert (yPow < map->map->numRows);
     695            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     696            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
     697        } else {
     698            psPolynomial2D *poly = trend->poly;
     699            assert (poly);
     700            assert (xPow >= 0);
     701            assert (yPow >= 0);
     702            assert (xPow <= poly->nX);
     703            assert (yPow <= poly->nY);
     704            poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
     705            poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
     706            poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
     707        }
    703708    }
    704709    psFree (header);
Note: See TracChangeset for help on using the changeset viewer.