IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34181


Ignore:
Timestamp:
Jul 19, 2012, 1:08:31 PM (14 years ago)
Author:
eugene
Message:

adding DV1,DV2,SV1 to pmSourceIO_CMF.c.in

Location:
branches/eam_branches/ipp-20120627/psModules/src/objects
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120627/psModules/src/objects/mksource.pl

    r32633 r34181  
    1616%cmfmodes = ("PS1_V1", 1,
    1717             "PS1_V2", 2,
     18             "PS1_DV1", 2.5,
    1819             "PS1_V3", 3,
    19              "PS1_V4", 4);
     20             "PS1_V4", 4,
     21             "PS1_DV2", 6,
     22             "PS1_SV1", 6,
     23    );
    2024
    2125print "1: $cmfmodes{1}\n";
    … …  
    3438# @>MODE@ : remove and keep if cmfmode > MODE
    3539# @<MODE@ : remove and keep if cmfmode > MODE
     40# XXX need to add features: split @foo,bar,baz@ by commas
     41# treat each chunk as a rule
     42#
    3643
    3744foreach $line (@list) {

  • branches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF.c.in

    r34144 r34181  
    106106        pmSourceOutputsSetMoments (&moments, source);
    107107
     108        @DV1,DV2@ pmSourceDiffStats diffStats;
     109        @DV1,DV2@ pmSourceDiffStatsInit(&diffStats);
     110        @DV1,DV2@ if (source->diffStats) {
     111        @DV1,DV2@     diffStats = *source->diffStats;
     112        @DV1,DV2@ }
     113
    108114        row = psMetadataAlloc ();
    109         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
    110         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
    111         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
    112         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
    113         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
     115        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
     116        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
     117        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
     118        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
     119        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
    114120
    115121        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
    116         @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
    117         @=PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
    118 
    119         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
    120         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
    121         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
    122         @ALL@    psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
    123 
    124         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
    125         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
    126 
    127         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
    128         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
    129         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
    130         @<PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    131 
    132         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
    133         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
     122        @=PS1_V1@                 psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
     123        @=PS1_V1@                 psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
     124
     125        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
     126        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
     127        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
     128        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
     129
     130        @>PS1_V2,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
     131        @>PS1_V2,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
     132
     133        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
     134        @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
     135        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              outputs.apRadius);
     136        @=PS1_DV2@                psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
     137        @=PS1_DV2@                psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
     138
     139        @<PS1_V3,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     140
     141        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
     142        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
    134143       
    135144        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
    136         @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
    137         @>PS1_V1@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
    138 
    139         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
    140         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
    141         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
    142 
    143         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
    144         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
    145         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
    146 
    147         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
    148         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
    149         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
    150         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
    151         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
    152         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
    153         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
    154 
    155         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
    156         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
    157         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
    158 
    159         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
    160         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
    161         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
    162         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
    163 
    164         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
    165         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
    166         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
    167         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
    168         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
    169         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
    170 
    171         @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
    172         @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
    173         @>PS1_V3@ psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
    174 
    175         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
    176         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
    177         @>PS1_V2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
     145        @>PS1_V1,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
     146        @>PS1_V1,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
     147
     148        @>PS1_V2@                 psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
     149        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
     150        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
     151
     152        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
     153        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
     154        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
     155
     156        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
     157        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
     158        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
     159        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
     160        @>PS1_V2,PS1_SV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
     161        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
     162        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
     163
     164        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
     165        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
     166        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
     167
     168        @>PS1_V2,PS1_SV1@                 psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
     169        @>PS1_V2,PS1_SV1@                 psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
     170        @>PS1_V2,PS1_SV1@                 psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
     171        @>PS1_V2,PS1_SV1@                 psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
     172
     173        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",                        moments.Mrf);
     174        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",                         moments.Mrh);
     175        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Krf);
     176        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",                            moments.dKrf);
     177        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kinner);
     178        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 2.5 R1)",                      moments.Kouter);
     179
     180        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
     181        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
     182        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
     183
     184        @PS1_DV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
     185        @PS1_DV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
     186        @PS1_DV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
     187        @PS1_DV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
     188        @PS1_DV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
     189
     190        @PS1_DV2@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
     191        @PS1_DV2@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",           diffStats.SNp);
     192        @PS1_DV2@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",          diffStats.Rm);
     193        @PS1_DV2@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",           diffStats.SNm);
     194
     195        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
     196        @>PS1_V2,PS1_SV1,PS1_DV2@         psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
     197        @>PS1_V2,PS1_SV1,PS1_DV1,PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
    178198
    179199        // XXX not sure how to get this : need to load Nimages with weight?
    180         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
    181         @ALL@     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
     200        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
     201        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
    182202
    183203        psArrayAdd (table, 100, row);
    … …  
    285305
    286306        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
    287         @>PS1_V2@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
    288         @>PS1_V2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
     307        // XXX add these to PS1_SV1?
     308        @>PS1_V2,PS1_DV1,PS1_DV2@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
     309        @>PS1_V2,PS1_DV1,PS1_DV2@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
    289310
    290311        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
    … …  
    307328
    308329        @ALL@     source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
    309         @>PS1_V2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
     330        @>PS1_V2,PS1_SV1,PS1_DV2@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
    310331        @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
    311332        @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
    … …  
    325346        @ALL@     source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
    326347
    327         @>PS1_V2@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
    328         @>PS1_V2@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
    329         @>PS1_V2@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
    330         @>PS1_V2@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
    331 
    332         @>PS1_V2@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
    333         @>PS1_V2@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
     348        @>PS1_V2,PS1_SV1@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
     349        @>PS1_V2,PS1_SV1@ source->moments->Mrh         = psMetadataLookupF32 (&status, row, "MOMENTS_RH");
     350        @>PS1_V2,PS1_SV1@ source->moments->KronFlux    = psMetadataLookupF32 (&status, row, "KRON_FLUX");
     351        @>PS1_V2,PS1_SV1@ source->moments->KronFluxErr = psMetadataLookupF32 (&status, row, "KRON_FLUX_ERR");
     352
     353        @>PS1_V2,PS1_SV1@ source->moments->KronFinner  = psMetadataLookupF32 (&status, row, "KRON_FLUX_INNER");
     354        @>PS1_V2,PS1_SV1@ source->moments->KronFouter  = psMetadataLookupF32 (&status, row, "KRON_FLUX_OUTER");
    334355
    335356        @>PS1_V3@ source->skyRadius            = psMetadataLookupF32 (&status, row, "SKY_LIMIT_RAD");
    … …  
    339360        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
    340361        // we are storing enough information so the output will be consistent with the input
    341         @>PS1_V2@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
    342         @>PS1_V2@ source->moments->Mxxy = 0.0;
    343         @>PS1_V2@ source->moments->Mxyy = 0.0;
    344         @>PS1_V2@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
    345 
    346         @>PS1_V2@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
    347         @>PS1_V2@ source->moments->Mxxxy = 0.0;
    348         @>PS1_V2@ source->moments->Mxxyy = 0.0;
    349         @>PS1_V2@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
    350         @>PS1_V2@ source->moments->Myyyy = 0.0;
     362        @>PS1_V2,PS1_SV1@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
     363        @>PS1_V2,PS1_SV1@ source->moments->Mxxy = 0.0;
     364        @>PS1_V2,PS1_SV1@ source->moments->Mxyy = 0.0;
     365        @>PS1_V2,PS1_SV1@ source->moments->Myyy = -1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3S");
     366
     367        @>PS1_V2,PS1_SV1@ source->moments->Mxxxx = +1.00 * psMetadataLookupF32 (&status, row, "MOMENTS_M4C");
     368        @>PS1_V2,PS1_SV1@ source->moments->Mxxxy = 0.0;
     369        @>PS1_V2,PS1_SV1@ source->moments->Mxxyy = 0.0;
     370        @>PS1_V2,PS1_SV1@ source->moments->Mxyyy = -0.25 * psMetadataLookupF32 (&status, row, "MOMENTS_M4S");
     371        @>PS1_V2,PS1_SV1@ source->moments->Myyyy = 0.0;
    351372
    352373        @ALL@     source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
    353         @>PS1_V2@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
     374        @>PS1_V2,PS1_SV1@ source->mode2 = psMetadataLookupU32 (&status, row, "FLAGS2");
    354375        assert (status);
     376
     377        @PS1_DV1,PS1_DV2@  int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
     378        @PS1_DV1,PS1_DV2@  if (nPos) {
     379        @PS1_DV1,PS1_DV2@      source->diffStats = pmSourceDiffStatsAlloc();
     380        @PS1_DV1,PS1_DV2@      source->diffStats->nGood      = nPos;
     381        @PS1_DV1,PS1_DV2@      source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
     382        @PS1_DV1,PS1_DV2@      source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
     383        @PS1_DV1,PS1_DV2@      source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
     384        @PS1_DV1,PS1_DV2@      source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
     385       
     386        @PS1_DV2@      source->diffStats->Rp         = psMetadataLookupF32 (&status, row, "DIFF_R_P");
     387        @PS1_DV2@      source->diffStats->SNp        = psMetadataLookupF32 (&status, row, "DIFF_SN_P");
     388        @PS1_DV2@      source->diffStats->Rm         = psMetadataLookupF32 (&status, row, "DIFF_R_M");
     389        @PS1_DV2@      source->diffStats->SNm        = psMetadataLookupF32 (&status, row, "DIFF_SN_M");
     390        @DV1,DV2@  }
    355391
    356392        sources->data[i] = source;
    … …  
    409445    // write the radial profile apertures to header
    410446    for (int i = 0; i < radMax->n; i++) {
    411       sprintf (keyword1, "RMIN_%02d", i);
    412       sprintf (keyword2, "RMAX_%02d", i);
    413       psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
    414       psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     447        sprintf (keyword1, "RMIN_%02d", i);
     448        sprintf (keyword2, "RMAX_%02d", i);
     449        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     450        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
    415451    }
    416452
    417453    // we write out all sources, regardless of quality.  the source flags tell us the state
    418454    for (int i = 0; i < sources->n; i++) {
    419         pmSource *source = sources->data[i];
     455        // this is the source associated with this image
     456        pmSource *thisSource = sources->data[i];
     457
     458        // this is the "real" version of this source
     459        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    420460
    421461        // skip sources without measurements
    … …  
    470510                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", extpars->petrosianR90);
    471511                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)", extpars->petrosianR90Err);
     512                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", extpars->petrosianFill);
    472513            } else {
    473514                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
    … …  
    479520                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90",     PS_DATA_F32, "Petrosian R90 (pix)", NAN);
    480521                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_90_ERR", PS_DATA_F32, "Petrosian R90 Error (pix)",NAN);
     522                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_FILL",          PS_DATA_F32, "Petrosian Fill Factor", NAN);
    481523            }
    482524        }
    … …  
    539581        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
    540582        psFree (outhead);
    541     psFree(table);
    542     return false;
     583        psFree(table);
     584        return false;
    543585    }
    544586    psFree (outhead);
    545587    psFree (table);
    546588   
     589    return true;
     590}
     591
     592bool pmSourcesRead_CMF_PS1_@CMFMODE@_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
     593{
     594    PS_ASSERT_PTR_NON_NULL(fits, false);
     595    PS_ASSERT_PTR_NON_NULL(sources, false);
     596
     597    bool status;
     598    long numSources = psFitsTableSize(fits); // Number of sources in table
     599    if (numSources == 0) {
     600        psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n");
     601        return false;
     602    }
     603
     604    // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point
     605    float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP");
     606    float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME");
     607    float magOffset = zeropt + 2.5*log10(exptime);
     608
     609    for (long i = 0; i < numSources; i++) {
     610        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     611        if (!row) {
     612            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     613            psFree(row);
     614            return false;
     615        }
     616        // Find the source with this sequence number.
     617        // XXX: I am assuming that sources is sorted in order of seq
     618        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     619        pmSource *source = NULL;
     620#ifndef ASSUME_SORTED
     621        long j = seq < sources->n ? seq : sources->n - 1;
     622        for (; j >= 0; j--) {
     623            source = sources->data[j];
     624            if (source->seq == seq) {
     625                break;
     626            }
     627        }
     628#else
     629        long j = sourceIndex[seq];
     630        psAssert(j >= 0 && j < sources->n, "invalid sourceIndex");
     631        source = sources->data[j];
     632#endif
     633        if (!source) {
     634            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     635            psFree(row);
     636            return false;
     637        }
     638
     639        if (!source->extpars) {
     640            source->extpars = pmSourceExtendedParsAlloc ();
     641        }
     642        pmSourceExtendedPars *extpars = source->extpars;
     643
     644        // Assume that X_EXT Y_EXT and sigmas match the psf src so skip
     645
     646        // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by
     647        // psf size or something?
     648        extpars->axes.major = 1.0;
     649        extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO");
     650        extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
     651
     652        float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
     653        float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
     654        if (isfinite(mag)) {
     655            extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
     656            if (isfinite(magErr)) {
     657                extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
     658            }
     659        }
     660
     661        extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
     662        extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR");
     663        extpars->petrosianR50      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50");
     664        extpars->petrosianR50Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR");
     665        extpars->petrosianR90      = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90");
     666        extpars->petrosianR90Err   = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR");
     667        extpars->petrosianFill     = psMetadataLookupF32(&status, row, "PETRO_FILL");
     668
     669        psVector *radSB   = psMetadataLookupVector(&status, row, "PROF_SB");
     670        psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX");
     671        psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL");
     672
     673        if (radSB && radSB->n > 0) {
     674            extpars->radProfile = pmSourceRadialProfileAlloc();
     675            extpars->radProfile->binSB   = psMemIncrRefCounter(radSB);
     676            extpars->radProfile->binSum   = psMemIncrRefCounter(radFlux);
     677            extpars->radProfile->binFill = psMemIncrRefCounter(radFill);
     678        }
     679
     680        psFree(row);
     681    }
     682
    547683    return true;
    548684}
    … …  
    572708    int nParamMax = 0;
    573709    for (int i = 0; i < sources->n; i++) {
    574         pmSource *source = sources->data[i];
     710        // this is the source associated with this image
     711        pmSource *thisSource = sources->data[i];
     712
     713        // this is the "real" version of this source
     714        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     715
    575716        if (source->modelFits == NULL) continue;
    576717        for (int j = 0; j < source->modelFits->n; j++) {
    … …  
    586727    for (int i = 0; i < sources->n; i++) {
    587728
    588         pmSource *source = sources->data[i];
     729        pmSource *thisSource = sources->data[i];
     730
     731        // this is the "real" version of this source
     732        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
    589733
    590734        // XXX if no model fits are saved, write out modelEXT?
    … …  
    642786            // XXX these should be major and minor, not 'x' and 'y'
    643787            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
    644               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
    645               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]);
    646               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
    647               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
    648               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SIGMA]);
    649               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
     788                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
     789                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]);
     790                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
     791                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
     792                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SIGMA]);
     793                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
    650794            } else {
    651               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
    652               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
    653               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
    654               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
    655               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SIGMA]);
    656               psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
     795                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
     796                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
     797                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
     798                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
     799                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SIGMA]);
     800                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
    657801            }
    658802
    … …  
    676820            }
    677821
    678             // XXX other parameters which may be set.
    679             // XXX flag / value to define the model
    680             // XXX write out the model type, fit status flags
    681 
     822            // optionally, write out the covariance matrix values
     823            // XXX do I need to pad this to match the biggest covar matrix?
     824            if (model->covar) {
     825                for (int iy = 0; iy < model->covar->numCols; iy++) {
     826                    for (int ix = iy; ix < model->covar->numCols; ix++) {
     827                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
     828                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
     829
     830                    }
     831                }                   
     832            }
    682833            psArrayAdd (table, 100, row);
    683834            psFree (row);
    … …  
    709860}
    710861
    711 bool pmSourcesWrite_CMF_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     862bool pmSourcesRead_CMF_PS1_@CMFMODE@_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
    712863{
     864    PS_ASSERT_PTR_NON_NULL(fits, false);
     865    PS_ASSERT_PTR_NON_NULL(sources, false);
     866
     867    bool status;
     868    long numSources = psFitsTableSize(fits); // Number of sources in table
     869    if (numSources == 0) {
     870        psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n");
     871        return false;
     872    }
     873
     874    for (long i = 0; i < numSources; i++) {
     875        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     876        if (!row) {
     877            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     878            psFree(row);
     879            return false;
     880        }
     881        // Find the source with this sequence number.
     882        // XXX: I am assuming that sources is sorted in order of seq.
     883        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     884        long j = seq < sources->n ? seq : sources->n - 1;
     885        pmSource *source = NULL;
     886        for (; j >= 0; j--) {
     887            source = sources->data[j];
     888            if (source->seq == seq) {
     889                break;
     890            }
     891        }
     892        if (!source) {
     893            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     894            psFree(row);
     895            return false;
     896        }
     897        if (!source->modelFits) {
     898            // XXX: where to find the number of models to expect?
     899            source->modelFits = psArrayAllocEmpty(5);
     900        }
     901        psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE");
     902        if (!modelName) {
     903            psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i);
     904            psFree(row);
     905            return false;
     906        }
     907        pmModelType modelType = pmModelClassGetType(modelName);
     908        if (modelType < 0) {
     909            psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName);
     910            psFree(row);
     911            return false;
     912        }
     913        pmModel *model = pmModelAlloc(modelType);
     914
     915        psF32 *PAR = model->params->data.F32;
     916        psF32 *dPAR = model->dparams->data.F32;
     917
     918        PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT");
     919        PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT");
     920        dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG");
     921        dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG");
     922
     923        model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG");
     924        model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG");
     925
     926        psEllipseAxes axes;
     927        axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ");
     928        axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN");
     929        axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA");
     930        if (!pmPSF_AxesToModel(PAR, axes, modelType)) {
     931            // Do we need to fail here or can this happen?
     932            psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model");
     933            psFree(model);
     934            psFree(row);
     935            return false;
     936        }
     937        // XXX: clean this up
     938        if (model->params->n > 7) {
     939            PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
     940        }
     941        // read the covariance matrix
     942        int nparams = model->params->n;
     943        psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
     944        for (int y = 0; y < nparams; y++) {
     945            for (int x = 0; x < nparams; x++) {
     946                char name[64];
     947                snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
     948                covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
     949            }
     950        }
     951        model->covar = covar;
     952
     953        psArrayAdd(source->modelFits, 1, model);
     954        psFree(model);
     955
     956        psFree(row);
     957    }
     958
    713959    return true;
    714960}
     961
     962// **** write out the radial flux values for the sources for a given matched-PSF image
     963// **** how do we distinguish the matched-PSF images from the non-matched version
     964bool pmSourcesWrite_CMF_PS1_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname, psMetadata *recipe)
     965{
     966    bool status = false;
     967    psArray *table;
     968    psMetadata *row;
     969    psF32 xPos, yPos;
     970    char keyword1[80], keyword2[80];
     971
     972    // perform full non-linear fits / extended source analysis?
     973    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
     974        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
     975        return true;
     976    }
     977
     978    // create a header to hold the output data
     979    psMetadata *outhead = psMetadataAlloc ();
     980
     981    // write the links to the image header
     982    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "radial flux table extension", extname);
     983
     984    // we use this just to define the output vectors (which must be present for all objects)
     985    psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     986    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     987    psAssert (radMax, "this must have been defined and tested earlier!");
     988    psAssert (radMax->n, "this must have been defined and tested earlier!");
     989    psAssert (radMin->n == radMax->n, "inconsistent annular bins");
     990
     991    // write the radial profile apertures to header
     992    for (int i = 0; i < radMax->n; i++) {
     993      sprintf (keyword1, "RMIN_%02d", i);
     994      sprintf (keyword2, "RMAX_%02d", i);
     995      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
     996      psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "min radius for SB profile", radMax->data.F32[i]);
     997    }
     998
     999    // the FWHM values are available if we measured a psf-matched convolved set
     1000    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
     1001
     1002    // let's write these out in S/N order
     1003    sources = psArraySort (sources, pmSourceSortByFlux);
     1004
     1005    table = psArrayAllocEmpty (sources->n);
     1006
     1007    // we write out all sources, regardless of quality.  the source flags tell us the state
     1008    for (int i = 0; i < sources->n; i++) {
     1009
     1010        // this is the source associated with this image
     1011        pmSource *thisSource = sources->data[i];
     1012
     1013        // this is the "real" version of this source
     1014        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
     1015
     1016        // skip sources without radial aper measurements (or insufficient)
     1017        if (source->radialAper == NULL) continue;
     1018
     1019        // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
     1020
     1021        for (int entry = 0; entry < source->radialAper->n; entry++) {
     1022
     1023            // choose the convolved EXT model, if available, otherwise the simple one
     1024            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
     1025            assert (radialAper);
     1026
     1027            if (pmSourcePositionUseMoments(source)) {
     1028                xPos = source->moments->Mx;
     1029                yPos = source->moments->My;
     1030            } else {
     1031                xPos = source->peak->xf;
     1032                yPos = source->peak->yf;
     1033            }
     1034
     1035            row = psMetadataAlloc ();
     1036
     1037            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
     1038            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
     1039            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
     1040            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
     1041            if (fwhmValues) {
     1042                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
     1043            } else {
     1044                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
     1045            }
     1046
     1047            // XXX if we have raw radial apertures, write them out here
     1048            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     1049            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
     1050            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
     1051            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
     1052            psVectorInit (radFlux,    NAN);
     1053            psVectorInit (radFluxErr, NAN);
     1054            psVectorInit (radFill,    NAN);
     1055            if (!radialAper->flux) goto write_annuli;
     1056            if (!radialAper->fill) goto write_annuli;
     1057            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
     1058            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
     1059
     1060            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
     1061            for (int j = 0; j < radialAper->flux->n; j++) {
     1062                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
     1063                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
     1064                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
     1065                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
     1066            }
     1067
     1068        write_annuli:
     1069            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX",       PS_DATA_VECTOR, "flux within annuli",       radFlux);
     1070            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_DATA_VECTOR, "flux error in annuli",     radFluxErr);
     1071            psMetadataAdd (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_DATA_VECTOR, "flux standard deviation",  radFluxStdev);
     1072            psMetadataAdd (row, PS_LIST_TAIL, "APER_FILL",       PS_DATA_VECTOR, "fill factor of annuli",    radFill);
     1073            psFree (radFlux);
     1074            psFree (radFluxErr);
     1075            psFree (radFluxStdev);
     1076            psFree (radFill);
     1077
     1078            psArrayAdd (table, 100, row);
     1079            psFree (row);
     1080        }
     1081    }
     1082
     1083    if (table->n == 0) {
     1084        if (!psFitsWriteBlank (fits, outhead, extname)) {
     1085            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
     1086            psFree(outhead);
     1087            psFree(table);
     1088            return false;
     1089        }
     1090        psFree (outhead);
     1091        psFree (table);
     1092        return true;
     1093    }
     1094
     1095    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     1096    if (!psFitsWriteTable (fits, outhead, table, extname)) {
     1097        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
     1098        psFree (outhead);
     1099        psFree(table);
     1100        return false;
     1101    }
     1102    psFree (outhead);
     1103    psFree (table);
     1104    return true;
     1105}
     1106
     1107bool pmSourcesRead_CMF_PS1_@CMFMODE@_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex)
     1108{
     1109    PS_ASSERT_PTR_NON_NULL(fits, false);
     1110    PS_ASSERT_PTR_NON_NULL(sources, false);
     1111
     1112    bool status;
     1113    long numSources = psFitsTableSize(fits); // Number of sources in table
     1114    if (numSources == 0) {
     1115        psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n");
     1116        return false;
     1117    }
     1118
     1119    long       seq_first = -1;
     1120    long       seq_last = -1;
     1121    psVector   *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
     1122    long       max_entries = -1;
     1123    long       num_entries = -1;
     1124
     1125    for (long i = 0; i < numSources; i++) {
     1126        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
     1127        if (!row) {
     1128            psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i);
     1129            psFree(row);
     1130            return false;
     1131        }
     1132        // Find the source with this sequence number.
     1133        // XXX: I am assuming that sources is sorted in order of seq.
     1134        long seq = psMetadataLookupU32 (&status, row, "IPP_IDET");
     1135        long j = seq < sources->n ? seq : sources->n - 1;
     1136        pmSource *source = NULL;
     1137        for (; j >= 0; j--) {
     1138            source = sources->data[j];
     1139            if (source->seq == seq) {
     1140                break;
     1141            }
     1142        }
     1143        if (!source) {
     1144            psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq);
     1145            psFree(row);
     1146            return false;
     1147        }
     1148        if (seq_first == -1) {
     1149            seq_first = seq;
     1150        }
     1151        if (seq == seq_first) {
     1152            psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM");
     1153            psVectorAppend(fwhmValues, value);
     1154        }
     1155        if (seq == seq_last) {
     1156            num_entries++;
     1157        } else {
     1158            num_entries = 1;
     1159            seq_last = seq;
     1160        }
     1161        if (num_entries > max_entries) {
     1162            max_entries = num_entries;
     1163        }
     1164
     1165        if (!source->radialAper) {
     1166            // XXX: where to find the number of models to expect?
     1167            source->radialAper = psArrayAllocEmpty(5);
     1168        }
     1169        pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc();
     1170
     1171        radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX"));
     1172        radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV"));
     1173        radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR"));
     1174        radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL"));
     1175
     1176        psArrayAdd(source->radialAper, 1, radialAper);
     1177
     1178        psFree(radialAper);
     1179        psFree(row);
     1180    }
     1181
     1182    // check for consistency between the length of fwhmValues and the maximum number of entries for each row
     1183    if (fwhmValues->n != max_entries) {
     1184        psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n",
     1185            fwhmValues->n, max_entries);
     1186        psAssert(0, "fixme");
     1187    }
     1188
     1189    if (!readout->analysis) {
     1190        readout->analysis = psMetadataAlloc();
     1191    }
     1192
     1193    psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
     1194    psFree(fwhmValues);
     1195
     1196    return true;
     1197}

  • branches/eam_branches/ipp-20120627/psModules/src/objects/pmSourceIO_CMF_PS1_DV2.c

    r33690 r34181  
    7373        pmSource *source = (pmSource *) sources->data[0];
    7474        if (source->seq == -1) {
    75           // let's write these out in S/N order
    76           sources = psArraySort (sources, pmSourceSortByFlux);
     75            // let's write these out in S/N order
     76            sources = psArraySort (sources, pmSourceSortByFlux);
    7777        } else {
    78           sources = psArraySort (sources, pmSourceSortBySeq);
     78            sources = psArraySort (sources, pmSourceSortBySeq);
    7979        }
    8080    }
    … …  
    9898        // generated on Alloc, and would thus be wrong for read in sources.
    9999        if (source->seq == -1) {
    100           source->seq = i;
     100            source->seq = i;
    101101        }
    102102
Note: See TracChangeset for help on using the changeset viewer.