IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34289


Ignore:
Timestamp:
Aug 3, 2012, 5:28:50 PM (14 years ago)
Author:
Serge CHASTEL
Message:

Fitting trails for MOPS

Location:
tags/ipp-20120802/ppTranslate/src
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • tags/ipp-20120802/ppTranslate/src/ppMops.c

    r34281 r34289  
    6767    }
    6868
    69 
    7069    psArray *detections = ppMopsRead(args); // Detections from each input
    7170    if (!detections) {
    … …  
    9291    psLibFinalize();
    9392
    94     fprintf (stderr, "found %d leaks at %s\n",
    95         psMemCheckLeaks2 (0,
    96                 NULL, stdout, false, 500), "ppMops");
     93/*     fprintf (stderr, "found %d leaks at %s\n",  */
     94/*      psMemCheckLeaks2 (0, */
     95/*              NULL, stdout, false, 500), "ppMops"); */
    9796
    9897    return PS_EXIT_SUCCESS;
    9998}
    10099
    101 
    102 #if 0
    103     ps
    104 
    105 
    106     psArray *detections = NULL;         // Detections
    107     psMetadata *header = NULL;          // Header for detections
    108     {
    109         psFits *fits = psFitsOpen(data->detections, "r"); // FITS file
    110         header = psFitsReadHeader(NULL, fits);
    111         if (!header) {
    112             psErrorStackPrint(stderr, "Unable to read header");
    113             psFitsClose(fits);
    114             psFree(data);
    115             exit(PS_EXIT_DATA_ERROR);
    116         }
    117         if (!psFitsMoveExtName(fits, IN_EXTNAME)) {
    118             psErrorStackPrint(stderr, "Unable to move to extension %s", IN_EXTNAME);
    119             psFitsClose(fits);
    120             psFree(header);
    121             psFree(data);
    122             exit(PS_EXIT_DATA_ERROR);
    123         }
    124         detections = psFitsReadTable(fits);
    125         psFitsClose(fits);
    126         if (!detections) {
    127             psErrorStackPrint(stderr, "Unable to read detections");
    128             psFree(data);
    129             psFree(header);
    130             exit(PS_EXIT_DATA_ERROR);
    131         }
    132     }
    133 
    134     // Translate the columns
    135     int numIn = detections->n;        // Number of rows in input
    136     int numOut = 0;                   // Number of rows in output
    137     psArray *output = psArrayAllocEmpty(numIn); // Output table
    138     double plateScale = 0.0;                 // Average plate scale
    139     for (int i = 0; i < numIn; i++) {
    140         psMetadata *inRow = detections->data[i]; // Input row
    141 
    142         double ra = psMetadataLookupF64(NULL, inRow, "RA_PSF");
    143         double dec = psMetadataLookupF64(NULL, inRow, "DEC_PSF");
    144         double mag = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG") + data->zp;
    145         double magErr = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG_SIG");
    146         float ext = psMetadataLookupF32(NULL, inRow, "EXT_NSIGMA");
    147         double xErr = psMetadataLookupF64(NULL, inRow, "X_PSF_SIG");
    148         double yErr = psMetadataLookupF64(NULL, inRow, "Y_PSF_SIG");
    149         double scale = psMetadataLookupF64(NULL, inRow, "PLTSCALE");
    150         double angle = psMetadataLookupF64(NULL, inRow, "POSANGLE");
    151         psU32 flags = psMetadataLookupU32(NULL, inRow, "FLAGS");
    152 
    153         if (!isfinite(mag) || !isfinite(magErr) || (flags & SOURCE_MASK) ||
    154             ((flags & PM_SOURCE_MODE_DEFECT) && !(flags & PM_SOURCE_MODE_MOMENTS_FAILURE))
    155             ) {
    156             // DEFECT can be due to bad moments
    157             continue;
    158         }
    159 
    160         psMetadata *outRow = output->data[numOut] = psMetadataAlloc(); // Output row
    161 
    162         numOut++;
    163         plateScale += scale;
    164 
    165         // XXX Not at all sure I've got the angles around the right way here...
    166         double cosAngle = cos(angle), sinAngle = sin(angle);
    167         double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
    168         double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    169         double errScale = scale / 3600.0;
    170         double raErr = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    171         double decErr = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    172 
    173         psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", ra);
    174         psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", dec);
    175         psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", raErr);
    176         psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", decErr);
    177         psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", mag);
    178         psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", magErr);
    179         psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", flags);
    180 
    181         // The below need fixing
    182         psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", ext);
    183         psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", 0.0);
    184         psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", 0.0);
    185         psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", 0.0);
    186         psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", 0.0);
    187     }
    188     output->n = numOut;
    189     plateScale /= numOut;
    190     psFree(detections);
    191 
    192     // Translate the header
    193     psMetadata *outHeader = psMetadataAlloc(); // Output header
    194     ppMopsVersionHeader(outHeader);
    195     {
    196         const char *ra = psMetadataLookupStr(NULL, header, "FPA.RA");
    197         const char *dec = psMetadataLookupStr(NULL, header, "FPA.DEC");
    198         const char *filter = psMetadataLookupStr(NULL, header, "FPA.FILTER");
    199         float airmass = psMetadataLookupF32(NULL, header, "FPA.AIRMASS");
    200         float exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
    201         double angle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
    202         double alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
    203         double az = psMetadataLookupF64(NULL, header, "FPA.AZ");
    204         psS64 imageid = psMetadataLookupS64(NULL, header, "IMAGEID");
    205         double mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + exptime / 2.0 / 3600 / 24;
    206 
    207         float psf = plateScale * 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
    208                                         psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    209 
    210         psMetadataAddStr(outHeader, PS_LIST_TAIL, "RA", 0, "Right ascension of boresight", ra);
    211         psMetadataAddStr(outHeader, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", dec);
    212         psMetadataAddF32(outHeader, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", airmass);
    213         psMetadataAddF64(outHeader, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", mjd);
    214         psMetadataAddStr(outHeader, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter);
    215         psMetadataAddF64(outHeader, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", exptime);
    216         psMetadataAddF64(outHeader, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", angle);
    217         psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", alt);
    218         psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", az);
    219         psMetadataAddS64(outHeader, PS_LIST_TAIL, "DIFFIMID", 0, "Difference image identifier", imageid);
    220         psMetadataAddStr(outHeader, PS_LIST_TAIL, "FPA_ID", 0, "Exposure name", data->exp_name);
    221         psMetadataAddS64(outHeader, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", data->exp_id);
    222         psMetadataAddBool(outHeader, PS_LIST_TAIL, "POSITIVE", 0, "Positive subtraction?", data->direction);
    223         psMetadataAddStr(outHeader, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    224         psMetadataAddF32(outHeader, PS_LIST_TAIL, "STARPSF", 0, "Stellar PSF (arcsec)", psf);
    225 
    226         // These are completely fake
    227         psMetadataAddF32(outHeader, PS_LIST_TAIL, "LIMITMAG", 0, "Limiting magnitude (FAKE)", 99.0);
    228         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE1", 0, "Detection efficiency (FAKE)", 0.0);
    229         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE2", 0, "Detection efficiency (FAKE)", 0.0);
    230         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE3", 0, "Detection efficiency (FAKE)", 0.0);
    231         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE4", 0, "Detection efficiency (FAKE)", 0.0);
    232         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE5", 0, "Detection efficiency (FAKE)", 0.0);
    233         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE6", 0, "Detection efficiency (FAKE)", 0.0);
    234         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE7", 0, "Detection efficiency (FAKE)", 0.0);
    235         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE8", 0, "Detection efficiency (FAKE)", 0.0);
    236         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE9", 0, "Detection efficiency (FAKE)", 0.0);
    237         psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE10", 0, "Detection efficiency (FAKE)", 0.0);
    238     }
    239     psFree(header);
    240 
    241     // Write the new table
    242     {
    243         psFits *fits = psFitsOpen(data->output, "w"); // FITS file
    244         if (!fits) {
    245             psErrorStackPrint(stderr, "Unable to open %s for writing", data->output);
    246             psFree(outHeader);
    247             psFree(output);
    248             psFree(data);
    249             exit(PS_EXIT_SYS_ERROR);
    250         }
    251 
    252         if (numOut == 0) {
    253             // Write dummy table
    254             psMetadata *outRow = psMetadataAlloc(); // Dummy output row
    255             psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", NAN);
    256             psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", NAN);
    257             psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", NAN);
    258             psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", NAN);
    259             psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
    260             psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", NAN);
    261             psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
    262 
    263             // The below need fixing
    264             psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN);
    265             psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", NAN);
    266             psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", NAN);
    267             psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", NAN);
    268             psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", NAN);
    269             if (!psFitsWriteTableEmpty(fits, outHeader, outRow, OUT_EXTNAME)) {
    270                 psErrorStackPrint(stderr, "Unable to write table.");
    271                 psFree(outHeader);
    272                 psFree(output);
    273                 psFree(outRow);
    274                 psFree(data);
    275                 exit(PS_EXIT_SYS_ERROR);
    276             }
    277             psFree(outRow);
    278         } else if (!psFitsWriteTable(fits, outHeader, output, OUT_EXTNAME)) {
    279             psErrorStackPrint(stderr, "Unable to write table.");
    280             psFree(outHeader);
    281             psFree(output);
    282             psFree(data);
    283             exit(PS_EXIT_SYS_ERROR);
    284         }
    285         psFitsClose(fits);
    286     }
    287 
    288     psFree(outHeader);
    289     psFree(output);
    290     psFree(data);
    291 
    292 #endif
    293 

  • tags/ipp-20120802/ppTranslate/src/ppMops.h

    r34281 r34289  
    3030} ppMopsArguments;
    3131
    32 #if 0
    33 #warning "IS THERE ANYTHING TO BE MODIFIED HERE?"
    34 TTYPE19 = 'PSF_CHISQ'          / label for field  19
    35 TFORM19 = '1E      '           / data format of field: 4-byte REAL
    36 TTYPE20 = 'CR_NSIGMA'          / label for field  20
    37 TFORM20 = '1E      '           / data format of field: 4-byte REAL
    38 TTYPE21 = 'EXT_NSIGMA'         / label for field  21
    39 TFORM21 = '1E      '           / data format of field: 4-byte REAL
    40 TTYPE22 = 'PSF_MAJOR'          / label for field  22
    41 TFORM22 = '1E      '           / data format of field: 4-byte REAL
    42 TTYPE23 = 'PSF_MINOR'          / label for field  23
    43 TFORM23 = '1E      '           / data format of field: 4-byte REAL
    44 TTYPE24 = 'PSF_THETA'          / label for field  24
    45 TFORM24 = '1E      '           / data format of field: 4-byte REAL
    46 TTYPE25 = 'PSF_QF  '           / label for field  25
    47 TFORM25 = '1E      '           / data format of field: 4-byte REAL
    48 TTYPE26 = 'PSF_NDOF'           / label for field  26
    49 TFORM26 = '1J      '           / data format of field: 4-byte INTEGER
    50 TTYPE27 = 'PSF_NPIX'           / label for field  27
    51 TFORM27 = '1J      '           / data format of field: 4-byte INTEGER
    52 TTYPE28 = 'MOMENTS_XX'         / label for field  28
    53 TFORM28 = '1E      '           / data format of field: 4-byte REAL
    54 TTYPE29 = 'MOMENTS_XY'         / label for field  29
    55 TFORM29 = '1E      '           / data format of field: 4-byte REAL
    56 TTYPE30 = 'MOMENTS_YY'         / label for field  30
    57 TFORM30 = '1E      '           / data format of field: 4-byte REAL
    58 #endif
    59 
    6032/// Parse arguments
    6133ppMopsArguments *ppMopsArgumentsParse(int argc, char *argv[]);
    … …  
    7547  long numGood;                       // Number of "good" detections
    7648  psS64 diffSkyfileId;                // unique id for input skyfile
    77   psMetadata *table;                  // Columns from the input file
     49  psMetadata *table;                  // Columns from the input file (SkyChip.psf extension)
    7850  psVector *x, *y;                    // Image coordinates
    7951  psVector *ra, *dec;                 // Sky coordinates
    … …  
    10375/// @returns 1 if EXTTYPE of "SkyChip.psf" is PS1_DV1
    10476/// @returns 2 if EXTTYPE of "SkyChip.psf" is PS1_DV2
     77/// @returns 3 if EXTTYPE of "SkyChip.psf" is PS1_DV3
    10578/// @returns 0 otherwise
    10679int ppMopsGetSkyChipPsfVersion(const psFits* fits);

  • tags/ipp-20120802/ppTranslate/src/ppMopsArguments.c

    r34281 r34289  
    3434    fprintf(stderr, "\n");
    3535    psArgumentHelp(arguments);
    36     fprintf(stderr, "\t\tCMF file version can be set to either 1 for PS1_DV1 or 2 for PS1_DV1\n");
     36    fprintf(stderr, "\t\tCMF file version can be set to either 1 for PS1_DV1, 2 for PS1_DV2, or 3 for PS1_DV3\n");
    3737    fprintf(stderr, "\t\tSee IPP-MOPS ICD for details\n");
    3838    psLibFinalize();

  • tags/ipp-20120802/ppTranslate/src/ppMopsGetSkyChipPsfVersion.c

    r34281 r34289  
    1414    psFree(headerSkyChip);
    1515    return 2;
     16  } else if (strcmp(version, "PS1_DV3") == 0) {
     17    psFree(headerSkyChip);
     18    return 2;
    1619  }
    1720  psWarning("Unsupported EXTTYPE in SkyChip.psf table: [%s]", version);

  • tags/ipp-20120802/ppTranslate/src/ppMopsMerge.c

    r34281 r34289  
    7171
    7272    for (int i = 0; i < numInputs; i++) {
     73      psTrace("ppMops.merge", 1, "Processing batch %d\n", i);
    7374        ppMopsDetections *det = detections->data[i]; // Detections of interest
    7475        if (!det) {

  • tags/ipp-20120802/ppTranslate/src/ppMopsRead.c

    r34281 r34289  
    88#include "ppMops.h"
    99
     10static void addDummyValues(psMetadata* md, long size);
     11static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes);
     12
    1013/*
    1114  ppMopsRead possibly modifies the args->version if the user did not
    1215  set it explicitely.
    13  */
     16*/
    1417psArray *ppMopsRead(ppMopsArguments *args)
    1518{
    16     psTrace("ppMops.read", 1, "Reading input detections\n");
    17 
    18     psArray *inNames = args->input;          // Input names
    19     long num = inNames->n;                   // Number of inputs
    20     psArray *detections = psArrayAlloc(num); // Array of detections, to return
    21     for (int i = 0; i < num; i++) {
    22         const char *name = inNames->data[i];
    23 
    24         psFits *fits = psFitsOpen(name,  "r"); // FITS file
    25         if (!fits) {
    26             psError(PS_ERR_IO, false, "Unable to open input %d", i);
    27             return false;
    28         }
    29 
    30         psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
    31         if (!header) {
    32             psError(PS_ERR_IO, false, "Unable to read header %d", i);
    33             return false;
    34         }
    35 
    36         psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
    37         if (diffSkyfileId == 0) {
    38             psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
    39             return false;
    40         }
    41 
    42         if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
    43             psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
    44             return false;
    45         }
    46         int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
    47         if (args->version == 0) {
    48           psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
    49           args->version = (unsigned short) skyChipPsfVersion;
     19  psTrace("ppMops.read", 1, "Reading input detections\n");
     20
     21  psArray *inNames = args->input;          // Input names
     22  long num = inNames->n;                   // Number of inputs
     23  psArray *detections = psArrayAlloc(num); // Array of detections, to return
     24  for (int i = 0; i < num; i++) {
     25    const char *name = inNames->data[i];
     26
     27    psFits *fits = psFitsOpen(name,  "r"); // FITS file
     28    if (!fits) {
     29      psError(PS_ERR_IO, false, "Unable to open input %d", i);
     30      return false;
     31    }
     32
     33    psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
     34    if (!header) {
     35      psError(PS_ERR_IO, false, "Unable to read header %d", i);
     36      return false;
     37    }
     38
     39    psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
     40    if (diffSkyfileId == 0) {
     41      psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
     42      return false;
     43    }
     44
     45    if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
     46      psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
     47      return false;
     48    }
     49    int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
     50    if (args->version == 0) {
     51      psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
     52      args->version = (unsigned short) skyChipPsfVersion;
     53    }
     54    if (skyChipPsfVersion == 0) {
     55      // Try to read with the user specified version?
     56      skyChipPsfVersion = args->version;
     57    }
     58    /* Display a warning message if there are version
     59       inconsistencies between the file and the flag (note that
     60       those inconsistencies might be wanted) */
     61    if (skyChipPsfVersion != args->version) {
     62      if (skyChipPsfVersion > args->version) {
     63        psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
     64                  skyChipPsfVersion, args->version);
     65      } else { // Necessarily: skyChipPsfVersion > args->version
     66        psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
     67                  skyChipPsfVersion, args->version);
     68      }
     69    }
     70
     71    long size = psFitsTableSize(fits); // Size of table
     72    if (size <= 0) {
     73      psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
     74      psErrorClear();
     75      psWarning("Ignoring input %d", i);
     76      psFree(header);
     77      psFitsClose(fits);
     78      continue;
     79    }
     80    ppMopsDetections *det = ppMopsDetectionsAlloc(size);
     81    detections->data[i] = det;
     82    det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
     83    det->num = size;
     84    det->diffSkyfileId = diffSkyfileId;
     85
     86    psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
     87
     88    det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
     89    det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
     90    det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
     91    det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
     92    det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
     93    det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
     94    det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
     95    det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
     96    det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
     97
     98    det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
     99                                 psMetadataLookupF32(NULL, header, "FWHM_MIN"));
     100
     101    det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     102    det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     103
     104    psFree(header);
     105
     106    psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
     107    if (!table) {
     108      psError(PS_ERR_IO, false, "Unable to read table %d", i);
     109      return NULL;
     110    }
     111    det->table = table;
     112    psTrace("ppMops.read", 19, "PSF table:\n%s\n", psMetadataConfigFormat(det->table));
     113    psTrace("ppMops.read", 19, "End of PSF table\n");
     114
     115    //fittedTrail extension (Supported for releases 2 and above)
     116    if (skyChipPsfVersion >= 2) {
     117      //First: append all the new columns that we want to the existing table
     118      addDummyValues(det->table, size);
     119      if (!psFitsMoveExtName(fits, "SkyChip.xfit")) {
     120        psTrace("ppMops.read", 3, "No fitted trails extension");
     121      } else {
     122        psTrace("ppMops.read", 3, "Fitted trails extension found\n");
     123        int fittedTrailsSize = psFitsTableSize(fits);
     124        if (fittedTrailsSize <= 0) {
     125          psErrorStackPrint(stderr, "Unable to determine size of fitted trails extension table %d", i);
     126          psTrace("ppMops.read", 3, "No entry in fitted trails extension!!!!\n");
     127          psErrorClear();
     128        } else {
     129          psMetadata* fittedTrails = psFitsReadTableAllColumns(fits); // Table of interest
     130          if (!fittedTrails) {
     131            psError(PS_ERR_IO, false, "Unable to read fittedTrails table in file %d", i);
     132            return NULL;
     133          }
     134          //Iterate on the different names and types expected in the fittedTrails parameters
     135          psVector* idet = psMetadataLookupVector(NULL, fittedTrails, "IPP_IDET");
     136          replaceDummyValuesF32("X_EXT", fittedTrails, det->table, idet);
     137          replaceDummyValuesF32("Y_EXT", fittedTrails, det->table, idet);
     138          replaceDummyValuesF32("X_EXT_SIG", fittedTrails, det->table, idet);
     139          replaceDummyValuesF32("Y_EXT_SIG", fittedTrails, det->table, idet);
     140          replaceDummyValuesF32("EXT_INST_MAG",  fittedTrails, det->table, idet);
     141          replaceDummyValuesF32("EXT_INST_MAG_SIG",  fittedTrails, det->table, idet);
     142          replaceDummyValuesF32("NPARAMS",  fittedTrails, det->table, idet);
     143          replaceDummyValuesF32("EXT_WIDTH_MAJ",  fittedTrails, det->table, idet);
     144          replaceDummyValuesF32("EXT_WIDTH_MIN",  fittedTrails, det->table, idet);
     145          replaceDummyValuesF32("EXT_THETA",  fittedTrails, det->table, idet);
     146          replaceDummyValuesF32("EXT_WIDTH_MAJ_ERR",  fittedTrails, det->table, idet);
     147          replaceDummyValuesF32("EXT_WIDTH_MIN_ERR",  fittedTrails, det->table, idet);
     148          replaceDummyValuesF32("EXT_THETA_ERR",  fittedTrails, det->table, idet);
    50149        }
    51         if (skyChipPsfVersion == 0) {
    52           // Try to read with the user specified version?
    53           skyChipPsfVersion = args->version;
    54         }
    55         /* Display a warning message if there are version
    56            inconsistencies between the file and the flag (note that
    57            those inconsistencies might be wanted) */
    58         if (skyChipPsfVersion != args->version) {
    59           if (skyChipPsfVersion > args->version) {
    60             psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
    61                       skyChipPsfVersion, args->version);
    62           } else { // Necessarily: skyChipPsfVersion > args->version
    63             psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
    64                       skyChipPsfVersion, args->version);           
    65           }
    66         }
    67 
    68         long size = psFitsTableSize(fits); // Size of table
    69         if (size <= 0) {
    70             psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
    71             psErrorClear();
    72             psWarning("Ignoring input %d", i);
    73             psFree(header);
    74             psFitsClose(fits);
    75             continue;
    76         }
    77         ppMopsDetections *det = ppMopsDetectionsAlloc(size);
    78         detections->data[i] = det;
    79         det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
    80         det->num = size;
    81         det->diffSkyfileId = diffSkyfileId;
    82 
    83         psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);
    84 
    85         det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
    86         det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
    87         det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
    88         det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
    89         det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
    90         det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
    91         det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
    92         det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
    93         det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;
    94 
    95         det->seeing = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
    96                                      psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    97 
    98         det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
    99         det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
    100 
    101         psFree(header);
    102 
    103         psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
    104         if (!table) {
    105             psError(PS_ERR_IO, false, "Unable to read table %d", i);
    106             return NULL;
    107         }
    108         det->table = table;
    109         psFitsClose(fits);
    110         if (args->version == 0) {
    111           if (skyChipPsfVersion < 2) {
    112             // XXX: TODO: Do we need to add dummy vectors for the missing columns?
    113            }
    114         }
    115 
    116         psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
    117         psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
    118 
    119         det->raErr = psVectorAlloc(size, PS_TYPE_F64);
    120         det->decErr = psVectorAlloc(size, PS_TYPE_F64);
    121         det->mask = psVectorAlloc(size, PS_TYPE_U8);
    122 
    123         // convert ra and dec to radians for use in the purge duplicates function
    124         det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    125         det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
    126         det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
    127         det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
    128         if (!det->ra || !det->dec || !det->x || !det->y) {
    129             psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
    130             return NULL;
    131         }
    132 
    133         // Add our new vectors to the table so that duplicates and masked items may be purged
    134         psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
    135         psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
    136 
    137         psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
    138 
    139         psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
    140         psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
    141         psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
    142         psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
    143         psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
    144         psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
    145         psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
    146 
    147         double plateScale = 0.0;        // Plate scale
    148         long numGood = 0;               // Number of good rows
    149         for (long row = 0; row < size; row++) {
    150 
    151             psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
    152             if (flags & SOURCE_MASK) {
    153                 psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
    154                 det->mask->data.U8[row] = 0xFF;
    155                 continue;
    156             }
    157 
    158             // Calculate error in RA, Dec
    159 
    160             double xErr = xErrV->data.F32[row];
    161             double yErr = yErrV->data.F32[row];
    162             double scale = scaleV->data.F32[row];
    163             double angle = angleV->data.F32[row];
    164 
    165             if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
    166                 !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
    167                 !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
    168                 !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
    169                 psTrace("ppMops.read", 10,
    170                         "Discarding row %ld from input %d because of non-finite values: "
    171                         "%f %f %lf %lf %f %f %f %f %f %f",
    172                         row, i,
    173                         det->x->data.F32[row], det->y->data.F32[row],
    174                         det->ra->data.F64[row], det->dec->data.F64[row],
    175                         mag->data.F32[row], magErr->data.F32[row],
    176                         xErr, yErr, scale, angle);
    177                 det->mask->data.U8[row] = 0xFF;
    178                 continue;
    179             }
    180 
    181             // XXX Not at all sure I've got the angles around the right way here...
    182             double cosAngle = cos(angle), sinAngle = sin(angle);
    183             double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
    184             double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
    185             double errScale = scale / 3600.0;
    186             det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
    187             det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
    188 
    189             det->mask->data.U8[row] = 0;
    190             plateScale += scale;
    191             numGood++;
    192         }
    193         det->seeing *= ((float) plateScale) / ((float) numGood);
    194 
    195         // Are we using numGood for anything outside of this function?
    196         det->numGood = numGood;
    197 
    198         if (isfinite(args->zp) && numGood > 0) {
    199             psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
    200         }
    201 
    202         psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
    203     }
    204 
    205     psTrace("ppMops.read", 1, "Done reading input detections\n");
    206 
    207     return detections;
    208 }
     150      }
     151    }
     152
     153    psFitsClose(fits);
     154
     155    if (args->version == 0) {
     156      if (skyChipPsfVersion < 2) {
     157        // XXX: TODO: Do we need to add dummy vectors for the missing columns?
     158      }
     159    }
     160
     161    psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
     162    psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");
     163
     164    det->raErr = psVectorAlloc(size, PS_TYPE_F64);
     165    det->decErr = psVectorAlloc(size, PS_TYPE_F64);
     166    det->mask = psVectorAlloc(size, PS_TYPE_U8);
     167
     168    // convert ra and dec to radians for use in the purge duplicates function
     169    det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     170    det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
     171    det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
     172    det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
     173    if (!det->ra || !det->dec || !det->x || !det->y) {
     174      psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
     175      return NULL;
     176    }
     177
     178    // Add our new vectors to the table so that duplicates and masked items may be purged
     179    psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
     180    psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);
     181
     182    psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);
     183
     184    psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
     185    psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
     186    psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
     187    psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
     188    psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
     189    psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
     190    psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");
     191
     192    double plateScale = 0.0;        // Plate scale
     193    long numGood = 0;               // Number of good rows
     194    for (long row = 0; row < size; row++) {
     195
     196      psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
     197      if (flags & SOURCE_MASK) {
     198        psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
     199        det->mask->data.U8[row] = 0xFF;
     200        continue;
     201      }
     202
     203      // Calculate error in RA, Dec
     204      double xErr = xErrV->data.F32[row];
     205      double yErr = yErrV->data.F32[row];
     206      double scale = scaleV->data.F32[row];
     207      double angle = angleV->data.F32[row];
     208
     209      if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
     210          !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
     211          !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
     212          !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
     213        psTrace("ppMops.read", 10,
     214                "Discarding row %ld from input %d because of non-finite values: "
     215                "%f %f %lf %lf %f %f %f %f %f %f",
     216                row, i,
     217                det->x->data.F32[row], det->y->data.F32[row],
     218                det->ra->data.F64[row], det->dec->data.F64[row],
     219                mag->data.F32[row], magErr->data.F32[row],
     220                xErr, yErr, scale, angle);
     221        det->mask->data.U8[row] = 0xFF;
     222        continue;
     223      }
     224
     225      // XXX Not at all sure I've got the angles around the right way here...
     226      double cosAngle = cos(angle), sinAngle = sin(angle);
     227      double cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
     228      double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
     229      double errScale = scale / 3600.0;
     230      det->raErr->data.F64[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
     231      det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);
     232      det->mask->data.U8[row] = 0;
     233      plateScale += scale;
     234      numGood++;
     235    }
     236    det->seeing *= ((float) plateScale) / ((float) numGood);
     237
     238    // Are we using numGood for anything outside of this function?
     239    det->numGood = numGood;
     240
     241    if (isfinite(args->zp) && numGood > 0) {
     242      psBinaryOp(mag, mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
     243    }
     244
     245    psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
     246  }
     247
     248  psTrace("ppMops.read", 1, "Done reading input detections\n");
     249
     250  return detections;
     251}
     252
     253static psVector* createDummyF32(long size) {
     254  psVector* dummy = psVectorAlloc(size, PS_TYPE_F32);
     255  psVectorInit(dummy, NAN);
     256  return dummy;
     257}
     258
     259static void addDummyValues(psMetadata* md, long size) {
     260  psMetadataAdd(md, PS_LIST_TAIL, "X_EXT", PS_DATA_VECTOR, "EXT model x coordinate", createDummyF32(size));
     261  psMetadataAdd(md, PS_LIST_TAIL, "Y_EXT", PS_DATA_VECTOR, "EXT model y coordinate", createDummyF32(size));
     262  psMetadataAdd(md, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_VECTOR, "Sigma in EXT x coordinate", createDummyF32(size));
     263  psMetadataAdd(md, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_VECTOR, "Sigma in EXT y coordinate", createDummyF32(size));
     264  psMetadataAdd(md, PS_LIST_TAIL, "EXT_INST_MAG",  PS_DATA_VECTOR, "EXT fit instrumental magnitude", createDummyF32(size));
     265  psMetadataAdd(md, PS_LIST_TAIL, "EXT_INST_MAG_SIG",  PS_DATA_VECTOR, "Sigma of PSF instrumental magnitude", createDummyF32(size));
     266  psMetadataAdd(md, PS_LIST_TAIL, "NPARAMS",  PS_DATA_VECTOR, "Number of model parameters", createDummyF32(size));
     267  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MAJ",  PS_DATA_VECTOR, "EXT width (major axis), length for trail", createDummyF32(size));
     268  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN",  PS_DATA_VECTOR, "EXT width (minor axis), sigma for trail", createDummyF32(size));
     269  psMetadataAdd(md, PS_LIST_TAIL, "EXT_THETA",  PS_DATA_VECTOR, "EXT orientation angle", createDummyF32(size));
     270  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",  PS_DATA_VECTOR, "EXT width error (major axis)", createDummyF32(size));
     271  psMetadataAdd(md, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",  PS_DATA_VECTOR, "EXT width error (minor axis)", createDummyF32(size));
     272  psMetadataAdd(md, PS_LIST_TAIL, "EXT_THETA_ERR",  PS_DATA_VECTOR, "EXT orientation angle (error)", createDummyF32(size));
     273}
     274
     275static void replaceDummyValuesF32(const char* colName, psMetadata* source, psMetadata* target, psVector* indexes) {
     276  psVector* source_vector = psMetadataLookupVector(NULL, source, colName);
     277  psVector* target_vector = psMetadataLookupVector(NULL, target, colName);
     278  for (long index = 0; index<indexes->n; ++index) {
     279    target_vector->data.F32[indexes->data.S64[index]] = source_vector->data.F32[index];
     280  }
     281}

  • tags/ipp-20120802/ppTranslate/src/ppMopsWrite.c

    r34281 r34289  
    1414bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
    1515{
    16 
    17     psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    18     if (!fits) {
    19         psError(PS_ERR_IO, false, "Unable to open output file.");
    20         return false;
    21     }
    22 
    23     psMetadata *header = psMetadataAlloc(); // Header to write
    24     psString source = ppTranslateSource(), version = ppTranslateVersion();
    25     psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
    26     psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
    27     ppTranslateVersionHeader(header);
    28     psFree(source);
    29     psFree(version);
    30 
    31     psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
    32     psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
    33     psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
    34     psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
    35     psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
    36     psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
    37     psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
    38     psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
    39 
    40     // Get these header words from the first input non null input
    41     ppMopsDetections *det = NULL;
    42     for (int d = 0; d < detections->n && det == NULL; d++) {
    43         det = detections->data[d];
    44     }
    45     if (det != NULL) {
    46         psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
    47         psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
    48         psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
    49         psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
    50         psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
    51         psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
    52         psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
    53         psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
    54         psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
    55         psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
    56     } else {
    57         psWarning("no inputs with surviving detections. output header will be incomplete");
    58     }
    59     psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    60     psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
    61     psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
    62     psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
    63 
    64     //field in header that tells about the CMF version
    65     char cmfVersion[8];
    66     sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
    67     psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
    68 
    69     // Find the total number of detections
    70 
    71     long total = 0;
    72     for (long i=0; i<detections->n; i++) {
    73         ppMopsDetections *det = detections->data[i];
    74         if (!det) {
    75             continue;
    76         }
    77         total += det->num;
    78     }
    79 
    80     psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
    81 
    82     if (total == 0) {
    83         // Write dummy table
    84         psMetadata *row = psMetadataAlloc(); // Output row
    85         psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
    86         psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
    87         psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
    88         psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
    89         psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
    90         psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
    91         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
    92         psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
    93         psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
    94         psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
    95         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
    96         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
    97         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
    98         psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
    99         psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
    100         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
    101         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
    102         psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
    103         psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
    104         psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
    105 
    106         psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
    107         psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
    108         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
    109         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
    110         psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
    111 
    112         if (args->version == 2) {
    113               // Write data of version 2 (see ICD)
    114               psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
    115                              NAN);
    116               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
    117                              NAN);
    118               psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
    119                              NAN);
    120               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
    121                              NAN);
    122               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
    123                              NAN);
    124               psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
    125                              NAN);
    126               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
    127                              NAN);
    128               psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
    129                              NAN);
    130               psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
    131                              NAN);
    132               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
    133                              NAN);
    134               psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
    135                              NAN);
    136               psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
    137                              NAN);
    138               psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
    139                              NAN);
    140               psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
    141                              NAN);
    142               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
    143                              NAN);
    144               psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
    145                              NAN);
    146               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
    147                              NAN);
    148               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
    149                              NAN);
    150               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
    151                              NAN);
    152               psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
    153                              NAN);
    154               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
    155                              NAN);
    156               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
    157                              NAN);
    158               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
    159                              NAN);
    160               psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
    161                              NAN);
    162               psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
    163                              0);
    164               psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
    165                              0);
    166               psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
    167                              0);
    168         }
    169 
    170         if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
    171             psErrorStackPrint(stderr, "Unable to write empty table.");
    172             psFree(header);
    173             psFree(row);
    174             return false;
    175         }
    176         psFree(row);
    177     } else {
    178 
    179 #define addColumn(_outName, _inName, _convertTo32) \
    180         if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
    181             psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
    182             return false; \
    183         }
    184 
    185         // Allocate the output table
    186         psMetadata *table = psMetadataAlloc();
    187         addColumn("RA", "RA_PSF", 0);
    188         addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
    189         addColumn("DEC", "DEC_PSF", 0);
    190         addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
    191         addColumn("MAG", "PSF_INST_MAG", 0);
    192         addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
    193         addColumn("PSF_CHI2", "PSF_CHISQ", 0);
    194         addColumn("PSF_DOF", "PSF_NDOF", 1);
    195         addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
    196         addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
    197         addColumn("PSF_MAJOR", NULL, 0);
    198         addColumn("PSF_MINOR", NULL, 0);     
    199         addColumn("PSF_THETA", NULL, 0);   
    200         addColumn("PSF_QUALITY", "PSF_QF", 0);
    201         addColumn("PSF_NPIX", NULL, 1);       
    202         addColumn("MOMENTS_XX", NULL, 0);
    203         addColumn("MOMENTS_XY", NULL, 0);
    204         addColumn("MOMENTS_YY", NULL, 0);
    205         addColumn("N_POS", "DIFF_NPOS", 1);
    206         addColumn("F_POS", "DIFF_FRATIO", 0);
    207         addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
    208         addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
    209         addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
    210         addColumn("FLAGS", "FLAGS", 1);
    211         addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
    212         if (args->version == 2) {
    213             addColumn("IPP_IDET", NULL, 1);
    214             addColumn("PSF_INST_FLUX", NULL, 0);
    215             addColumn("PSF_INST_FLUX_SIG", NULL, 0);
    216             addColumn("AP_MAG", NULL, 0);
    217             addColumn("AP_MAG_RAW", NULL, 0);
    218             addColumn("AP_MAG_RADIUS", NULL, 0);
    219             addColumn("AP_FLUX", NULL, 0);
    220             addColumn("AP_FLUX_SIG", NULL, 0);
    221             addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
    222             addColumn("CAL_PSF_MAG", NULL, 0);
    223             addColumn("CAL_PSF_MAG_SIG", NULL, 0);
    224             addColumn("SKY", NULL, 0);
    225             addColumn("SKY_SIGMA", NULL, 0);
    226             addColumn("PSF_QF_PERFECT", NULL, 0);
    227             addColumn("MOMENTS_R1", NULL, 0);
    228             addColumn("MOMENTS_RH", NULL, 0);
    229             addColumn("KRON_FLUX", NULL, 0);
    230             addColumn("KRON_FLUX_ERR", NULL, 0);
    231             addColumn("KRON_FLUX_INNER", NULL, 0);
    232             addColumn("KRON_FLUX_OUTER", NULL, 0);
    233             addColumn("DIFF_R_P", NULL, 0);
    234             addColumn("DIFF_SN_P", NULL, 0);
    235             addColumn("DIFF_R_M", NULL, 0);
    236             addColumn("DIFF_SN_M", NULL, 0);
    237             addColumn("FLAGS2", NULL, 1);
    238             addColumn("IPP_IDET", NULL, 0);
    239             addColumn("N_FRAMES", NULL, 0);
    240             addColumn("PADDING", NULL, 0);
    241         }
    242         if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
    243             psError(psErrorCodeLast(), false, "Unable to write table");
    244             return false;
    245         }
    246         psFree(table);
    247     }
    248 
    249     psFree(header);
    250     psFitsClose(fits);
    251 
    252     psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
    253 
    254     return true;
     16  psFits *fits = psFitsOpen(args->output, "w"); // FITS file
     17  if (!fits) {
     18    psError(PS_ERR_IO, false, "Unable to open output file.");
     19    return false;
     20  }
     21
     22  psMetadata *header = psMetadataAlloc(); // Header to write
     23  psString source = ppTranslateSource(), version = ppTranslateVersion();
     24  psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
     25  psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
     26  ppTranslateVersionHeader(header);
     27  psFree(source);
     28  psFree(version);
     29
     30  psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
     31  psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
     32  psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
     33  psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
     34  psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
     35  psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
     36  psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
     37  psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
     38
     39  // Get these header words from the first input non null input
     40  ppMopsDetections *det = NULL;
     41  for (int d = 0; d < detections->n && det == NULL; d++) {
     42    det = detections->data[d];
     43  }
     44  if (det != NULL) {
     45    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
     46    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
     47    psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
     48    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
     49    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
     50    psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
     51    psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
     52    psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
     53    psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     54    psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
     55  } else {
     56    psWarning("no inputs with surviving detections. output header will be incomplete");
     57  }
     58  psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
     59  psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
     60  psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
     61  psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
     62
     63  //field in header that tells about the CMF version
     64  char cmfVersion[8];
     65  sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
     66  psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);
     67
     68  // Find the total number of detections
     69
     70  long total = 0;
     71  for (long i=0; i<detections->n; i++) {
     72    ppMopsDetections *det = detections->data[i];
     73    if (!det) {
     74      continue;
     75    }
     76    total += det->num;
     77  }
     78
     79  psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);
     80
     81  if (total == 0) {
     82    // Write dummy table
     83    psMetadata *row = psMetadataAlloc(); // Output row
     84    psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
     85    psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
     86    psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
     87    psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
     88    psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
     89    psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
     90    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
     91    psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
     92    psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
     93    psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
     94    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
     95    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
     96    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
     97    psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
     98    psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
     99    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
     100    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
     101    psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
     102    psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
     103    psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
     104
     105    psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
     106    psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
     107    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
     108    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
     109    psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);
     110
     111    if (args->version == 2) {
     112      // Write data of version 2 (see ICD)
     113      psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
     114                     NAN);
     115      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
     116                     NAN);
     117      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
     118                     NAN);
     119      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
     120                     NAN);
     121      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
     122                     NAN);
     123      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
     124                     NAN);
     125      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
     126                     NAN);
     127      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
     128                     NAN);
     129      psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
     130                     NAN);
     131      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
     132                     NAN);
     133      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
     134                     NAN);
     135      psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
     136                     NAN);
     137      psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
     138                     NAN);
     139      psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
     140                     NAN);
     141      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
     142                     NAN);
     143      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
     144                     NAN);
     145      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
     146                     NAN);
     147      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
     148                     NAN);
     149      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)",
     150                     NAN);
     151      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
     152                     NAN);
     153      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
     154                     NAN);
     155      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
     156                     NAN);
     157      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
     158                     NAN);
     159      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
     160                     NAN);
     161      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
     162                     0);
     163      psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center",
     164                     0);
     165      psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding",
     166                     0);
     167      if (args->version == 3) {
     168        // TODO: Write data of version 3 (see ICD)
     169      }
     170    }
     171
     172    if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
     173      psErrorStackPrint(stderr, "Unable to write empty table.");
     174      psFree(header);
     175      psFree(row);
     176      return false;
     177    }
     178    psFree(row);
     179  } else {
     180
     181#define addColumn(_outName, _inName, _convertTo32)                      \
     182    if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
     183      psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
     184      return false;                                                     \
     185    }
     186
     187    // Allocate the output table
     188    psMetadata *table = psMetadataAlloc();
     189    addColumn("RA", "RA_PSF", 0);
     190    addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     191    addColumn("DEC", "DEC_PSF", 0);
     192    addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
     193    addColumn("MAG", "PSF_INST_MAG", 0);
     194    addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
     195    addColumn("PSF_CHI2", "PSF_CHISQ", 0);
     196    addColumn("PSF_DOF", "PSF_NDOF", 1);
     197    addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
     198    addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
     199    addColumn("PSF_MAJOR", NULL, 0);
     200    addColumn("PSF_MINOR", NULL, 0);     
     201    addColumn("PSF_THETA", NULL, 0);   
     202    addColumn("PSF_QUALITY", "PSF_QF", 0);
     203    addColumn("PSF_NPIX", NULL, 1);       
     204    addColumn("MOMENTS_XX", NULL, 0);
     205    addColumn("MOMENTS_XY", NULL, 0);
     206    addColumn("MOMENTS_YY", NULL, 0);
     207    addColumn("N_POS", "DIFF_NPOS", 1);
     208    addColumn("F_POS", "DIFF_FRATIO", 0);
     209    addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
     210    addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
     211    addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
     212    addColumn("FLAGS", "FLAGS", 1);
     213    addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
     214    if (args->version >= 2) {
     215      addColumn("IPP_IDET", NULL, 1);
     216      addColumn("PSF_INST_FLUX", NULL, 0);
     217      addColumn("PSF_INST_FLUX_SIG", NULL, 0);
     218      addColumn("AP_MAG", NULL, 0);
     219      addColumn("AP_MAG_RAW", NULL, 0);
     220      addColumn("AP_MAG_RADIUS", NULL, 0);
     221      addColumn("AP_FLUX", NULL, 0);
     222      addColumn("AP_FLUX_SIG", NULL, 0);
     223      addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
     224      addColumn("CAL_PSF_MAG", NULL, 0);
     225      addColumn("CAL_PSF_MAG_SIG", NULL, 0);
     226      addColumn("SKY", NULL, 0);
     227      addColumn("SKY_SIGMA", NULL, 0);
     228      addColumn("PSF_QF_PERFECT", NULL, 0);
     229      addColumn("MOMENTS_R1", NULL, 0);
     230      addColumn("MOMENTS_RH", NULL, 0);
     231      addColumn("KRON_FLUX", NULL, 0);
     232      addColumn("KRON_FLUX_ERR", NULL, 0);
     233      addColumn("KRON_FLUX_INNER", NULL, 0);
     234      addColumn("KRON_FLUX_OUTER", NULL, 0);
     235      addColumn("DIFF_R_P", NULL, 0);
     236      addColumn("DIFF_SN_P", NULL, 0);
     237      addColumn("DIFF_R_M", NULL, 0);
     238      addColumn("DIFF_SN_M", NULL, 0);
     239      addColumn("FLAGS2", NULL, 1);
     240      addColumn("IPP_IDET", NULL, 0);
     241      addColumn("N_FRAMES", NULL, 0);
     242      addColumn("PADDING", NULL, 0);
     243      addColumn("X_EXT", NULL, 0);
     244      addColumn("Y_EXT", NULL, 0);
     245      addColumn("X_EXT_SIG", NULL, 0);
     246      addColumn("Y_EXT_SIG", NULL, 0);
     247      addColumn("EXT_INST_MAG", NULL, 0);
     248      addColumn("EXT_INST_MAG_SIG", NULL, 0);
     249      addColumn("NPARAMS", NULL, 0);
     250      addColumn("EXT_WIDTH_MAJ", NULL, 0);
     251      addColumn("EXT_WIDTH_MIN", NULL, 0);
     252      addColumn("EXT_THETA", NULL, 0);
     253      addColumn("EXT_WIDTH_MAJ_ERR", NULL, 0);
     254      addColumn("EXT_WIDTH_MIN_ERR", NULL, 0);
     255      addColumn("EXT_THETA_ERR", NULL, 0);
     256    }
     257    if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
     258      psError(psErrorCodeLast(), false, "Unable to write table");
     259      return false;
     260    }
     261    psFree(table);
     262  }
     263
     264  psFree(header);
     265  psFitsClose(fits);
     266
     267  psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);
     268
     269  return true;
    255270}
    256271
     272// extension parameter values:
     273// 0: SkyChip.psf
     274// 1: SkyChip.xfit
     275// Any other value is ignored
    257276static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
    258277{
    259     if (inColumnName == NULL) {
    260         inColumnName = outColumnName;
    261     }
    262 
    263     psVector *out = NULL;
    264     if (convertTo32) {
    265         // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
    266         // We want to write 32 bits to the output.
    267         int next = 0;
    268         for (long i=0; i<detections->n; i++) {
    269             ppMopsDetections *det = detections->data[i];
    270             if (!det || det->num == 0) {
    271                 // no detections survived for this input
    272                 continue;
    273             }
    274             psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
    275             if (!in) {
    276                 psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
    277                 return false;
    278             }
    279             if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
    280                 psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
    281                     inColumnName, in->type.type);
    282                 return false;
    283             }
    284             if (out == NULL) {
    285                 // First time through set up the output vector and the copy parameters
    286                 if (in->type.type == PS_TYPE_S64) {
    287                     out = psVectorAlloc(outputSize, PS_TYPE_S32);
    288                 } else {
    289                     out = psVectorAlloc(outputSize, PS_TYPE_U32);
    290                 }
    291             }
    292             for (long d=0; d < det->num; d++) {
    293                 if (in->type.type == PS_TYPE_S64) {
    294                     out->data.S32[next++] = in->data.S64[d];
    295                 } else {
    296                     out->data.U32[next++] = in->data.U64[d];
    297                 }
    298             }
    299         }
    300     } else {
    301         void *next = NULL;
    302         int elementSize = 0;    // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy
    303         for (long i=0; i<detections->n; i++) {
    304             ppMopsDetections *det = detections->data[i];
    305             if (!det || det->num == 0) {
    306                 // no detections survived for this input
    307                 continue;
    308             }
    309             psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
    310             if (!in) {
    311                 psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
    312                 return false;
    313             }
    314             if (out == NULL) {
    315                 // First time through set up the output vector and the copy parameters
    316                 out = psVectorAlloc(outputSize, in->type.type);
    317                 next = (void *) out->data.U8;
    318                 switch (in->type.type) {
    319                     case PS_TYPE_S8:
    320                     elementSize = sizeof(psS8);
    321                     break;
    322                 case PS_TYPE_U8:
    323                     elementSize = sizeof(psU8);
    324                     break;
    325                 case PS_TYPE_S16:
    326                     elementSize = sizeof(psS16);
    327                     break;
    328                 case PS_TYPE_U16:
    329                     elementSize = sizeof(psU16);
    330                     break;
    331                 case PS_TYPE_S32:
    332                     elementSize = sizeof(psS32);
    333                     break;
    334                 case PS_TYPE_U32:
    335                     elementSize = sizeof(psU32);
    336                     break;
    337                 case PS_TYPE_S64:
    338                     elementSize = sizeof(psS64);
    339                     break;
    340                 case PS_TYPE_U64:
    341                     elementSize = sizeof(psU64);
    342                     break;
    343                 case PS_TYPE_F32:
    344                     elementSize = sizeof(psF32);
    345                     break;
    346                 case PS_TYPE_F64:
    347                     elementSize = sizeof(psF64);
    348                     break;
    349                 default:
    350                     psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
    351                     return false;
    352 
    353                 }
    354             }
    355             // We are doing nasty things here so we can use memcpy.
    356             // It would be safer to do a proper loop over the elements.
    357             long toCopy = det->num * elementSize;
    358             memcpy(next, in->data.U8, toCopy);
    359             next += toCopy;
    360         }
    361     }
    362 
    363     // Finally add the new column to the output table
    364     psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
    365     psFree(out);    // drop reference
    366 
    367     return true;
     278  if (inColumnName == NULL) {
     279    inColumnName = outColumnName;
     280  }
     281  psVector *out = NULL;
     282  if (convertTo32) {
     283    // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
     284    // We want to write 32 bits to the output.
     285    int next = 0;
     286    for (long i=0; i<detections->n; i++) {
     287      ppMopsDetections *det = detections->data[i];
     288      if (!det || det->num == 0) {
     289        // no detections survived for this input
     290        continue;
     291      }
     292      psVector *in = NULL;
     293      in = psMetadataLookupVector(NULL, det->table, inColumnName);
     294      if (!in) {
     295        psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is true)", inColumnName);
     296        return false;
     297      }
     298      if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
     299        psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
     300                inColumnName, in->type.type);
     301        return false;
     302      }
     303      if (out == NULL) {
     304        // First time through set up the output vector and the copy parameters
     305        if (in->type.type == PS_TYPE_S64) {
     306          out = psVectorAlloc(outputSize, PS_TYPE_S32);
     307        } else {
     308          out = psVectorAlloc(outputSize, PS_TYPE_U32);
     309        }
     310      }
     311      for (long d=0; d < det->num; d++) {
     312        if (in->type.type == PS_TYPE_S64) {
     313          out->data.S32[next++] = in->data.S64[d];
     314        } else {
     315          out->data.U32[next++] = in->data.U64[d];
     316        }
     317      }
     318    }
     319  } else {
     320    void *next = NULL;
     321    int elementSize = 0;    // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy
     322    for (long i=0; i<detections->n; i++) {
     323      ppMopsDetections *det = detections->data[i];
     324      if (!det || det->num == 0) {
     325        // no detections survived for this input
     326        continue;
     327      }
     328      psVector *in = NULL;
     329      in = psMetadataLookupVector(NULL, det->table, inColumnName);
     330      if (!in) {
     331        psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s (convertTo32 is false)", inColumnName);
     332        out = psVectorAlloc(outputSize, PS_TYPE_F32);
     333        psVectorInit(out, NAN);
     334        psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     335        return false;
     336      }
     337      if (out == NULL) {
     338        // First time through set up the output vector and the copy parameters
     339        out = psVectorAlloc(outputSize, in->type.type);
     340        next = (void *) out->data.U8;
     341        switch (in->type.type) {
     342        case PS_TYPE_S8:
     343          elementSize = sizeof(psS8);
     344          break;
     345        case PS_TYPE_U8:
     346          elementSize = sizeof(psU8);
     347          break;
     348        case PS_TYPE_S16:
     349          elementSize = sizeof(psS16);
     350          break;
     351        case PS_TYPE_U16:
     352          elementSize = sizeof(psU16);
     353          break;
     354        case PS_TYPE_S32:
     355          elementSize = sizeof(psS32);
     356          break;
     357        case PS_TYPE_U32:
     358          elementSize = sizeof(psU32);
     359          break;
     360        case PS_TYPE_S64:
     361          elementSize = sizeof(psS64);
     362          break;
     363        case PS_TYPE_U64:
     364          elementSize = sizeof(psU64);
     365          break;
     366        case PS_TYPE_F32:
     367          elementSize = sizeof(psF32);
     368          break;
     369        case PS_TYPE_F64:
     370          elementSize = sizeof(psF64);
     371          break;
     372        default:
     373          psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
     374          return false;
     375        }
     376      }
     377      // We are doing nasty things here so we can use memcpy.
     378      // It would be safer to do a proper loop over the elements.
     379      long toCopy = det->num * elementSize;
     380      memcpy(next, in->data.U8, toCopy);
     381      next += toCopy;
     382    }
     383  }
     384  // Finally add the new column to the output table
     385  psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
     386  psFree(out);    // drop reference
     387  return true;
    368388}
     389
    369390static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
    370391{
    371     psVector *out = psVectorAlloc(total, PS_TYPE_S64);
    372     long next = 0;
    373     for (long i = 0; i<detections->n; i++) {
    374         ppMopsDetections *det = detections->data[i];
    375         if (!det) {
    376             continue;
    377         }
    378         psS64 diffSkyfileId = det->diffSkyfileId;
    379         for (long j = 0; j < det->num; j++) {
    380             out->data.S64[next++] = diffSkyfileId;
    381         }
    382     }
    383     psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
    384     psFree(out);
    385     return true;
     392  psVector *out = psVectorAlloc(total, PS_TYPE_S64);
     393  long next = 0;
     394  for (long i = 0; i<detections->n; i++) {
     395    ppMopsDetections *det = detections->data[i];
     396    if (!det) {
     397      continue;
     398    }
     399    psS64 diffSkyfileId = det->diffSkyfileId;
     400    for (long j = 0; j < det->num; j++) {
     401      out->data.S64[next++] = diffSkyfileId;
     402    }
     403  }
     404  psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
     405  psFree(out);
     406  return true;
    386407}
Note: See TracChangeset for help on using the changeset viewer.