IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 4, 2021, 6:05:18 PM (5 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-dev-20210817 (pmPattern updates, new inverse transform extra orders api, forward transform uses ORD, pmSourceIO_CMF.c.in conversion to pmFitsTableNew)

Location:
trunk/psModules
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules

  • trunk/psModules/src/objects/pmSourceIO_CMF.c.in

    r38114 r41892  
    5757// followed by a zero-size matrix, followed by the table data
    5858
    59 bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
     59/** this file is used to generate c-code for a number of different formats.  the perl program
     60 ** 'mksource.pl' reads this template file and generates an output file for one of the named
     61 ** formats, e.g., PS1_V1 or PS1_DV4 (see list in mksource.pl). 
     62
     63 ** Any instances of the word @CMFMODE@ are replaced by the named version of the format.
     64
     65 ** Any line which starts with a command of the form: @COMMAND[,COMMAND...]@ is kept or removed
     66 ** from the output c-code depending on the commands.  These may be: ALL (lines are kept for
     67 ** all formats) or the name of a format.  A line with just a plain format will be kept only
     68 ** for that format.  The format name may be optionally proceded by >, <, <=, >=, ! or =.  In
     69 ** these cases, the line is kept if the format is greater / lesser / not the same, etc as the
     70 ** listed format.  The sequence of formats is defined for series (V, SV, DV) in mksource.pl so
     71 ** that, e.g., >PS1_V3 would include e.g., PS1_V4 and PS1_V5, but not PS1_DV4.  The number in
     72 ** a series may be replaced with ? to stand for the full set.  The format
     73 ** commands may be grouped on a line separated by commas.  Note that only ! (not) is
     74 ** restrictive: only this rule will reduce the number of matches.  So for example,
     75 ** @ALL,<PS1_DV3@ will have the effect of @ALL@ since any format >= PS1_DV3 will be matched by
     76 ** ALL.  But @ALL,!PS1_DV3@ will generate a line for any format except PS1_DV3.
     77
     78 **/
     79
     80bool pmSourcesWrite_CMF_@CMFMODE@_New (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
    6081{
     82    // fprintf (stderr, "writing with %s\n", __func__);
    6183    PS_ASSERT_PTR_NON_NULL(fits, false);
    6284    PS_ASSERT_PTR_NON_NULL(sources, false);
    6385    PS_ASSERT_PTR_NON_NULL(extname, false);
    64 
    65     psArray *table;
    66     psMetadata *row;
    6786
    6887    pmChip *chip = readout->parent->parent;
    … …  
    7392        pmSource *source = (pmSource *) sources->data[0];
    7493        if (source->seq == -1) {
    75             // let's write these out in S/N order
    76             sources = psArraySort (sources, pmSourceSortByFlux);
     94            // let's write these out in S/N order
     95            sources = psArraySort (sources, pmSourceSortByFlux);
    7796        } else {
    78             sources = psArraySort (sources, pmSourceSortBySeq);
    79         }
    80     }
    81 
    82     table = psArrayAllocEmpty (sources->n);
     97            sources = psArraySort (sources, pmSourceSortBySeq);
     98        }
     99    }
    83100
    84101    float magOffset;
    … …  
    88105    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
    89106
     107    // before we can generate the table structure, we need to know if certain fields were measured.
     108    // this information is not available a priori, so we need to check all sources
     109   
     110    bool haveLensingOBJsmear = false;
     111    bool haveLensingOBJshear = false;
     112    bool haveLensingPSFsmear = false;
     113    bool haveLensingPSFshear = false;
     114    bool haveLensingPSF      = false;
     115
     116    for (int i = 0; (i < sources->n) && !haveLensingOBJsmear && !haveLensingOBJshear && !haveLensingPSFsmear && !haveLensingPSFshear && !haveLensingPSF; i++) {
     117        // this is the source associated with this image
     118        pmSource *thisSource = sources->data[i];
     119
     120        // this is the "real" version of this source
     121        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     122
     123        if (source->lensingOBJ && source->lensingOBJ->smear) haveLensingOBJsmear = true;
     124        if (source->lensingOBJ && source->lensingOBJ->shear) haveLensingOBJshear = true;
     125        if (source->lensingPSF && source->lensingPSF->smear) haveLensingPSFsmear = true;
     126        if (source->lensingPSF && source->lensingPSF->shear) haveLensingPSFshear = true;
     127        if (source->lensingPSF                             ) haveLensingPSF      = true;
     128    }
     129
     130    /************ generate the table columns *****************/
     131
     132    // before we allocate the table, generate an empty array of table columns and generate them
     133    psArray *tableColumns = psArrayAllocEmpty (100);
     134
     135    // add the named / typed columns to the collection of columns
     136    @ALL@                       psFitsTableColumnAdd (tableColumns, "IPP_IDET",         PS_DATA_U32); // "IPP detection identifier index"
     137    @ALL@                       psFitsTableColumnAdd (tableColumns, "X_PSF",            PS_DATA_F32); // "PSF x coordinate"
     138    @ALL@                       psFitsTableColumnAdd (tableColumns, "Y_PSF",            PS_DATA_F32); // "PSF y coordinate"
     139    @ALL@                       psFitsTableColumnAdd (tableColumns, "X_PSF_SIG",        PS_DATA_F32); // "Sigma in PSF x coordinate"
     140    @ALL@                       psFitsTableColumnAdd (tableColumns, "Y_PSF_SIG",        PS_DATA_F32); // "Sigma in PSF y coordinate"
     141
     142    // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
     143    @PS1_V1@                    psFitsTableColumnAdd (tableColumns, "RA_PSF",           PS_DATA_F32); // "PSF RA coordinate (degrees)"
     144    @PS1_V1@                    psFitsTableColumnAdd (tableColumns, "DEC_PSF",          PS_DATA_F32); // "PSF DEC coordinate (degrees)"
     145
     146    @ALL@                       psFitsTableColumnAdd (tableColumns, "POSANGLE",         PS_DATA_F32); // "position angle at source (degrees)"
     147    @ALL@                       psFitsTableColumnAdd (tableColumns, "PLTSCALE",         PS_DATA_F32); // "plate scale at source (arcsec/pixel)"
     148    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG",     PS_DATA_F32); // "PSF fit instrumental magnitude"
     149    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG_SIG", PS_DATA_F32); // "Sigma of PSF instrumental magnitude"
     150
     151    @ALL,!PS1_V1,!PS1_V2@       psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX",    PS_DATA_F32); // "PSF fit instrumental flux (counts)"
     152    @ALL,!PS1_V1,!PS1_V2@       psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX_SIG", PS_DATA_F32); // "Sigma of PSF instrumental flux"
     153
     154    @ALL@                       psFitsTableColumnAdd (tableColumns, "AP_MAG",           PS_DATA_F32); // "magnitude in standard aperture"
     155    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "AP_MAG_RAW",       PS_DATA_F32); // "magnitude in reported aperture"
     156    @ALL@                       psFitsTableColumnAdd (tableColumns, "AP_MAG_RADIUS",    PS_DATA_F32); // "radius used for aperture mags"
     157    @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX",          PS_DATA_F32); // "instrumental flux in standard aperture"
     158    @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX_SIG",      PS_DATA_F32); // "aperture flux error"
     159    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "AP_NPIX",          PS_DATA_S32); // "aperture unmasked pixels"
     160
     161    @<PS1_V3,PS1_SV1,PS1_DV?@   psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude"
     162
     163    @ALL@                       psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG",      PS_DATA_F32); // "PSF Magnitude using supplied calibration"
     164    @ALL@                       psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG_SIG",  PS_DATA_F32); // "measured scatter of zero point calibration"
     165
     166    // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
     167    @ALL,!PS1_V1@               psFitsTableColumnAdd (tableColumns, "RA_PSF",           PS_DATA_F64); // "PSF RA coordinate (degrees)"
     168    @ALL,!PS1_V1@               psFitsTableColumnAdd (tableColumns, "DEC_PSF",          PS_DATA_F64); // "PSF DEC coordinate (degrees)"
     169
     170    @>=PS1_V3,>PS1_SV1@         psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude"
     171    @ALL@                       psFitsTableColumnAdd (tableColumns, "SKY",              PS_DATA_F32); // "Sky level"
     172    @ALL@                       psFitsTableColumnAdd (tableColumns, "SKY_SIGMA",        PS_DATA_F32); // "Sigma of sky level"
     173    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_CHISQ",        PS_DATA_F32); // "Chisq of PSF-fit"
     174    @ALL@                       psFitsTableColumnAdd (tableColumns, "CR_NSIGMA",        PS_DATA_F32); // "Nsigma deviations from PSF to CF"
     175    @ALL@                       psFitsTableColumnAdd (tableColumns, "EXT_NSIGMA",       PS_DATA_F32); // "Nsigma deviations from PSF to EXT"
     176
     177    // PSF shape parameters:
     178    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_MAJOR",        PS_DATA_F32); // "PSF width (major axis)"
     179    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_MINOR",        PS_DATA_F32); // "PSF width (minor axis)"
     180    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_THETA",        PS_DATA_F32); // "PSF orientation angle"
     181    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_CORE",         PS_DATA_F32); // "k term if defined"
     182    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MAJ",     PS_DATA_F32); // "PSF FWHM (major axis)"
     183    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MIN",     PS_DATA_F32); // "PSF FWHM (minor axis)"
     184
     185    // psf data quality
     186    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_QF",           PS_DATA_F32); // "PSF coverage/quality factor (bad)"
     187    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "PSF_QF_PERFECT",   PS_DATA_F32); // "PSF coverage/quality factor (poor)"
     188    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_NDOF",         PS_DATA_S32); // "degrees of freedom"
     189    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_NPIX",         PS_DATA_S32); // "number of pixels in fit"
     190
     191    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_XX",       PS_DATA_F32); // "second moments (X^2)"
     192    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_XY",       PS_DATA_F32); // "second moments (X*Y)"
     193    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_YY",       PS_DATA_F32); // "second moments (Y*Y)"
     194    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M3C",      PS_DATA_F32); // "third momemt cos theta"
     195    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M3S",      PS_DATA_F32); // "third momemt sin theta"
     196    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M4C",      PS_DATA_F32); // "fourth momemt cos theta"
     197    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M4S",      PS_DATA_F32); // "fourth momemt sin theta"
     198 
     199    // Lensing parameters are only written if they are measured
     200    if (haveLensingOBJsmear) {
     201        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
     202        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
     203        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
     204        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_OBJ",   PS_DATA_F32); // "smear polarizability element (objects)"
     205        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_OBJ",   PS_DATA_F32); // "smear polarizability element (objects)"
     206    }
     207    if (haveLensingOBJshear) {
     208        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
     209        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
     210        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
     211        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_OBJ",   PS_DATA_F32); // "shear polarizability element (objects)"
     212        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_OBJ",   PS_DATA_F32); // "shear polarizability element (objects)"
     213    }
     214    if (false) {
     215        // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
     216        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_OBJ",      PS_DATA_F32); // "shear polarizability element (objects)"
     217        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_OBJ",      PS_DATA_F32); // "shear polarizability element (objects)"
     218    }
     219
     220    if (haveLensingPSFsmear) {
     221        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
     222        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
     223        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
     224        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_PSF",   PS_DATA_F32); // "smear polarizability element (PSFs)"
     225        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_PSF",   PS_DATA_F32); // "smear polarizability element (PSFs)"
     226    }
     227    if (haveLensingPSFshear) {
     228        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
     229        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
     230        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
     231        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_PSF",   PS_DATA_F32); // "shear polarizability element (PSFs)"
     232        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_PSF",   PS_DATA_F32); // "shear polarizability element (PSFs)"
     233    }
     234    if (haveLensingPSF) {
     235        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_PSF",      PS_DATA_F32); // "shear polarizability element (PSFs)"
     236        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_PSF",      PS_DATA_F32); // "shear polarizability element (PSFs)"
     237    }
     238
     239    // if lensing params exist also include the backmapped chipID and chip coordinates
     240    if (haveLensingPSF) {
     241        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM",     PS_DATA_S16); // "id of warp input chip"
     242        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X",       PS_DATA_S16); // "x coord in warp input chip"
     243        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y",       PS_DATA_S16); // "y coord in warp input chip"
     244        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "PADDING3",         PS_DATA_S16); // "more padding"
     245    }
     246
     247    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "MOMENTS_R1",       PS_DATA_F32); // "first radial moment"
     248    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "MOMENTS_RH",       PS_DATA_F32); // "half radial moment"
     249    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX",        PS_DATA_F32); // "Kron Flux (in 2.5 R1)"
     250    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_ERR",    PS_DATA_F32); // "Kron Flux Error"
     251    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_INNER",  PS_DATA_F32); // "Kron Flux (in 1.0 R1)"
     252    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_OUTER",  PS_DATA_F32); // "Kron Flux (in 4.0 R1)"
     253
     254    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_RAD",    PS_DATA_F32); // "Radius where object hits sky"
     255    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_FLUX",   PS_DATA_F32); // "Flux / pix where object hits sky"
     256    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_SLOPE",  PS_DATA_F32); // "d(Flux/pix)/dRadius where object hits sky"
     257
     258    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM",     PS_DATA_S16); // "id of warp input chip"
     259    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X",       PS_DATA_S16); // "x coord in warp input chip"
     260    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y",       PS_DATA_S16); // "y coord in warp input chip"
     261    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "PADDING3",         PS_DATA_S16); // "more padding"
     262
     263    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NPOS",        PS_DATA_S32); // "nPos (n pix > 3 sigma)"
     264    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_FRATIO",      PS_DATA_F32); // "fPos / (fPos + fNeg)"
     265    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_BAD",  PS_DATA_F32); // "nPos / (nPos + nNeg)"
     266    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_MASK", PS_DATA_F32); // "nPos / (nPos + nMask)"
     267    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_ALL",  PS_DATA_F32); // "nPos / (nGood + nMask + nBad)"
     268
     269    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_R_P",         PS_DATA_F32); // "distance to positive match source"
     270    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_SN_P",        PS_DATA_F32); // "signal-to-noise of pos match src"
     271    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_R_M",         PS_DATA_F32); // "distance to negative match source"
     272    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_SN_M",        PS_DATA_F32); // "signal-to-noise of neg match src"
     273
     274    @ALL@                       psFitsTableColumnAdd (tableColumns, "FLAGS",            PS_DATA_U32); // "psphot analysis flags"
     275    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "FLAGS2",           PS_DATA_U32); // "psphot analysis flags"
     276    @PS1_V3,PS1_SV2@            psFitsTableColumnAdd (tableColumns, "PADDING2",         PS_DATA_S32); // "more padding"
     277    @ALL@                       psFitsTableColumnAdd (tableColumns, "N_FRAMES",         PS_DATA_U16); // "Number of frames overlapping source center"
     278
     279    @ALL@                       psFitsTableColumnAdd (tableColumns, "PADDING",          PS_DATA_S16); // "more padding"
     280
     281    // note that this definition is inconsistent with the definition in
     282    // Ohana/src/libautocode/def.  This version creates a table with data not
     283    // properly aligned with the 8-byte boundaries. The structure defined by
     284    // libautocode does this, but has a different order of elements (and adds
     285    // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
     286    // I'll leave it. But in future a PS1_SV2 should be forced to match
     287    // libautocode. Note that addstar knows to detect the alternate version of
     288    // PS1_SV1 and correctly interpret its fields.
     289
     290    // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
     291    // subtracted
     292
     293    // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
     294    // edge; 3) any pixels in the 3x3 peak region are masked;
     295
     296    // generate an FITS table using the array of columns defined above, with space for all rows
     297    psFitsTable *table = psFitsTableCreate (tableColumns, sources->n);
     298
     299    /************ write the data to the table *****************/
     300
    90301    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
    91302    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
    92303    for (int i = 0; i < sources->n; i++) {
    93         // this is the source associated with this image
     304        // this is the source associated with this image
    94305        pmSource *thisSource = sources->data[i];
    95306
    96         // this is the "real" version of this source
    97         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     307        // this is the "real" version of this source
     308        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    98309
    99310        // If source->seq is -1, source was generated in this analysis.  If source->seq is
    … …  
    102313        // generated on Alloc, and would thus be wrong for read in sources.
    103314        if (source->seq == -1) {
    104             source->seq = i;
    105         }
    106 
    107         // set the 'best' values for various output fields:
    108         pmSourceOutputs outputs;
    109         pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
    110 
    111         pmSourceOutputsMoments moments;
    112         pmSourceOutputsSetMoments (&moments, source);
    113 
    114         @PS1_DV?@ pmSourceDiffStats diffStats;
    115         @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
    116         @PS1_DV?@ if (source->diffStats) {
    117         @PS1_DV?@     diffStats = *source->diffStats;
    118         @PS1_DV?@ }
    119 
    120         row = psMetadataAlloc ();
     315            source->seq = i;
     316        }
     317
     318        // set the 'best' values for various output fields:
     319        pmSourceOutputs outputs;
     320        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     321
     322        pmSourceOutputsMoments moments;
     323        pmSourceOutputsSetMoments (&moments, source);
     324
     325        @PS1_DV?@ pmSourceDiffStats diffStats;
     326        @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
     327        @PS1_DV?@ if (source->diffStats) { diffStats = *source->diffStats;}
    121328
    122329        // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
    123330        // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment"
    124         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    125         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
    126         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
    127         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
    128         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
    129 
    130         // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
    131         @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
    132         @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
    133 
    134         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
    135         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    136         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    137         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
    138 
    139         @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
    140         @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
    141 
    142         @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    143         @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    144         @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              source->apRadius);
    145         @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
    146         @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
    147         @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX",          PS_DATA_S32, "aperture unmasked pixels",                   source->apNpixels);
    148 
    149         @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    150 
    151         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    152         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    153        
    154         // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
    155         @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
    156         @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    157 
    158         @>=PS1_V3,>PS1_SV1@       psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    159         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    160         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    161 
    162         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    163         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    164         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
     331        @ALL@                       psFitsTableSetU32 (table, i, "IPP_IDET",          source->seq);                     // "IPP detection identifier index",           
     332        @ALL@                       psFitsTableSetF32 (table, i, "X_PSF",             outputs.xPos);                    // "PSF x coordinate",                         
     333        @ALL@                       psFitsTableSetF32 (table, i, "Y_PSF",             outputs.yPos);                    // "PSF y coordinate",                         
     334        @ALL@                       psFitsTableSetF32 (table, i, "X_PSF_SIG",         outputs.xErr);                    // "Sigma in PSF x coordinate",                 
     335        @ALL@                       psFitsTableSetF32 (table, i, "Y_PSF_SIG",         outputs.yErr);                    // "Sigma in PSF y coordinate",                 
     336
     337        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision                           
     338        @PS1_V1@                    psFitsTableSetF32 (table, i, "RA_PSF",            outputs.ra);                      // "PSF RA coordinate (degrees)",               
     339        @PS1_V1@                    psFitsTableSetF32 (table, i, "DEC_PSF",           outputs.dec);                     // "PSF DEC coordinate (degrees)",               
     340
     341        @ALL@                       psFitsTableSetF32 (table, i, "POSANGLE",          outputs.posAngle);                // "position angle at source (degrees)",       
     342        @ALL@                       psFitsTableSetF32 (table, i, "PLTSCALE",          outputs.pltScale);                // "plate scale at source (arcsec/pixel)",     
     343        @ALL@                       psFitsTableSetF32 (table, i, "PSF_INST_MAG",      source->psfMag);                  // "PSF fit instrumental magnitude",           
     344        @ALL@                       psFitsTableSetF32 (table, i, "PSF_INST_MAG_SIG",  source->psfMagErr);               // "Sigma of PSF instrumental magnitude",       
     345
     346        @ALL,!PS1_V1,!PS1_V2@       psFitsTableSetF32 (table, i, "PSF_INST_FLUX",     source->psfFlux);                 // "PSF fit instrumental flux (counts)",       
     347        @ALL,!PS1_V1,!PS1_V2@       psFitsTableSetF32 (table, i, "PSF_INST_FLUX_SIG", source->psfFluxErr);              // "Sigma of PSF instrumental flux",           
     348
     349        @ALL@                       psFitsTableSetF32 (table, i, "AP_MAG",            source->apMag);                   // "magnitude in standard aperture",           
     350        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "AP_MAG_RAW",        source->apMagRaw);                // "magnitude in reported aperture",           
     351        @ALL@                       psFitsTableSetF32 (table, i, "AP_MAG_RADIUS",     source->apRadius);                // "radius used for aperture mags",             
     352        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX",           source->apFlux);                  // "instrumental flux in standard aperture",   
     353        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX_SIG",       source->apFluxErr);               // "aperture flux error",                       
     354        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetS32 (table, i, "AP_NPIX",           source->apNpixels);               // "aperture unmasked pixels",                 
     355
     356        @<PS1_V3,PS1_SV1,PS1_DV?@   psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG",  outputs.peakMag);                 // "Peak flux expressed as magnitude"
     357
     358        @ALL@                       psFitsTableSetF32 (table, i, "CAL_PSF_MAG",       outputs.calMag);                  // "PSF Magnitude using supplied calibration", 
     359        @ALL@                       psFitsTableSetF32 (table, i, "CAL_PSF_MAG_SIG",   zeroptErr);                       // "measured scatter of zero point calibration",
     360
     361        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...                                     
     362        @ALL,!PS1_V1@               psFitsTableSetF64 (table, i, "RA_PSF",            outputs.ra);                      // "PSF RA coordinate (degrees)",               
     363        @ALL,!PS1_V1@               psFitsTableSetF64 (table, i, "DEC_PSF",           outputs.dec);                     // "PSF DEC coordinate (degrees)",             
     364
     365        @>=PS1_V3,>PS1_SV1@         psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG",  outputs.peakMag);                 // "Peak flux expressed as magnitude",         
     366        @ALL@                       psFitsTableSetF32 (table, i, "SKY",               source->sky);                     // "Sky level",                                 
     367        @ALL@                       psFitsTableSetF32 (table, i, "SKY_SIGMA",         source->skyErr);                  // "Sigma of sky level",                       
     368
     369        @ALL@                       psFitsTableSetF32 (table, i, "PSF_CHISQ",         outputs.chisq);                   // "Chisq of PSF-fit",                         
     370        @ALL@                       psFitsTableSetF32 (table, i, "CR_NSIGMA",         source->crNsigma);                // "Nsigma deviations from PSF to CF",         
     371        @ALL@                       psFitsTableSetF32 (table, i, "EXT_NSIGMA",        source->extNsigma);               // "Nsigma deviations from PSF to EXT",         
    165372
    166373        // PSF shape parameters:
    167         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
    168         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
    169         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
     374        @ALL@                       psFitsTableSetF32 (table, i, "PSF_MAJOR",         outputs.psfMajor);                // "PSF width (major axis)",                   
     375        @ALL@                       psFitsTableSetF32 (table, i, "PSF_MINOR",         outputs.psfMinor);                // "PSF width (minor axis)",                   
     376        @ALL@                       psFitsTableSetF32 (table, i, "PSF_THETA",         outputs.psfTheta);                // "PSF orientation angle",                     
     377        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_CORE",          outputs.psfCore);                 // "k term if defined",                         
     378
     379        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
     380        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MAJ",      outputs.psfMajorFWHM);            // "PSF FWHM (major axis)",                     
     381        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MIN",      outputs.psfMinorFWHM);            // "PSF FWHM (minor axis)",                     
     382
     383        // psf data quality
     384        @ALL@                       psFitsTableSetF32 (table, i, "PSF_QF",            source->pixWeightNotBad);         // "PSF coverage/quality factor (bad)",         
     385        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "PSF_QF_PERFECT",    source->pixWeightNotPoor);        // "PSF coverage/quality factor (poor)",       
     386        @ALL@                       psFitsTableSetS32 (table, i, "PSF_NDOF",          outputs.nDOF);                    // "degrees of freedom",                       
     387        @ALL@                       psFitsTableSetS32 (table, i, "PSF_NPIX",          outputs.nPix);                    // "number of pixels in fit",                   
     388
     389        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_XX",        moments.Mxx);                     // "second moments (X^2)",                     
     390        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_XY",        moments.Mxy);                     // "second moments (X*Y)",                     
     391        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_YY",        moments.Myy);                     // "second moments (Y*Y)",                     
     392
     393        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M3C",       moments.M_c3);                    // "third momemt cos theta",                   
     394        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M3S",       moments.M_s3);                    // "third momemt sin theta",                   
     395        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M4C",       moments.M_c4);                    // "fourth momemt cos theta",                   
     396        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M4S",       moments.M_s4);                    // "fourth momemt sin theta",                   
     397
     398        // Lensing parameters:
     399        if (source->lensingOBJ && source->lensingOBJ->smear) {
     400            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SM_OBJ",   source->lensingOBJ->smear->X11);  // "smear polarizability element (objects)",     
     401            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SM_OBJ",   source->lensingOBJ->smear->X12);  // "smear polarizability element (objects)",     
     402            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SM_OBJ",   source->lensingOBJ->smear->X22);  // "smear polarizability element (objects)",     
     403            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SM_OBJ",    source->lensingOBJ->smear->e1);   // "smear polarizability element (objects)",     
     404            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SM_OBJ",    source->lensingOBJ->smear->e2);   // "smear polarizability element (objects)",     
     405        }
     406        if (source->lensingOBJ && source->lensingOBJ->shear) {
     407            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SH_OBJ",   source->lensingOBJ->shear->X11);  // "shear polarizability element (objects)",     
     408            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SH_OBJ",   source->lensingOBJ->shear->X12);  // "shear polarizability element (objects)",     
     409            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SH_OBJ",   source->lensingOBJ->shear->X22);  // "shear polarizability element (objects)",     
     410            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SH_OBJ",    source->lensingOBJ->shear->e1);   // "shear polarizability element (objects)",     
     411            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SH_OBJ",    source->lensingOBJ->shear->e2);   // "shear polarizability element (objects)",     
     412        }
     413        if (false && source->lensingOBJ) {
     414            // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
     415            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_PSF",       source->lensingOBJ->e1);          // "shear polarizability element (objects)",     
     416            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_PSF",       source->lensingOBJ->e2);          // "shear polarizability element (objects)",     
     417        }
     418
     419        if (source->lensingPSF && source->lensingPSF->smear) {
     420            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SM_PSF",   source->lensingPSF->smear->X11);  // "smear polarizability element (objects)",     
     421            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SM_PSF",   source->lensingPSF->smear->X12);  // "smear polarizability element (objects)",     
     422            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SM_PSF",   source->lensingPSF->smear->X22);  // "smear polarizability element (objects)",     
     423            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SM_PSF",    source->lensingPSF->smear->e1);   // "smear polarizability element (objects)",     
     424            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SM_PSF",    source->lensingPSF->smear->e2);   // "smear polarizability element (objects)",     
     425        }
     426        if (source->lensingPSF && source->lensingPSF->shear) {
     427            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SH_PSF",   source->lensingPSF->shear->X11);  // "shear polarizability element (objects)",     
     428            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SH_PSF",   source->lensingPSF->shear->X12);  // "shear polarizability element (objects)",     
     429            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SH_PSF",   source->lensingPSF->shear->X22);  // "shear polarizability element (objects)",     
     430            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SH_PSF",    source->lensingPSF->shear->e1);   // "shear polarizability element (objects)",     
     431            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SH_PSF",    source->lensingPSF->shear->e2);   // "shear polarizability element (objects)",     
     432        }
     433        if (source->lensingPSF) {
     434            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_PSF",       source->lensingPSF->e1);          // "shear polarizability element (objects)",     
     435            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_PSF",       source->lensingPSF->e2);          // "shear polarizability element (objects)",     
     436        }
     437
     438        // if lensing params exist also include the backmapped chipID and chip coordinates
     439        if (source->lensingPSF && source->lensingPSF->shear) {
     440            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_NUM",      source->chipNum);                 // "id of warp input chip",       
     441            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_X",        source->chipX);                   // "x coord in warp input chip", 
     442            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_Y",        source->chipY);                   // "y coord in warp input chip", 
     443            @>PS1_V4@               psFitsTableSetS16 (table, i, "PADDING3",          0);                               // "more padding",               
     444        }
     445
     446        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "MOMENTS_R1",        moments.Mrf);                     // "first radial moment",                       
     447        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "MOMENTS_RH",        moments.Mrh);                     // "half radial moment",                         
     448        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX",         moments.Krf);                     // "Kron Flux (in 2.5 R1)",                     
     449        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_ERR",     moments.dKrf);                    // "Kron Flux Error",                           
     450        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_INNER",   moments.Kinner);                  // "Kron Flux (in 1.0 R1)",                     
     451        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_OUTER",   moments.Kouter);                  // "Kron Flux (in 4.0 R1)",                     
     452        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_RAD",     source->skyRadius);               // "Radius where object hits sky",               
     453        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_FLUX",    source->skyFlux);                 // "Flux / pix where object hits sky",           
     454        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_SLOPE",   source->skySlope);                // "d(Flux/pix)/dRadius where object hits sky", 
     455
     456        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_NUM",      source->chipNum);                 // "id of warp input chip",       
     457        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_X",        source->chipX);                   // "x coord in warp input chip", 
     458        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_Y",        source->chipY);                   // "y coord in warp input chip", 
     459        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "PADDING3",          0);                               // "more padding",               
     460
     461        @PS1_DV?@                   psFitsTableSetS32 (table, i, "DIFF_NPOS",         diffStats.nGood);                 // "nPos (n pix > 3 sigma)",                     
     462        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_FRATIO",       diffStats.fRatio);                // "fPos / (fPos + fNeg)",                       
     463        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_BAD",   diffStats.nRatioBad);             // "nPos / (nPos + nNeg)",                       
     464        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_MASK",  diffStats.nRatioMask);            // "nPos / (nPos + nMask)",                     
     465        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_ALL",   diffStats.nRatioAll);             // "nPos / (nGood + nMask + nBad)",             
     466
     467        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_R_P",          diffStats.Rp);                    // "distance to positive match source", 
     468        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_SN_P",         diffStats.SNp);                   // "signal-to-noise of pos match src",   
     469        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_R_M",          diffStats.Rm);                    // "distance to negative match source", 
     470        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_SN_M",         diffStats.SNm);                   // "signal-to-noise of neg match src",   
     471
     472        @ALL@                      psFitsTableSetU32 (table, i, "FLAGS",              source->mode);                    // "psphot analysis flags",                     
     473        @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetU32 (table, i, "FLAGS2",             source->mode2);                   // "psphot analysis flags",                     
     474        @PS1_V3,PS1_SV2@           psFitsTableSetS32 (table, i, "PADDING2",           0);                               // "padding",                                   
     475
     476        @ALL@                      psFitsTableSetU16 (table, i, "N_FRAMES",           source->nFrames);                 // "Number of frames overlapping source center",
     477        @ALL@                      psFitsTableSetS16 (table, i, "PADDING",            0);                               // "padding",                                   
     478    }
     479
     480    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
     481    psMetadata *header = psMetadataCopy(NULL, tableHeader);
     482    pmSourceMasksHeader(header);
     483
     484    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     485    if (!psFitsWriteTableNew(fits, header, table, extname)) {
     486        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     487        psFree(table);
     488        psFree(header);
     489        return false;
     490    }
     491    psFree(table);
     492    psFree(header);
     493
     494    return true;
     495}
     496
     497bool pmSourcesWrite_CMF_@CMFMODE@_Old (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
     498{
     499    // fprintf (stderr, "writing with %s\n", __func__);
     500    PS_ASSERT_PTR_NON_NULL(fits, false);
     501    PS_ASSERT_PTR_NON_NULL(sources, false);
     502    PS_ASSERT_PTR_NON_NULL(extname, false);
     503
     504    psArray *table;
     505    psMetadata *row;
     506
     507    pmChip *chip = readout->parent->parent;
     508
     509    // if the sequence is defined, write these in seq order; otherwise
     510    // write them in S/N order:
     511    if (sources->n > 0) {
     512        pmSource *source = (pmSource *) sources->data[0];
     513        if (source->seq == -1) {
     514            // let's write these out in S/N order
     515            sources = psArraySort (sources, pmSourceSortByFlux);
     516        } else {
     517            sources = psArraySort (sources, pmSourceSortBySeq);
     518        }
     519    }
     520
     521    table = psArrayAllocEmpty (sources->n);
     522
     523    float magOffset;
     524    float zeroptErr;
     525    float fwhmMajor;
     526    float fwhmMinor;
     527    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
     528
     529    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     530    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
     531    for (int i = 0; i < sources->n; i++) {
     532        // this is the source associated with this image
     533        pmSource *thisSource = sources->data[i];
     534
     535        // this is the "real" version of this source
     536        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     537
     538        // If source->seq is -1, source was generated in this analysis.  If source->seq is
     539        // not -1, source was read from elsewhere: in the latter case, preserve the source
     540        // ID.  source.seq is used instead of source.id since the latter is a const
     541        // generated on Alloc, and would thus be wrong for read in sources.
     542        if (source->seq == -1) {
     543            source->seq = i;
     544        }
     545
     546        // set the 'best' values for various output fields:
     547        pmSourceOutputs outputs;
     548        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
     549
     550        pmSourceOutputsMoments moments;
     551        pmSourceOutputsSetMoments (&moments, source);
     552
     553        @PS1_DV?@ pmSourceDiffStats diffStats;
     554        @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
     555        @PS1_DV?@ if (source->diffStats) {
     556        @PS1_DV?@     diffStats = *source->diffStats;
     557        @PS1_DV?@ }
     558
     559        row = psMetadataAlloc ();
     560
     561        // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
     562        // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment"
     563        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
     564        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     565        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     566        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
     567        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     568
     569        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
     570        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
     571        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
     572
     573        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     574        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
     575        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
     576        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
     577
     578        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
     579        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
     580
     581        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
     582        @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
     583        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              source->apRadius);
     584        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
     585        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
     586        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX",          PS_DATA_S32, "aperture unmasked pixels",                   source->apNpixels);
     587
     588        @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     589
     590        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
     591        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
     592       
     593        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
     594        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     595        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
     596
     597        @>=PS1_V3,>PS1_SV1@       psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     598        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
     599        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
     600
     601        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
     602        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
     603        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
     604
     605        // PSF shape parameters:
     606        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     607        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     608        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    170609        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE",         PS_DATA_F32, "k term if defined",                          outputs.psfCore);
    171610
    172         // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
     611        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
    173612        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ",        PS_DATA_F32, "PSF FWHM (major axis)",                   outputs.psfMajorFWHM);
    174613        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN",        PS_DATA_F32, "PSF FWHM (minor axis)",                   outputs.psfMinorFWHM);
    175614
    176615        // psf data quality
    177         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    178         @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    179         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
    180         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
    181 
    182         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
    183         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
    184         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
    185 
    186         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
    187         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
    188         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
    189         @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
    190 
    191         // Lensing parameters:
    192         if (source->lensingOBJ && source->lensingOBJ->smear) {
    193           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
    194           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
    195           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
    196           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
    197           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
    198         }
    199         if (source->lensingOBJ && source->lensingOBJ->shear) {
    200           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
    201           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
    202           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
    203           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1);
    204           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2);
    205         }
    206         if (false && source->lensingOBJ) {
    207           // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
    208           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e1);
    209           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e2);
    210         }
    211 
    212         if (source->lensingPSF && source->lensingPSF->smear) {
    213           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
    214           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
    215           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
    216           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1);
    217           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2);
    218         }
    219         if (source->lensingPSF && source->lensingPSF->shear) {
    220           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
    221           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
    222           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
    223           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1);
    224           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2);
    225         }
    226         if (source->lensingPSF) {
    227           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e1);
    228           @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e2);
    229         }
     616        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
     617        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
     618        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     619        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     620
     621        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     622        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     623        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     624
     625        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     626        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     627        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     628        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     629
     630        // Lensing parameters:
     631        if (source->lensingOBJ && source->lensingOBJ->smear) {
     632          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
     633          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
     634          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
     635          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
     636          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
     637        }
     638        if (source->lensingOBJ && source->lensingOBJ->shear) {
     639          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
     640          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
     641          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
     642          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1);
     643          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2);
     644        }
     645        if (false && source->lensingOBJ) {
     646          // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
     647          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e1);
     648          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e2);
     649        }
     650
     651        if (source->lensingPSF && source->lensingPSF->smear) {
     652          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
     653          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
     654          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
     655          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1);
     656          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2);
     657        }
     658        if (source->lensingPSF && source->lensingPSF->shear) {
     659          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
     660          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
     661          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
     662          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1);
     663          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2);
     664        }
     665        if (source->lensingPSF) {
     666          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e1);
     667          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e2);
     668        }
    230669
    231670        // if lensing params exist also include the backmapped chipID and chip coordinates
    232         if (source->lensingPSF && source->lensingPSF->shear) {
     671        if (source->lensingPSF && source->lensingPSF->shear) {
    233672            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",   PS_DATA_S16, "id of warp input chip",     source->chipNum);
    234673            @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",    PS_DATA_S16, "x coord in warp input chip",     source->chipX);
    … …  
    244683        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                      moments.Kouter);
    245684
    246         @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
    247         @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
    248         @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
    249 
    250         @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",     PS_DATA_S16, "id of warp input chip",                      source->chipNum);
    251         @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",       PS_DATA_S16, "x coord in warp input chip",                 source->chipX);
    252         @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y",       PS_DATA_S16, "y coord in warp input chip",                 source->chipY);
    253         @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "PADDING3",         PS_DATA_S16, "more padding", 0);
    254 
    255         @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
    256         @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
    257         @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
    258         @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
    259         @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
    260 
    261         @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
    262         @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",           diffStats.SNp);
    263         @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",          diffStats.Rm);
    264         @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",           diffStats.SNm);
    265 
    266         @ALL@                      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
    267         @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
    268         @PS1_V3,PS1_SV2@           psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
    269         @PS1_SV?@
    270 
    271           // note that this definition is inconsistent with the definition in
    272           // Ohana/src/libautocode/def.  This version creates a table with data not
    273           // properly aligned with the 8-byte boundaries. The structure defined by
    274           // libautocode does this, but has a different order of elements (and adds
    275           // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
    276           // I'll leave it. But in future a PS1_SV2 should be forced to match
    277           // libautocode. Note that addstar knows to detect the alternate version of
    278           // PS1_SV1 and correctly interpret its fields.
    279 
    280         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    281         @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
     685        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
     686        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
     687        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
     688
     689        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",     PS_DATA_S16, "id of warp input chip",                      source->chipNum);
     690        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",       PS_DATA_S16, "x coord in warp input chip",                 source->chipX);
     691        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y",       PS_DATA_S16, "y coord in warp input chip",                 source->chipY);
     692        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "PADDING3",         PS_DATA_S16, "more padding", 0);
     693
     694        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
     695        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
     696        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
     697        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
     698        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
     699
     700        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
     701        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",           diffStats.SNp);
     702        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",          diffStats.Rm);
     703        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",           diffStats.SNm);
     704
     705        @ALL@                      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     706        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
     707        @PS1_V3,PS1_SV2@           psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
     708        @PS1_SV?@
     709
     710          // note that this definition is inconsistent with the definition in
     711          // Ohana/src/libautocode/def.  This version creates a table with data not
     712          // properly aligned with the 8-byte boundaries. The structure defined by
     713          // libautocode does this, but has a different order of elements (and adds
     714          // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
     715          // I'll leave it. But in future a PS1_SV2 should be forced to match
     716          // libautocode. Note that addstar knows to detect the alternate version of
     717          // PS1_SV1 and correctly interpret its fields.
     718
     719        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
     720        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    282721
    283722        psArrayAdd (table, 100, row);
    … …  
    320759}
    321760
     761bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
     762{
     763  // bool status = pmSourcesWrite_CMF_@CMFMODE@_Old (fits, readout, sources, imageHeader, tableHeader, extname, recipe);
     764     bool status = pmSourcesWrite_CMF_@CMFMODE@_New (fits, readout, sources, imageHeader, tableHeader, extname, recipe);
     765    return status;
     766}
     767
    322768// read in a readout from the fits file
    323 psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header)
     769psArray *pmSourcesRead_CMF_@CMFMODE@_New (psFits *fits, psMetadata *header)
    324770{
     771    // fprintf (stderr, "reading with %s\n", __func__);
     772
    325773    PS_ASSERT_PTR_NON_NULL(fits, false);
    326774    PS_ASSERT_PTR_NON_NULL(header, false);
    … …  
    355803    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
    356804
    357     // convert the table to the pmSource entriesa
     805    // reads full table into memory
     806    psFitsTable *table = psFitsReadTableNew (fits);
     807
     808    // convert the table to the pmSource entries
     809    for (int i = 0; i < numSources; i++) {
     810        pmSource *source = pmSourceAlloc ();
     811        pmModel *model = pmModelAlloc (modelType);
     812        source->modelPSF  = model;
     813        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
     814
     815        // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
     816        PAR = model->params->data.F32;
     817        dPAR = model->dparams->data.F32;
     818
     819        @ALL@                       source->seq                 = psFitsTableGetU32 (&status, table, i, "IPP_IDET");
     820        @ALL@                       PAR[PM_PAR_XPOS]            = psFitsTableGetF32 (&status, table, i, "X_PSF");
     821        @ALL@                       PAR[PM_PAR_YPOS]            = psFitsTableGetF32 (&status, table, i, "Y_PSF");
     822        @ALL@                       dPAR[PM_PAR_XPOS]           = psFitsTableGetF32 (&status, table, i, "X_PSF_SIG");
     823        @ALL@                       dPAR[PM_PAR_YPOS]           = psFitsTableGetF32 (&status, table, i, "Y_PSF_SIG");
     824        @ALL@                       axes.major                  = psFitsTableGetF32 (&status, table, i, "PSF_MAJOR");
     825        @ALL@                       axes.minor                  = psFitsTableGetF32 (&status, table, i, "PSF_MINOR");
     826        @ALL@                       axes.theta                  = psFitsTableGetF32 (&status, table, i, "PSF_THETA");
     827        @ALL@                       axes.theta                  = axes.theta * PS_RAD_DEG;
     828       
     829        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
     830        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7]           = psFitsTableGetF32 (&status, table, i, "PSF_CORE");
     831        @>PS1_V4,>PS1_SV2,>PS1_DV3@ }
     832
     833        @ALL@                       PAR[PM_PAR_SKY]             = psFitsTableGetF32 (&status, table, i, "SKY");
     834        @ALL@                       dPAR[PM_PAR_SKY]            = psFitsTableGetF32 (&status, table, i, "SKY_SIGMA");
     835        @ALL@                       source->sky                 = PAR[PM_PAR_SKY];
     836        @ALL@                       source->skyErr              = dPAR[PM_PAR_SKY];
     837
     838        // XXX use these to determine PAR[PM_PAR_I0]?
     839        @ALL@                       source->psfMag              = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG");
     840        @ALL@                       source->psfMagErr           = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG_SIG");
     841        @ALL@                       source->apMag               = psFitsTableGetF32 (&status, table, i, "AP_MAG");
     842        @>PS1_V2,PS1_SV?,>PS1_DV1@  source->apMagRaw            = psFitsTableGetF32 (&status, table, i, "AP_MAG_RAW");
     843        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux              = psFitsTableGetF32 (&status, table, i, "AP_FLUX");
     844        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr           = psFitsTableGetF32 (&status, table, i, "AP_FLUX_SIG");
     845
     846        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
     847        // XXX add these to PS1_SV1?
     848        @>PS1_V2,PS1_SV?,PS1_DV?@   source->psfFlux             = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX");
     849        @>PS1_V2,PS1_SV?,PS1_DV?@   source->psfFluxErr          = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX_SIG");
     850
     851        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
     852        @ALL@     PAR[PM_PAR_I0]                                = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
     853        @ALL@     dPAR[PM_PAR_I0]                               = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
     854
     855        pmPSF_AxesToModel (PAR, axes, model->class->useReff);
     856
     857        @ALL@     float peakMag                                 = psFitsTableGetF32 (&status, table, i, "PEAK_FLUX_AS_MAG");
     858        @ALL@     float peakFlux                                = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN;
     859
     860        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
     861        @ALL@     source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
     862        @ALL@     source->peak->rawFlux = peakFlux;
     863        @ALL@     source->peak->smoothFlux = peakFlux;
     864        @ALL@     source->peak->xf                              = PAR[PM_PAR_XPOS]; // more accurate position
     865        @ALL@     source->peak->yf                              = PAR[PM_PAR_YPOS]; // more accurate position
     866        @ALL@     source->peak->dx                              = dPAR[PM_PAR_XPOS];
     867        @ALL@     source->peak->dy                              = dPAR[PM_PAR_YPOS];
     868
     869        @ALL@     source->pixWeightNotBad                       = psFitsTableGetF32 (&status, table, i, "PSF_QF");
     870        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor     = psFitsTableGetF32 (&status, table, i, "PSF_QF_PERFECT");
     871        @ALL@     source->crNsigma                              = psFitsTableGetF32 (&status, table, i, "CR_NSIGMA");
     872        @ALL@     source->extNsigma                             = psFitsTableGetF32 (&status, table, i, "EXT_NSIGMA");
     873        @ALL@     source->apRadius                              = psFitsTableGetF32 (&status, table, i, "AP_MAG_RADIUS");
     874        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels           = psFitsTableGetS32 (&status, table, i, "AP_NPIX");
     875
     876        // note that some older versions used PSF_PROBABILITY: this was not well defined.
     877        @ALL@     model->chisq                                  = psFitsTableGetF32 (&status, table, i, "PSF_CHISQ");
     878        @ALL@     model->nDOF                                   = psFitsTableGetS32 (&status, table, i, "PSF_NDOF");
     879        @ALL@     model->nPix                                   = psFitsTableGetS32 (&status, table, i, "PSF_NPIX");
     880
     881        @ALL@     source->moments = pmMomentsAlloc ();
     882        @ALL@     source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
     883        @ALL@     source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
     884
     885        @ALL@     source->moments->Mxx                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_XX");
     886        @ALL@     source->moments->Mxy                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_XY");
     887        @ALL@     source->moments->Myy                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_YY");
     888
     889        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
     890        // we are storing enough information so the output will be consistent with the input
     891        @>PS1_V2,PS1_SV?@ source->moments->Mxxx         = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3C");
     892        @>PS1_V2,PS1_SV?@ source->moments->Mxxy         =  0.00;
     893        @>PS1_V2,PS1_SV?@ source->moments->Mxyy         =  0.00;
     894        @>PS1_V2,PS1_SV?@ source->moments->Myyy         = -1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3S");
     895        @>PS1_V2,PS1_SV?@ source->moments->Mxxxx        = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4C");
     896        @>PS1_V2,PS1_SV?@ source->moments->Mxxxy        =  0.00;
     897        @>PS1_V2,PS1_SV?@ source->moments->Mxxyy        =  0.00;
     898        @>PS1_V2,PS1_SV?@ source->moments->Mxyyy        = -0.25 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4S");
     899        @>PS1_V2,PS1_SV?@ source->moments->Myyyy        =  0.00;
     900
     901        // Lensing parameters (on read if PS1_V5+)
     902        if (haveLensOBJ) {
     903            source->lensingOBJ = pmSourceLensingAlloc ();
     904            source->lensingOBJ->smear = pmLensingParsAlloc();
     905            source->lensingOBJ->shear = pmLensingParsAlloc();
     906
     907            @>PS1_V4@ source->lensingOBJ->smear->X11            = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_OBJ");
     908            @>PS1_V4@ source->lensingOBJ->smear->X12            = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_OBJ");
     909            @>PS1_V4@ source->lensingOBJ->smear->X22            = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_OBJ");
     910            @>PS1_V4@ source->lensingOBJ->smear->e1             = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_OBJ");
     911            @>PS1_V4@ source->lensingOBJ->smear->e2             = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_OBJ");
     912            @>PS1_V4@ source->lensingOBJ->shear->X11            = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_OBJ");
     913            @>PS1_V4@ source->lensingOBJ->shear->X12            = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_OBJ");
     914            @>PS1_V4@ source->lensingOBJ->shear->X22            = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_OBJ");
     915            @>PS1_V4@ source->lensingOBJ->shear->e1             = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_OBJ");
     916            @>PS1_V4@ source->lensingOBJ->shear->e2             = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_OBJ");
     917        }
     918
     919        @>PS1_V4@ source->chipNum                               = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_NUM");
     920        @>PS1_V4@ source->chipX                                 = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_X");
     921        @>PS1_V4@ source->chipY                                 = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_Y");
     922
     923        if (haveLensPSF) {
     924            source->lensingPSF = pmSourceLensingAlloc ();
     925            source->lensingPSF->smear = pmLensingParsAlloc();
     926            source->lensingPSF->shear = pmLensingParsAlloc();
     927           
     928            @>PS1_V4@            source->lensingPSF->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_PSF");
     929            @>PS1_V4@            source->lensingPSF->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_PSF");
     930            @>PS1_V4@            source->lensingPSF->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_PSF");
     931            @>PS1_V4@            source->lensingPSF->smear->e1  = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_PSF");
     932            @>PS1_V4@            source->lensingPSF->smear->e2  = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_PSF");
     933            @>PS1_V4@            source->lensingPSF->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_PSF");
     934            @>PS1_V4@            source->lensingPSF->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_PSF");
     935            @>PS1_V4@            source->lensingPSF->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_PSF");
     936            @>PS1_V4@            source->lensingPSF->shear->e1  = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_PSF");
     937            @>PS1_V4@            source->lensingPSF->shear->e2  = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_PSF");
     938            @>PS1_V4@            source->lensingPSF->e1         = psFitsTableGetF32 (&status, table, i, "LENS_E1_PSF");
     939            @>PS1_V4@            source->lensingPSF->e2         = psFitsTableGetF32 (&status, table, i, "LENS_E2_PSF");
     940        }
     941
     942        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psFitsTableGetF32 (&status, table, i, "MOMENTS_R1");
     943        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrh         = psFitsTableGetF32 (&status, table, i, "MOMENTS_RH");
     944        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFlux    = psFitsTableGetF32 (&status, table, i, "KRON_FLUX");
     945        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFluxErr = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_ERR");
     946        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFinner  = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_INNER");
     947        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFouter  = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_OUTER");
     948
     949        @>PS1_V3@                  source->skyRadius            = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_RAD");
     950        @>PS1_V3@                  source->skyFlux              = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_FLUX");
     951        @>PS1_V3@                  source->skySlope             = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_SLOPE");
     952
     953        @PS1_DV?@                  int nPos                     = psFitsTableGetS32 (&status, table, i, "DIFF_NPOS");
     954        @PS1_DV?@  if (nPos) {
     955        @PS1_DV?@      source->diffStats                        = pmSourceDiffStatsAlloc();
     956        @PS1_DV?@      source->diffStats->nGood                 = nPos;
     957        @PS1_DV?@      source->diffStats->fRatio                = psFitsTableGetF32 (&status, table, i, "DIFF_FRATIO");
     958        @PS1_DV?@      source->diffStats->nRatioBad             = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_BAD");
     959        @PS1_DV?@      source->diffStats->nRatioMask            = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_MASK");
     960        @PS1_DV?@      source->diffStats->nRatioAll             = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_ALL");
     961       
     962        @>PS1_DV1@      source->diffStats->Rp                   = psFitsTableGetF32 (&status, table, i, "DIFF_R_P");
     963        @>PS1_DV1@      source->diffStats->SNp                  = psFitsTableGetF32 (&status, table, i, "DIFF_SN_P");
     964        @>PS1_DV1@      source->diffStats->Rm                   = psFitsTableGetF32 (&status, table, i, "DIFF_R_M");
     965        @>PS1_DV1@      source->diffStats->SNm                  = psFitsTableGetF32 (&status, table, i, "DIFF_SN_M");
     966        @PS1_DV?@  }
     967
     968        @ALL@                      source->mode                 = psFitsTableGetU32 (&status, table, i, "FLAGS");
     969        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->mode2                = psFitsTableGetU32 (&status, table, i, "FLAGS2");
     970        @ALL@                      source->nFrames              = psFitsTableGetU16 (&status, table, i, "N_FRAMES");
     971        assert (status);
     972
     973        sources->data[i] = source;
     974    }
     975    psFree (table);
     976    return sources;
     977}
     978
     979// read in a readout from the fits file
     980psArray *pmSourcesRead_CMF_@CMFMODE@_Old (psFits *fits, psMetadata *header)
     981{
     982    // fprintf (stderr, "reading with %s\n", __func__);
     983
     984    PS_ASSERT_PTR_NON_NULL(fits, false);
     985    PS_ASSERT_PTR_NON_NULL(header, false);
     986
     987    bool status;
     988    psF32 *PAR, *dPAR;
     989    psEllipseAxes axes;
     990
     991    // define PSF model type
     992    int defaultModelType = pmModelClassGetType ("PS_MODEL_GAUSS");
     993    int modelType = -1;
     994
     995    // if header does not define the model, default to a gaussian
     996    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
     997    if (PSF_NAME != NULL) {
     998        modelType = pmModelClassGetType (PSF_NAME);
     999    }
     1000    // work around bug in psphotFullForce
     1001    if (modelType < 0) {
     1002        modelType = defaultModelType;
     1003    }
     1004    // assert (modelType > -1);
     1005
     1006    // do we expect to find lensing parameters?
     1007    bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ");
     1008    bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF");
     1009
     1010    // We get the size of the table, and allocate the array of sources first because the table
     1011    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
     1012    // the table is read blocks the free.  In fact, it's better to read the table row by row.
     1013    long numSources = psFitsTableSize(fits); // Number of sources in table
     1014    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
     1015
     1016    // convert the table to the pmSource entries
    3581017    for (int i = 0; i < numSources; i++) {
    3591018        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
    … …  
    3821041        @ALL@     axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");
    3831042        @ALL@     axes.theta        = axes.theta * PS_RAD_DEG;
    384        
    385         @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
    386         @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
    387         @>PS1_V4,>PS1_SV2,>PS1_DV3@ }
     1043       
     1044        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
     1045        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
     1046        @>PS1_V4,>PS1_SV2,>PS1_DV3@ }
    3881047
    3891048        @ALL@     PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
    … …  
    3961055        @ALL@     source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
    3971056        @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
    398         @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
    399         @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX");
    400         @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
     1057        @>PS1_V2,PS1_SV?,>PS1_DV1@  source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
     1058        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux    = psMetadataLookupF32 (&status, row, "AP_FLUX");
     1059        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
    4011060
    4021061        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
    403         // XXX add these to PS1_SV1?
    404         @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
    405         @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
     1062        // XXX add these to PS1_SV1?
     1063        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
     1064        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
    4061065
    4071066        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
    … …  
    4241083
    4251084        @ALL@     source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    426         @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
     1085        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
    4271086        @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    4281087        @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    4291088        @ALL@     source->apRadius  = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS");
    430         @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
     1089        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
    4311090
    4321091        // note that some older versions used PSF_PROBABILITY: this was not well defined.
    … …  
    4431102        @ALL@     source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
    4441103
    445         // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
    446         // we are storing enough information so the output will be consistent with the input
     1104        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
     1105        // we are storing enough information so the output will be consistent with the input
    4471106        @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
    4481107        @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.0;
    … …  
    4561115        @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0;
    4571116
    458         // Lensing parameters (on read if PS1_V5+)
    459         if (haveLensOBJ) {
    460           source->lensingOBJ = pmSourceLensingAlloc ();
    461           source->lensingOBJ->smear = pmLensingParsAlloc();
    462           source->lensingOBJ->shear = pmLensingParsAlloc();
    463 
    464           @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
    465           @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
    466           @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
    467           @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
    468           @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
    469           @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
    470           @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
    471           @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
    472           @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
    473           @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
    474         }
    475 
    476         @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM");
    477         @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X");
    478         @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y");
    479 
    480         if (haveLensPSF) {
    481           source->lensingPSF = pmSourceLensingAlloc ();
    482           source->lensingPSF->smear = pmLensingParsAlloc();
    483           source->lensingPSF->shear = pmLensingParsAlloc();
    484 
    485           @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
    486           @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
    487           @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
    488           @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
    489           @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
    490           @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
    491           @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
    492           @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
    493           @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
    494           @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
    495           @>PS1_V4@ source->lensingPSF->e1         = psMetadataLookupF32 (&status, row, "LENS_E1_PSF");
    496           @>PS1_V4@ source->lensingPSF->e2         = psMetadataLookupF32 (&status, row, "LENS_E2_PSF");
    497         }
     1117        // Lensing parameters (on read if PS1_V5+)
     1118        if (haveLensOBJ) {
     1119          source->lensingOBJ = pmSourceLensingAlloc ();
     1120          source->lensingOBJ->smear = pmLensingParsAlloc();
     1121          source->lensingOBJ->shear = pmLensingParsAlloc();
     1122
     1123          @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
     1124          @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
     1125          @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
     1126          @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
     1127          @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
     1128          @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
     1129          @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
     1130          @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
     1131          @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
     1132          @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
     1133        }
     1134
     1135        @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM");
     1136        @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X");
     1137        @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y");
     1138
     1139        if (haveLensPSF) {
     1140          source->lensingPSF = pmSourceLensingAlloc ();
     1141          source->lensingPSF->smear = pmLensingParsAlloc();
     1142          source->lensingPSF->shear = pmLensingParsAlloc();
     1143
     1144          @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
     1145          @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
     1146          @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
     1147          @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
     1148          @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
     1149          @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
     1150          @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
     1151          @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
     1152          @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
     1153          @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
     1154          @>PS1_V4@ source->lensingPSF->e1         = psMetadataLookupF32 (&status, row, "LENS_E1_PSF");
     1155          @>PS1_V4@ source->lensingPSF->e2         = psMetadataLookupF32 (&status, row, "LENS_E2_PSF");
     1156        }
    4981157
    4991158        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
    … …  
    5081167        @>PS1_V3@ source->skySlope             = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE");
    5091168
    510         @PS1_DV?@  int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
    511         @PS1_DV?@  if (nPos) {
    512         @PS1_DV?@      source->diffStats = pmSourceDiffStatsAlloc();
    513         @PS1_DV?@      source->diffStats->nGood      = nPos;
    514         @PS1_DV?@      source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
    515         @PS1_DV?@      source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
    516         @PS1_DV?@      source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
    517         @PS1_DV?@      source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
    518        
    519         @>PS1_DV1@      source->diffStats->Rp         = psMetadataLookupF32 (&status, row, "DIFF_R_P");
    520         @>PS1_DV1@      source->diffStats->SNp        = psMetadataLookupF32 (&status, row, "DIFF_SN_P");
    521         @>PS1_DV1@      source->diffStats->Rm         = psMetadataLookupF32 (&status, row, "DIFF_R_M");
    522         @>PS1_DV1@      source->diffStats->SNm        = psMetadataLookupF32 (&status, row, "DIFF_SN_M");
    523         @PS1_DV?@  }
     1169        @PS1_DV?@  int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
     1170        @PS1_DV?@  if (nPos) {
     1171        @PS1_DV?@      source->diffStats = pmSourceDiffStatsAlloc();
     1172        @PS1_DV?@      source->diffStats->nGood      = nPos;
     1173        @PS1_DV?@      source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
     1174        @PS1_DV?@      source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
     1175        @PS1_DV?@      source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
     1176        @PS1_DV?@      source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
     1177       
     1178        @>PS1_DV1@      source->diffStats->Rp         = psMetadataLookupF32 (&status, row, "DIFF_R_P");
     1179        @>PS1_DV1@      source->diffStats->SNp        = psMetadataLookupF32 (&status, row, "DIFF_SN_P");
     1180        @>PS1_DV1@      source->diffStats->Rm         = psMetadataLookupF32 (&status, row, "DIFF_R_M");
     1181        @>PS1_DV1@      source->diffStats->SNm        = psMetadataLookupF32 (&status, row, "DIFF_SN_M");
     1182        @PS1_DV?@  }
    5241183
    5251184        @ALL@                      source->mode       = psMetadataLookupU32 (&status, row, "FLAGS");
    … …  
    5331192
    5341193    return sources;
     1194}
     1195
     1196psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) {
     1197  // psArray *array = pmSourcesRead_CMF_@CMFMODE@_Old (fits, header);
     1198      psArray *array = pmSourcesRead_CMF_@CMFMODE@_New (fits, header);
     1199    return array;
    5351200}
    5361201
    … …  
    5991264    // write the radial profile apertures to header
    6001265    for (int i = 0; i < radMax->n; i++) {
    601         sprintf (keyword1, "RMIN_%02d", i);
    602         sprintf (keyword2, "RMAX_%02d", i);
    603         psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
    604         psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]);
     1266        sprintf (keyword1, "RMIN_%02d", i);
     1267        sprintf (keyword2, "RMAX_%02d", i);
     1268        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     1269        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]);
    6051270    }
    6061271
    6071272    // we write out all sources, regardless of quality.  the source flags tell us the state
    6081273    for (int i = 0; i < sources->n; i++) {
    609         // this is the source associated with this image
     1274        // this is the source associated with this image
    6101275        pmSource *thisSource = sources->data[i];
    6111276
    612         // this is the "real" version of this source
    613         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     1277        // this is the "real" version of this source
     1278        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    6141279
    6151280        // skip sources without measurements
    … …  
    6391304        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
    6401305
    641         float AxialRatio = NAN;
    642         float AxialTheta = NAN;
    643         pmSourceExtendedPars *extpars = source->extpars;
    644         if (extpars) {
    645             AxialRatio = extpars->axes.minor / extpars->axes.major;
    646             AxialTheta = extpars->axes.theta;
    647         }
     1306        float AxialRatio = NAN;
     1307        float AxialTheta = NAN;
     1308        pmSourceExtendedPars *extpars = source->extpars;
     1309        if (extpars) {
     1310            AxialRatio = extpars->axes.minor / extpars->axes.major;
     1311            AxialTheta = extpars->axes.theta;
     1312        }
    6481313        psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
    6491314        psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
    … …  
    6511316        // Petrosian measurements
    6521317        // XXX insert header data: petrosian ref radius, flux ratio
    653         // XXX check flags to see if Pet was measured
     1318        // XXX check flags to see if Pet was measured
    6541319        if (doPetrosian) {
    655             pmSourceExtendedPars *extpars = source->extpars;
     1320            pmSourceExtendedPars *extpars = source->extpars;
    6561321            if (extpars) {
    657                 // XXX note that this mag is either calibrated or instrumental depending on existence of zero point
    658                 float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
    659                 // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired
    660                 float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
    661                 if (repairMagErrors) {
    662                   // I need to add the kron error in quadrature becasue pet_error ignores the object flux
    663                   float Krf  = source->moments ? source->moments->KronFlux : NAN;
    664                   float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    665                   if (isfinite (Krf) && isfinite (dKrf)) {
    666                     magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf));
    667                   } else {
    668                     magErr = 1.0 / magErr;
    669                   }
    670                 }
     1322                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point
     1323                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
     1324                // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired
     1325                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
     1326                if (repairMagErrors) {
     1327                  // I need to add the kron error in quadrature becasue pet_error ignores the object flux
     1328                  float Krf  = source->moments ? source->moments->KronFlux : NAN;
     1329                  float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
     1330                  if (isfinite (Krf) && isfinite (dKrf)) {
     1331                    magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf));
     1332                  } else {
     1333                    magErr = 1.0 / magErr;
     1334                  }
     1335                }
    6711336                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
    6721337                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
    … …  
    7051370        // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
    7061371        if (doAnnuli) {
    707             psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
    708             psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
    709             psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
    710             psVectorInit (radSB, NAN);
    711             psVectorInit (radFlux, NAN);
    712             psVectorInit (radFill, NAN);
    713             if (!source->extpars) goto empty_annuli;
    714             if (!source->extpars->radProfile) goto empty_annuli;
    715             if (!source->extpars->radProfile->binSB) goto empty_annuli;
    716             psAssert (source->extpars->radProfile->binSum, "programming error");
    717             psAssert (source->extpars->radProfile->binFill, "programming error");
    718             psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
    719             psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
    720             psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
    721 
    722             // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    723             for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
    724                 radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
    725                 radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
    726                 radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
    727             }
    728 
    729         empty_annuli:
    730             psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
    731             psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
    732             psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
    733             psFree (radSB);
    734             psFree (radFlux);
    735             psFree (radFill);
    736         }
    737         if (nRow < 0) {
    738             nRow = row->list->n;
    739         } else {
    740             psAssert (nRow == row->list->n, "inconsistent row lengths");
    741         }
    742         psArrayAdd (table, 100, row);
    743         psFree (row);
     1372            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
     1373            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
     1374            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
     1375            psVectorInit (radSB, NAN);
     1376            psVectorInit (radFlux, NAN);
     1377            psVectorInit (radFill, NAN);
     1378            if (!source->extpars) goto empty_annuli;
     1379            if (!source->extpars->radProfile) goto empty_annuli;
     1380            if (!source->extpars->radProfile->binSB) goto empty_annuli;
     1381            psAssert (source->extpars->radProfile->binSum, "programming error");
     1382            psAssert (source->extpars->radProfile->binFill, "programming error");
     1383            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
     1384            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
     1385            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
     1386
     1387            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     1388            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
     1389                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
     1390                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
     1391                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
     1392            }
     1393
     1394        empty_annuli:
     1395            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
     1396            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
     1397            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
     1398            psFree (radSB);
     1399            psFree (radFlux);
     1400            psFree (radFill);
     1401        }
     1402        if (nRow < 0) {
     1403            nRow = row->list->n;
     1404        } else {
     1405            psAssert (nRow == row->list->n, "inconsistent row lengths");
     1406        }
     1407        psArrayAdd (table, 100, row);
     1408        psFree (row);
    7441409    }
    7451410   
    7461411    if (table->n == 0) {
    747         if (!psFitsWriteBlank (fits, outhead, extname)) {
    748             psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
    749             psFree(outhead);
    750             psFree(table);
    751             return false;
    752         }
    753         psFree (outhead);
    754         psFree (table);
    755         return true;
     1412        if (!psFitsWriteBlank (fits, outhead, extname)) {
     1413            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     1414            psFree(outhead);
     1415            psFree(table);
     1416            return false;
     1417        }
     1418        psFree (outhead);
     1419        psFree (table);
     1420        return true;
    7561421    }
    7571422   
    7581423    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
    7591424    if (!psFitsWriteTable (fits, outhead, table, extname)) {
    760         psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
    761         psFree (outhead);
    762         psFree(table);
    763         return false;
     1425        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     1426        psFree (outhead);
     1427        psFree(table);
     1428        return false;
    7641429    }
    7651430    psFree (outhead);
    … …  
    8851550        extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
    8861551
    887         // magErr may have been saved in inverted form
     1552        // magErr may have been saved in inverted form
    8881553        float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
    8891554        float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
    8901555        if (isfinite(mag)) {
    891           extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
    892           if (isfinite(magErr)) {
    893             if (repairMagErrors) {
    894               // I need to add the kron error in quadrature becasue pet_error ignores the object flux
    895               float Krf  = source->moments ? source->moments->KronFlux : NAN;
    896               float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
    897               if (isfinite (Krf) && isfinite (dKrf)) {
    898                 magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf));
    899               } else {
    900                 magErr = 1.0 / magErr;
    901               }
    902               extpars->petrosianFluxErr = extpars->petrosianFlux * magErr;
    903             } else {
    904               extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
    905             }
    906           }
    907         }
     1556          extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
     1557          if (isfinite(magErr)) {
     1558            if (repairMagErrors) {
     1559              // I need to add the kron error in quadrature becasue pet_error ignores the object flux
     1560              float Krf  = source->moments ? source->moments->KronFlux : NAN;
     1561              float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
     1562              if (isfinite (Krf) && isfinite (dKrf)) {
     1563                magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf));
     1564              } else {
     1565                magErr = 1.0 / magErr;
     1566              }
     1567              extpars->petrosianFluxErr = extpars->petrosianFlux * magErr;
     1568            } else {
     1569              extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
     1570            }
     1571          }
     1572        }
    9081573
    9091574        extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
    … …  
    9731638    int nParamMax = 0;
    9741639    for (int i = 0; i < sources->n; i++) {
    975         // this is the source associated with this image
     1640        // this is the source associated with this image
    9761641        pmSource *thisSource = sources->data[i];
    9771642
    978         // this is the "real" version of this source
    979         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     1643        // this is the "real" version of this source
     1644        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    9801645
    9811646        if (source->modelFits == NULL) continue;
    … …  
    10041669        pmSource *thisSource = sources->data[i];
    10051670
    1006         // this is the "real" version of this source
    1007         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     1671        // this is the "real" version of this source
     1672        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    10081673
    10091674        // XXX if no model fits are saved, write out modelEXT?
    … …  
    10181683
    10191684            // pmSourceExtFitPars *extPars = source->extFitPars->data[j];
    1020             // assert (extPars);
    1021 
    1022             // skip models which were not actually fitted
    1023             // XXX
    1024             if (model->flags & badModel) continue;
     1685            // assert (extPars);
     1686
     1687            // skip models which were not actually fitted
     1688            // XXX
     1689            if (model->flags & badModel) continue;
    10251690
    10261691            PAR = model->params->data.F32;
    … …  
    10291694            yPos = PAR[PM_PAR_YPOS];
    10301695
    1031             // for the extended source models, we do not always fit the centroid in the non-linear fitting process
    1032             // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
    1033             // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
    1034             // TRAIL : X,Y are fitted
    1035             //
    1036            
    1037             // XXX this should be based on what happened, not on the model type
    1038             if (model->type == modelTypeTrail) {
    1039                 xErr = dPAR[PM_PAR_XPOS];
    1040                 yErr = dPAR[PM_PAR_YPOS];
    1041             } else {
    1042                 // this is definitely an underestimate since it does not
    1043                 // account for the extent of the source
    1044                 xErr = fwhmMajor * model->magErr / 2.35;
    1045                 yErr = fwhmMinor * model->magErr / 2.35;
    1046             }
    1047 
    1048             @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
    1049             @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0;
    1050             @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0;
    1051             @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
    1052             @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD;
    1053             @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD;
    1054             @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD;
    1055             @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0;
    1056 
    1057             float kronFlux = source->moments ? source->moments->KronFlux : NAN;
    1058             float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
     1696            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
     1697            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
     1698            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
     1699            // TRAIL : X,Y are fitted
     1700            //
     1701           
     1702            // XXX this should be based on what happened, not on the model type
     1703            if (model->type == modelTypeTrail) {
     1704                xErr = dPAR[PM_PAR_XPOS];
     1705                yErr = dPAR[PM_PAR_YPOS];
     1706            } else {
     1707                // this is definitely an underestimate since it does not
     1708                // account for the extent of the source
     1709                xErr = fwhmMajor * model->magErr / 2.35;
     1710                yErr = fwhmMinor * model->magErr / 2.35;
     1711            }
     1712
     1713            @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
     1714            @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0;
     1715            @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0;
     1716            @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
     1717            @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD;
     1718            @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD;
     1719            @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD;
     1720            @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0;
     1721
     1722            float kronFlux = source->moments ? source->moments->KronFlux : NAN;
     1723            float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
    10591724
    10601725            row = psMetadataAlloc ();
    10611726
    1062             // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
    1063             // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
     1727            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
     1728            // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
    10641729            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
    10651730            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
    … …  
    10701735            @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model ra coordinate",                    raPos);
    10711736            @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model dec coordinate",                   decPos);
    1072             @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
    1073             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
     1737            @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
     1738            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
    10741739
    10751740            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
    10761741            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
    10771742
    1078             @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
    1079             @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
    1080             @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
    1081             @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
    1082             @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
    1083 
    1084             // EAM : adding for PV2 outputs:
    1085             @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
     1743            @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
     1744            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
     1745            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
     1746            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
     1747            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
     1748
     1749            // EAM : adding for PV2 outputs:
     1750            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
    10861751
    10871752            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
    1088             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE",   0, "plate scale at source (arcsec/pixel)",       pltScale);
     1753            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE",   0, "plate scale at source (arcsec/pixel)",       pltScale);
    10891754
    10901755            // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX",       0, "second moment in x",                      extPars->Mxx);
    … …  
    11031768
    11041769            // XXX these should be major and minor, not 'x' and 'y'
    1105             if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
    1106                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
    1107                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]); // this is not fitted
    1108                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
    1109                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
    1110                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            NAN); // this is not fitted, so error is NAN
    1111                 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
    1112             } else {
    1113                 if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY])) {
    1114                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",     0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]);
    1115                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",     0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]);
    1116                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",         0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]);
    1117                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]);
    1118                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]);
    1119                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
    1120                 } else {
    1121                     psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
    1122                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
    1123                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
    1124                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
    1125                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_SXX]);
    1126                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SYY]);
    1127                     psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_SXY]);
    1128                 }
    1129             }
     1770            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
     1771                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
     1772                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]); // this is not fitted
     1773                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
     1774                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
     1775                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            NAN); // this is not fitted, so error is NAN
     1776                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
     1777            } else {
     1778                if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY])) {
     1779                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",     0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]);
     1780                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",     0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]);
     1781                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",         0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]);
     1782                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]);
     1783                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]);
     1784                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
     1785                } else {
     1786                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
     1787                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
     1788                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
     1789                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
     1790                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_SXX]);
     1791                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SYY]);
     1792                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_SXY]);
     1793                }
     1794            }
    11301795
    11311796            // write out the other generic parameters
    … …  
    11401805
    11411806                snprintf (name, 64, "EXT_PAR_%02d", k);
    1142                
     1807               
    11431808                if (k < model->params->n) {
    11441809                    psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
    … …  
    11481813            }
    11491814
    1150             // optionally, write out the covariance matrix values
    1151             // XXX do I need to pad this to match the biggest covar matrix?
    1152             if (false && model->covar) {
    1153                 for (int iy = 0; iy < model->covar->numCols; iy++) {
    1154                     for (int ix = iy; ix < model->covar->numCols; ix++) {
    1155                         snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
    1156                         psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
    1157 
    1158                     }
    1159                 }                  
    1160             }
     1815            // optionally, write out the covariance matrix values
     1816            // XXX do I need to pad this to match the biggest covar matrix?
     1817            if (false && model->covar) {
     1818                for (int iy = 0; iy < model->covar->numCols; iy++) {
     1819                    for (int ix = iy; ix < model->covar->numCols; ix++) {
     1820                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
     1821                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
     1822
     1823                    }
     1824                }                  
     1825            }
    11611826            psArrayAdd (table, 100, row);
    11621827            psFree (row);
    … …  
    12801945            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
    12811946            // XXX add an error:
    1282             // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
    1283         }
    1284 
    1285         // NOTE: we no longer write out the covariance matrix
    1286         if (false) {
    1287             // read the covariance matrix
    1288             int nparams = model->params->n;
    1289             psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
    1290             for (int y = 0; y < nparams; y++) {
    1291                 for (int x = 0; x < nparams; x++) {
    1292                     char name[64];
    1293                     snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
    1294                     covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
    1295                 }
    1296             }
    1297             model->covar = covar;
    1298         }
    1299 
    1300         // we are only saving the values stored in dPAR[SXX,etc]
     1947            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
     1948        }
     1949
     1950        // NOTE: we no longer write out the covariance matrix
     1951        if (false) {
     1952            // read the covariance matrix
     1953            int nparams = model->params->n;
     1954            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
     1955            for (int y = 0; y < nparams; y++) {
     1956                for (int x = 0; x < nparams; x++) {
     1957                    char name[64];
     1958                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
     1959                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
     1960                }
     1961            }
     1962            model->covar = covar;
     1963        }
     1964
     1965        // we are only saving the values stored in dPAR[SXX,etc]
    13011966        dPAR[PM_PAR_SXX] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ_ERR");
    13021967        dPAR[PM_PAR_SYY] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN_ERR");
    13031968        dPAR[PM_PAR_SXY] = psMetadataLookupF32(&status, row, "EXT_THETA_ERR");
    13041969
    1305         // other parameters that we need to read
     1970        // other parameters that we need to read
    13061971        PAR[PM_PAR_SKY] = psMetadataLookupF32(&status, row, "SKY_EXT");
    13071972
    … …  
    13412006    // perform full non-linear fits / extended source analysis?
    13422007    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
    1343         psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
    1344         return true;
     2008        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
     2009        return true;
    13452010    }
    13462011
    … …  
    13852050    for (int i = 0; i < sources->n; i++) {
    13862051
    1387         // this is the source associated with this image
     2052        // this is the source associated with this image
    13882053        pmSource *thisSource = sources->data[i];
    13892054
    1390         // this is the "real" version of this source
    1391         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     2055        // this is the "real" version of this source
     2056        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    13922057
    13932058        // skip sources without radial aper measurements (or insufficient)
    1394         if (source->radialAper == NULL) continue;
     2059        if (source->radialAper == NULL) continue;
    13952060
    13962061        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
    13972062
    1398         for (int entry = 0; entry < source->radialAper->n; entry++) {
    1399 
    1400             // choose the convolved EXT model, if available, otherwise the simple one
    1401             pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
    1402             assert (radialAper);
    1403 
    1404             if (pmSourcePositionUseMoments(source)) {
    1405                 xPos = source->moments->Mx;
    1406                 yPos = source->moments->My;
    1407             } else {
    1408                 xPos = source->peak->xf;
    1409                 yPos = source->peak->yf;
    1410             }
    1411 
    1412             row = psMetadataAlloc ();
    1413 
    1414             // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
    1415             // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
    1416             // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment"
    1417             psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
    1418             psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
    1419             psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
    1420             if (fwhmValues) {
    1421                 psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
    1422             } else {
    1423                 psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
    1424             }
    1425 
    1426             // XXX if we have raw radial apertures, write them out here
    1427             psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
    1428             psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
    1429             psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
    1430             psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
    1431             psVectorInit (radFlux,    NAN);
    1432             psVectorInit (radFluxErr, NAN);
    1433             psVectorInit (radFill,    NAN);
    1434             if (!radialAper->flux) goto write_annuli;
    1435             if (!radialAper->fill) goto write_annuli;
    1436             psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
    1437             psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
    1438 
    1439             // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
    1440             for (int j = 0; j < radialAper->flux->n; j++) {
    1441                 radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
    1442                 radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
    1443                 radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
    1444                 radFill->data.F32[j]      = radialAper->fill->data.F32[j];
    1445             }
    1446 
    1447         write_annuli:
    1448             psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
    1449             psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
    1450             psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
    1451             psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
    1452             psFree (radFlux);
    1453             psFree (radFluxErr);
    1454             psFree (radFluxStdev);
    1455             psFree (radFill);
    1456 
    1457             psArrayAdd (table, 100, row);
    1458             psFree (row);
    1459         }
     2063        for (int entry = 0; entry < source->radialAper->n; entry++) {
     2064
     2065            // choose the convolved EXT model, if available, otherwise the simple one
     2066            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
     2067            assert (radialAper);
     2068
     2069            if (pmSourcePositionUseMoments(source)) {
     2070                xPos = source->moments->Mx;
     2071                yPos = source->moments->My;
     2072            } else {
     2073                xPos = source->peak->xf;
     2074                yPos = source->peak->yf;
     2075            }
     2076
     2077            row = psMetadataAlloc ();
     2078
     2079            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     2080            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
     2081            // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment"
     2082            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     2083            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
     2084            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
     2085            if (fwhmValues) {
     2086                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
     2087            } else {
     2088                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
     2089            }
     2090
     2091            // XXX if we have raw radial apertures, write them out here
     2092            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     2093            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
     2094            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     2095            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
     2096            psVectorInit (radFlux,    NAN);
     2097            psVectorInit (radFluxErr, NAN);
     2098            psVectorInit (radFill,    NAN);
     2099            if (!radialAper->flux) goto write_annuli;
     2100            if (!radialAper->fill) goto write_annuli;
     2101            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
     2102            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
     2103
     2104            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     2105            for (int j = 0; j < radialAper->flux->n; j++) {
     2106                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
     2107                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
     2108                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
     2109                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
     2110            }
     2111
     2112        write_annuli:
     2113            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
     2114            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
     2115            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
     2116            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
     2117            psFree (radFlux);
     2118            psFree (radFluxErr);
     2119            psFree (radFluxStdev);
     2120            psFree (radFill);
     2121
     2122            psArrayAdd (table, 100, row);
     2123            psFree (row);
     2124        }
    14602125    }
    14612126
    … …  
    15812246    // perform full non-linear fits / extended source analysis?
    15822247    if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
    1583         psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
    1584         return true;
     2248        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
     2249        return true;
    15852250    }
    15862251
    … …  
    16092274        pmSource *thisSource = sources->data[i];
    16102275
    1611         // this is the "real" version of this source
    1612         pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    1613 
    1614         // if we did not fit the galaxy model, modelFits will be NULL
     2276        // this is the "real" version of this source
     2277        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     2278
     2279        // if we did not fit the galaxy model, modelFits will be NULL
    16152280        if (source->modelFits == NULL) continue;
    16162281
    1617         // if we did not fit the galaxy model, galaxyFits will also be NULL
     2282        // if we did not fit the galaxy model, galaxyFits will also be NULL
    16182283        if (source->galaxyFits == NULL) continue;
    16192284
Note: See TracChangeset for help on using the changeset viewer.