IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 15, 2006, 8:30:16 AM (21 years ago)
Author:
eugene
Message:

API cleanup, removed old test files, some re-organization, changed FLT/GAL to EXT

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotOutput.c

    r5986 r5993  
    1414
    1515    if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
    16     if (psfFile != NULL) psphotWritePSF (psf, psfFile);
    1716    psphotSamplePSFs (config, psf, imdata->image);
    1817
     18    if (psfFile != NULL) {
     19        psMetadata *psfData = pmPSFtoMD (NULL, psf);
     20        psMetadataConfigWrite (psfData, psfFile);
     21        psFree (psfData);
     22    }
    1923    if (outputFile == NULL) {
    2024        psLogMsg ("output", 3, "no data output file selected");
    … …  
    166170    float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
    167171    if (!status) ZERO_POINT = 25.0;
    168     char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE");
    169172
    170173    // write necessary information to output header
    171     psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS",     sources->n);
    172     psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT",  PS_DATA_F32 | PS_META_REPLACE, "zero point",          ZERO_POINT);
    173     psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    174     psMetadataAdd (imdata->header, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    175     psMetadataAdd (imdata->header, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
    176     psMetadataAdd (imdata->header, PS_LIST_TAIL, "NPSFSTAR", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nPSFstars);
    177     psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
    178     psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
    179     psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code",     PHOTCODE);
    180     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
    181     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
    182     psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
    183     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
    184     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
     174    psphotUpdateHeader (imdata->header, config);
     175    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n);
    185176
    186177    // create file, write-out header
    … …  
    227218        psLineAdd (line, "%2d ",   type);
    228219        psLineAdd (line, "%3.1f ", lsky);
    229         psLineAdd (line, "%6.3f ", 99.999);    // should be 'Mgal
     220        psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal
    230221        psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT));
    231         psLineAdd (line, "%6.2f ", PAR[4]);  // should be 'FHWM x'
    232         psLineAdd (line, "%6.2f ", PAR[5]);  // should be 'FHWM y'
    233         psLineAdd (line, "%5.1f\n", 0);   // should be theta
     222        psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x'
     223        psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y'
     224        psLineAdd (line, "%5.1f\n", 0); // should be theta
    234225        fwrite (line->line, 1, line->Nline, f);
    235226    }
    … …  
    256247    float RADIUS = psMetadataLookupF32 (&status, config, "AP_RADIUS");
    257248    float ZERO_POINT = psMetadataLookupF32 (&status, config, "ZERO_POINT");
    258     char *PHOTCODE = psMetadataLookupPtr (&status, config, "PHOTCODE");
    259249
    260250    // write necessary information to output header
    261     psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS",   PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS",     sources->n);
    262     psMetadataAdd (imdata->header, PS_LIST_TAIL, "ZERO_PT",  PS_DATA_F32 | PS_META_REPLACE, "zero point",          ZERO_POINT);
    263     psMetadataAdd (imdata->header, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    264     psMetadataAdd (imdata->header, PS_LIST_TAIL, "dAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    265     psMetadataAdd (imdata->header, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
    266     psMetadataAdd (imdata->header, PS_LIST_TAIL, "PHOTCODE", PS_DATA_STRING | PS_META_REPLACE, "photometry code",     PHOTCODE);
    267     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
    268     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
    269     psMetadataAdd (imdata->header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
    270     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
    271     psMetadataAdd (imdata->header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
    272 
    273     // set NAXIS to 0
    274     mdi = psMetadataLookup (imdata->header, "NAXIS");
    275     mdi->data.S32 = 0;
    276     mdi->type = PS_DATA_S32;
     251    psphotUpdateHeader (imdata->header, config);
    277252
    278253    table = psArrayAlloc (sources->n);
    … …  
    299274        psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", type);
    300275        psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP",  PS_DATA_F32, "", lsky);
    301         psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0);    // should be 'Mgal'
     276        psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", 32.0);
    302277        psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X",  PS_DATA_F32, "", source->apMag + ZERO_POINT);
    303         psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y",  PS_DATA_F32, "", PAR[4]);  // should be 'FHWM x'
    304         psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", PAR[5]);  // should be 'FHWM y'
    305         psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_TYPE_S32, "", 0);   // should be theta
    306         // psMetadataAdd (header, PS_LIST_HEAD, "DUMMY",   PS_DATA_STRING, "", NULL);
     278        psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y",  PS_DATA_F32, "", PAR[4]);
     279        psMetadataAdd (row, PS_LIST_TAIL, "THETA",   PS_DATA_F32, "", PAR[5]);
     280        psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT",  PS_DATA_STRING, "", "0");
     281        psMetadataAdd (row, PS_LIST_TAIL, "DUMMY",   PS_DATA_STRING, "", "123");
    307282   
    308283        psArrayAdd (table, 100, row);
    … …  
    310285    }
    311286
    312     // psAbort ("psphotOutput", "FITS table output not finished");
    313 
    314     // psImage *image = psImageAlloc (1, 1, PS_TYPE_F32);
    315     // XXX missing: theader = psFitsHeaderFromTable (NULL, table, "SMPFILE");
    316     // psFitsWriteImage (fits, imdata->header, image, 0);
    317 
     287    // write the correct number of sources in the header
     288    psMetadataAdd (imdata->header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n);
     289    psMetadataAdd (imdata->header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true);
     290
     291    // set NAXIS to 0 (we don't write out the data array)
     292    mdi = psMetadataLookup (imdata->header, "NAXIS");
     293    mdi->data.S32 = 0;
     294    mdi->type = PS_DATA_S32;
     295
     296    // create the basic table header
    318297    theader = psMetadataAlloc ();
    319298    psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");
    320299
    321     psImage *image;
    322 
    323     image = psImageAlloc (512, 512, PS_DATA_F32);
    324     for (int j = 0; j < 512; j++) {
    325         for (int i = 0; i < 512; i++) {
    326             image->data.F32[j][i] = i;
    327         }
    328     }
    329 
     300    // write out the header and table
    330301    psFits *fits = psFitsOpen (filename, "w");
    331     // psFitsWriteHeader (imdata->header, fits);
    332     // psFitsWriteTable (fits, theader, table);
    333 
    334     psFitsWriteImage (fits, NULL, image, 0);
    335     psFitsWriteTable (fits, NULL, table);
     302    psFitsWriteHeader (imdata->header, fits);
     303    psFitsWriteTable (fits, theader, table);
    336304    psFitsClose (fits);
    337305
    … …  
    348316    pmModelWritePSFs (sources, config, name, psf);
    349317
    350     sprintf (name, "%s.flt.dat", filename);
    351     pmModelWriteFLTs (sources, name);
     318    sprintf (name, "%s.ext.dat", filename);
     319    pmModelWriteEXTs (sources, name);
    352320
    353321    sprintf (name, "%s.nul.dat", filename);
    … …  
    417385
    418386// dump the sources to an output file
    419 bool pmModelWriteFLTs (psArray *sources, char *filename) {
     387bool pmModelWriteEXTs (psArray *sources, char *filename) {
    420388
    421389    double dPos, dMag;
    … …  
    427395    f = fopen (filename, "w");
    428396    if (f == NULL) {
    429         psLogMsg ("pmModelWriteFLTs", 3, "can't open output file for moments%s\n", filename);
     397        psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename);
    430398        return false;
    431399    }
    … …  
    435403        pmSource *source = (pmSource *) sources->data[i];
    436404
    437         if (source->type != PM_SOURCE_GALAXY) continue;
     405        if (source->type != PM_SOURCE_EXTENDED) continue;
    438406        model = pmSourceMagnitudes (source, NULL, 0.0);
    439407        if (model == NULL) continue;
    … …  
    491459        source = sources->data[i];
    492460
    493         // skip these sources (in PSF or FLT)
     461        // skip these sources (in PSF or EXT)
    494462        if (source->type == PM_SOURCE_STAR) continue;
    495         if (source->type == PM_SOURCE_GALAXY) continue;
     463        if (source->type == PM_SOURCE_EXTENDED) continue;
    496464           
    497465        if (source->moments == NULL) {
    … …  
    580548        return (1);
    581549
    582       case PM_SOURCE_GALAXY:
     550      case PM_SOURCE_EXTENDED:
    583551        return (2);
    584552
    … …  
    597565}
    598566
    599 psImage *pmModelPSFatXY (psImage *image, pmModel *modelFLT, pmPSF *psf, int x, int y, int dx, int dy) {
     567psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) {
    600568
    601569    psRegion region = {x - dx, x + dx, y - dy, y + dy};
    … …  
    603571    psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32);
    604572    psImageInit (sample, 0);
    605     modelFLT->params->data.F32[2] = x;
    606     modelFLT->params->data.F32[3] = y;
    607     pmModel *modelPSF = pmModelFromPSF (modelFLT, psf);
     573    modelEXT->params->data.F32[2] = x;
     574    modelEXT->params->data.F32[3] = y;
     575    pmModel *modelPSF = pmModelFromPSF (modelEXT, psf);
    608576    pmSourceAddModel (sample, NULL, modelPSF, false, false);
    609577    psFree (modelPSF);
    … …  
    625593    if (output[0] == 0) return false;
    626594
    627     pmModel *modelFLT = pmModelAlloc (psf->type);
    628     modelFLT->params->data.F32[0] = 0;
    629     modelFLT->params->data.F32[1] = 1;
     595    pmModel *modelEXT = pmModelAlloc (psf->type);
     596    modelEXT->params->data.F32[0] = 0;
     597    modelEXT->params->data.F32[1] = 1;
    630598
    631599    psFits *fits = psFitsOpen (output, "w");
    632600
    633     sample = pmModelPSFatXY (image, modelFLT, psf, 25, 25, 25, 25);
     601    sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25);
    634602    psFitsWriteImage (fits, NULL, sample, 0);
    635     sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
     603    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25);
    636604    psFitsWriteImage (fits, NULL, sample, 0);
    637     sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols - 25, 25, 25, 25);
     605    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25);
    638606    psFitsWriteImage (fits, NULL, sample, 0);
    639     sample = pmModelPSFatXY (image, modelFLT, psf, 25, image->numRows - 25, 25, 25);
     607    sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25);
    640608    psFitsWriteImage (fits, NULL, sample, 0);
    641     sample = pmModelPSFatXY (image, modelFLT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
     609    sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25);
    642610    psFitsWriteImage (fits, NULL, sample, 0);
    643611
    644612    psFitsClose (fits);
    645613    return (TRUE);
    646 }
    647 
    648 psPolynomial2D *psPolynomial2DfromMD (psMetadata *folder) {
    649 
    650     bool status;
    651     char keyword[80];
    652 
    653     // get polynomial orders
    654     // XXX add status failures tests
    655     int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
    656     int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
    657 
    658     psPolynomial2D *poly = psPolynomial2DAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
    659 
    660     for (int nx = 0; nx < poly->nX + 1; nx++) {
    661         for (int ny = 0; ny < poly->nY + 1; ny++) {
    662             sprintf (keyword, "VAL_X%02d_Y%02d", nx, ny);
    663             poly->coeff[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
    664             if (!status) poly->mask[nx][ny] = 1;
    665 
    666             sprintf (keyword, "ERR_X%02d_Y%02d", nx, ny);
    667             poly->coeffErr[nx][ny] = psMetadataLookupF32 (&status, folder, keyword);
    668         }
    669     }
    670     return (poly);
    671 }
    672 
    673 psPolynomial3D *psPolynomial3DfromMD (psMetadata *folder) {
    674 
    675     bool status;
    676     char keyword[80];
    677 
    678     // get polynomial orders
    679     // XXX add status failures tests
    680     int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
    681     int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
    682     int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
    683 
    684     psPolynomial3D *poly = psPolynomial3DAlloc (nXorder, nYorder, nZorder, PS_POLYNOMIAL_ORD);
    685 
    686     for (int nx = 0; nx < poly->nX + 1; nx++) {
    687         for (int ny = 0; ny < poly->nY + 1; ny++) {
    688             for (int nz = 0; nz < poly->nZ + 1; nz++) {
    689                 sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d", nx, ny, nz);
    690                 poly->coeff[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
    691                 if (!status) poly->mask[nx][ny][nz] = 1;
    692 
    693                 sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d", nx, ny, nz);
    694                 poly->coeffErr[nx][ny][nz] = psMetadataLookupF32 (&status, folder, keyword);
    695             }
    696         }
    697     }
    698     return (poly);
    699 }
    700 
    701 psPolynomial4D *psPolynomial4DfromMD (psMetadata *folder) {
    702 
    703     bool status;
    704     char keyword[80];
    705 
    706     // get polynomial orders
    707     // XXX add status failures tests
    708     int nXorder = psMetadataLookupS32 (&status, folder, "NORDER_X");
    709     int nYorder = psMetadataLookupS32 (&status, folder, "NORDER_Y");
    710     int nZorder = psMetadataLookupS32 (&status, folder, "NORDER_Z");
    711     int nTorder = psMetadataLookupS32 (&status, folder, "NORDER_T");
    712 
    713     psPolynomial4D *poly = psPolynomial4DAlloc (nXorder, nYorder, nZorder, nTorder, PS_POLYNOMIAL_ORD);
    714 
    715     for (int nx = 0; nx < poly->nX + 1; nx++) {
    716         for (int ny = 0; ny < poly->nY + 1; ny++) {
    717             for (int nz = 0; nz < poly->nZ + 1; nz++) {
    718                 for (int nt = 0; nt < poly->nT + 1; nt++) {
    719                     sprintf (keyword, "VAL_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
    720                     poly->coeff[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
    721                     if (!status) poly->mask[nx][ny][nz][nt] = 1;
    722                    
    723                     sprintf (keyword, "ERR_X%02d_Y%02d_Z%02d_T%02d", nx, ny, nz, nt);
    724                     poly->coeffErr[nx][ny][nz][nt] = psMetadataLookupF32 (&status, folder, keyword);
    725                 }
    726             }
    727         }
    728     }
    729     return (poly);
    730 }
    731 
    732 // XXX : these may need F64, or %g format for output
    733 bool psPolynomial2DtoMD (psMetadata *md, psPolynomial2D *poly, char *format, ...) {
    734 
    735     int Nbyte;
    736     char tmp;
    737     char *root;
    738     va_list argp; 
    739 
    740     va_start (argp, format);
    741     Nbyte = vsnprintf (&tmp, 0, format, argp);
    742     va_end (argp);
    743 
    744     if (!Nbyte) return false;
    745 
    746     va_start (argp, format);
    747     root = (char *) psAlloc (Nbyte + 1);
    748     memset (root, 0, Nbyte + 1);
    749     vsnprintf (root, Nbyte + 1, format, argp);
    750     va_end (argp);
    751 
    752     psMetadata *folder = psMetadataAlloc ();
    753     psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 2D polynomial", folder);
    754     psFree (root);
    755 
    756     // specify the polynomial orders
    757     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
    758     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
    759 
    760     // place polynomial entries on folder
    761     for (int nx = 0; nx < poly->nX + 1; nx++) {
    762         for (int ny = 0; ny < poly->nY + 1; ny++) {
    763             if (poly->mask[nx][ny]) continue;
    764             psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny], nx, ny);
    765             psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d", PS_DATA_F32, "polynomial coefficient error", poly->coeffErr[nx][ny], nx, ny);
    766         }
    767     }
    768     return true;
    769 }
    770 
    771 bool psPolynomial3DtoMD (psMetadata *md, psPolynomial3D *poly, char *format, ...) {
    772 
    773     int Nbyte;
    774     char tmp;
    775     char *root;
    776     va_list argp; 
    777 
    778     va_start (argp, format);
    779     Nbyte = vsnprintf (&tmp, 0, format, argp);
    780     va_end (argp);
    781 
    782     if (!Nbyte) return false;
    783 
    784     va_start (argp, format);
    785     root = (char *) psAlloc (Nbyte + 1);
    786     memset (root, 0, Nbyte + 1);
    787     vsnprintf (root, Nbyte + 1, format, argp);
    788     va_end (argp);
    789 
    790     psMetadata *folder = psMetadataAlloc ();
    791     psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 3D polynomial", folder);
    792     psFree (root);
    793 
    794     // specify the polynomial orders
    795     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
    796     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
    797     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
    798 
    799     // place polynomial entries on folder
    800     for (int nx = 0; nx < poly->nX + 1; nx++) {
    801         for (int ny = 0; ny < poly->nY + 1; ny++) {
    802             for (int nz = 0; nz < poly->nZ + 1; nz++) {
    803                 if (poly->mask[nx][ny][nz]) continue;
    804                 psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz], nx, ny, nz);
    805                 psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz], nx, ny, nz);
    806             }
    807         }
    808     }
    809     return true;
    810 }
    811 
    812 bool psPolynomial4DtoMD (psMetadata *md, psPolynomial4D *poly, char *format, ...) {
    813 
    814     int Nbyte;
    815     char tmp;
    816     char *root;
    817     va_list argp; 
    818 
    819     va_start (argp, format);
    820     Nbyte = vsnprintf (&tmp, 0, format, argp);
    821     va_end (argp);
    822 
    823     if (!Nbyte) return false;
    824 
    825     va_start (argp, format);
    826     root = (char *) psAlloc (Nbyte + 1);
    827     memset (root, 0, Nbyte + 1);
    828     vsnprintf (root, Nbyte + 1, format, argp);
    829     va_end (argp);
    830 
    831     psMetadata *folder = psMetadataAlloc ();
    832     psMetadataAdd (md, PS_LIST_TAIL, root, PS_DATA_METADATA, "folder for 4D polynomial", folder);
    833     psFree (root);
    834 
    835     // specify the polynomial orders
    836     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_X", PS_DATA_S32, "number of x orders", poly->nX);
    837     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Y", PS_DATA_S32, "number of y orders", poly->nY);
    838     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_Z", PS_DATA_S32, "number of z orders", poly->nZ);
    839     psMetadataAdd (folder, PS_LIST_TAIL, "NORDER_T", PS_DATA_S32, "number of z orders", poly->nT);
    840 
    841     // place polynomial entries on folder
    842     for (int nx = 0; nx < poly->nX + 1; nx++) {
    843         for (int ny = 0; ny < poly->nY + 1; ny++) {
    844             for (int nz = 0; nz < poly->nZ + 1; nz++) {
    845                 for (int nt = 0; nt < poly->nT + 1; nt++) {
    846                     if (poly->mask[nx][ny][nz][nt]) continue;
    847                     psMetadataAdd (folder, PS_LIST_TAIL, "VAL_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coefficient", poly->coeff[nx][ny][nz][nt], nx, ny, nz, nt);
    848                     psMetadataAdd (folder, PS_LIST_TAIL, "ERR_X%02d_Y%02d_Z%02d_T%02d", PS_DATA_F32, "polynomial coeffficient error", poly->coeffErr[nx][ny][nz][nt], nx, ny, nz, nt);
    849                 }
    850             }
    851         }
    852     }
    853     return true;
    854 }
    855 
    856 bool psphotWritePSF (pmPSF *psf, char *filename) {
    857 
    858     psMetadata *psfdata = psMetadataAlloc ();
    859 
    860     char *modelName = pmModelGetType (psf->type);
    861     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName);
    862 
    863     int nPar = pmModelParameterCount (psf->type)    ;
    864     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
    865 
    866     for (int i = 0; i < nPar - 4; i++) {
    867         psPolynomial2D *poly = psf->params->data[i];
    868         psPolynomial2DtoMD (psfdata, poly, "PSF_PAR%02d", i);
    869     }
    870     psPolynomial4DtoMD (psfdata, psf->ApTrend, "APTREND");
    871 
    872     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid);
    873     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid);
    874     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
    875 
    876     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq);
    877     psMetadataAdd (psfdata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars);
    878    
    879     psMetadataConfigWrite (psfdata, filename);
    880     psFree (psfdata);
    881 
    882     return true;
    883 }
    884 
    885 pmPSF *psphotReadPSF (char *filename) {
    886 
    887     bool status;
    888     unsigned int Nfail;
    889     char keyword[80];
    890 
    891     psMetadata *psfdata = psMetadataConfigParse (NULL, &Nfail, filename, FALSE);
    892 
    893     char *modelName = psMetadataLookupPtr (&status, psfdata, "PSF_MODEL_NAME");
    894     pmModelType type = pmModelSetType (modelName);
    895 
    896     pmPSF *psf = pmPSFAlloc (type);
    897 
    898     int nPar = psMetadataLookupS32 (&status, psfdata, "PSF_MODEL_NPAR");
    899     if (nPar != pmModelParameterCount (psf->type)) psAbort ("read PSF" , "mismatch model par count");
    900 
    901     for (int i = 0; i < nPar - 4; i++) {
    902         sprintf (keyword, "PSF_PAR%02d", i);
    903         psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
    904         psPolynomial2D *poly = psPolynomial2DfromMD (folder);
    905         psFree (psf->params->data[i]);
    906         psf->params->data[i] = poly;
    907     }
    908     sprintf (keyword, "APTREND");
    909     psMetadata *folder = psMetadataLookupPtr (&status, psfdata, keyword);
    910     psPolynomial4D *poly = psPolynomial4DfromMD (folder);
    911     psFree (psf->ApTrend);
    912     psf->ApTrend = poly;
    913 
    914     psf->ApResid = psMetadataLookupF32 (&status, psfdata, "PSF_AP_RESID");
    915     psf->dApResid = psMetadataLookupF32 (&status, psfdata, "PSF_dAP_RESID");
    916     psf->skyBias = psMetadataLookupF32 (&status, psfdata, "PSF_SKY_BIAS");
    917 
    918     psf->chisq = psMetadataLookupF32 (&status, psfdata, "PSF_CHISQ");
    919     psf->nPSFstars = psMetadataLookupS32 (&status, psfdata, "PSF_NSTARS");
    920 
    921     psFree (psfdata);
    922     return (psf);
    923614}
    924615
    … …  
    937628    return true;
    938629}
     630
     631bool psphotUpdateHeader (psMetadata *header, psMetadata *config) {
     632
     633    // write necessary information to output header
     634    psMetadataItemTransfer (header, config, "NSTARS");
     635    psMetadataItemTransfer (header, config, "ZERO_PT");
     636    psMetadataItemTransfer (header, config, "APMIFIT");
     637    psMetadataItemTransfer (header, config, "dAPMIFIT");
     638    psMetadataItemTransfer (header, config, "SKYBIAS");
     639    psMetadataItemTransfer (header, config, "PHOTCODE");
     640   
     641    // write necessary information to output header
     642    psMetadataItemTransfer (header, config, "NPSFSTAR");
     643    psMetadataItemTransfer (header, config, "SKYBIAS");
     644    psMetadataItemTransfer (header, config, "SKYSAT");
     645    psMetadataItemTransfer (header, config, "APMIFIT");
     646    psMetadataItemTransfer (header, config, "DAPMIFIT");
     647    psMetadataItemTransfer (header, config, "NAPMIFIT");
     648    psMetadataItemTransfer (header, config, "APLOSS");
     649
     650    // XXX these need to be defined from elsewhere
     651    psMetadataAdd (header, PS_LIST_TAIL, "FWHM_X",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM X",          0.0);
     652    psMetadataAdd (header, PS_LIST_TAIL, "FWHM_Y",   PS_DATA_F32 | PS_META_REPLACE, "PSF FWHM Y",          0.0);
     653    psMetadataAdd (header, PS_LIST_TAIL, "ANGLE",    PS_DATA_F32 | PS_META_REPLACE, "PSF ANGLE",           0.0);
     654    psMetadataAdd (header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
     655    psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
     656
     657    return true;
     658}
Note: See TracChangeset for help on using the changeset viewer.