IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 28, 2009, 1:36:37 PM (18 years ago)
Author:
beaumont
Message:

Added QA plots to ppSub. Migrated psastro QA plots to psModules/extras.

Location:
branches/cnb_branch_20090113/psastro/src
Files:
1 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branch_20090113/psastro/src/Makefile.am

    r20805 r21208  
    5353        psastroErrorCodes.c         \
    5454        psastroVersion.c            \
    55         psastroVisual.c             \
    5655        psastroDefineFiles.c        \
    5756        psastroAnalysis.c           \

  • branches/cnb_branch_20090113/psastro/src/psastro.h

    r20805 r21208  
    1414# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
    1515
     16#if 0
    1617// this structure represents a fit to the logN / logS curve for a set of stars
    1718// logN = offset + slope * logS
    … …  
    2526    int sPeak;                          // sum of stars to peak bin
    2627} pmLumFunc;
     28#endif
    2729
    2830bool              psastroDataSave (pmConfig *config);
    … …  
    7779psString          psastroVersionLong(void);
    7880
    79 // psastroVisual functions
    80 bool psastroSetVisual (bool mode);
    81 bool psastroVisualClose();
    82 bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
    83 bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
    84 bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe);
    85 bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit);
    86 bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld);
    87 bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    88 bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd);
    89 bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    90 bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale);
    91 bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe);
    92 
    9381// demo plots
    9482bool              psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);

  • branches/cnb_branch_20090113/psastro/src/psastroArguments.c

    r20805 r21208  
    8888    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    8989        psArgumentRemove (N, &argc, argv);
    90         psastroSetVisual (true);
    9190        pmAstromSetVisual (true);
    9291    }

  • branches/cnb_branch_20090113/psastro/src/psastroAstromGuess.c

    r20805 r21208  
    2929    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    3030    if (!recipe) {
    31         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
    32         return false;
     31        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
     32        return false;
    3333    }
    3434
    … …  
    3636    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
    3737    if (!status) {
    38         useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     38        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
    3939    }
    4040
    … …  
    4242    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    4343    if (!input) {
    44         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    45         return false;
     44        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     45        return false;
    4646    }
    4747
    … …  
    4949    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    5050    if (!status) {
    51         psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
    52         return false;
    53     } 
     51        psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
     52        return false;
     53    }
    5454
    5555    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    6767    bool bilevelAstrometry = false;
    6868    if (!useModel) {
    69         psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
     69        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    7070    }
    7171
    … …  
    7575        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    7676
    77         if (!useModel) {
    78             if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
    79         }
     77        if (!useModel) {
     78            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
     79        }
    8080
    8181        if (newFPA) {
    8282            newFPA = false;
    83             while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
    84             while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
     83            while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
     84            while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
    8585            RAminSky = fpa->toSky->R - M_PI;
    8686            RAmaxSky = fpa->toSky->R + M_PI;
    8787        }
    8888
    89         // report and save the current best guess for the chip 0,0 pixel coordinates
    90         {
    91             psPlane ptCH, ptFP, ptTP;
    92             psSphere ptSky;
    93 
    94             ptCH.x = 0;
    95             ptCH.y = 0;
    96             psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    97             psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    98             psDeproject (&ptSky, &ptTP, fpa->toSky);
    99             psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    100 
    101             psVectorAppend (cornerL, ptFP.x);
    102             psVectorAppend (cornerM, ptFP.y);
    103             psVectorAppend (cornerP, ptTP.x);
    104             psVectorAppend (cornerQ, ptTP.y);
    105             psVectorAppend (cornerR, ptSky.r);
    106             psVectorAppend (cornerD, ptSky.d);
    107         }
     89        // report and save the current best guess for the chip 0,0 pixel coordinates
     90        {
     91            psPlane ptCH, ptFP, ptTP;
     92            psSphere ptSky;
     93
     94            ptCH.x = 0;
     95            ptCH.y = 0;
     96            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     97            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     98            psDeproject (&ptSky, &ptTP, fpa->toSky);
     99            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     100
     101            psVectorAppend (cornerL, ptFP.x);
     102            psVectorAppend (cornerM, ptFP.y);
     103            psVectorAppend (cornerP, ptTP.x);
     104            psVectorAppend (cornerQ, ptTP.y);
     105            psVectorAppend (cornerR, ptSky.r);
     106            psVectorAppend (cornerD, ptSky.d);
     107        }
    108108
    109109        // apply the new WCS guess data to all of the data in the readouts
    … …  
    119119                if (rawstars == NULL) { continue; }
    120120
    121                 *nStars += rawstars->n;
     121                *nStars += rawstars->n;
    122122                for (int i = 0; i < rawstars->n; i++) {
    123123                    pmAstromObj *raw = rawstars->data[i];
    … …  
    138138                }
    139139
    140                 // dump or plot the resulting projected positions
    141                 if (psTraceGetLevel("psastro.dump") > 0) {
    142                     psastroDumpRawstars (rawstars, fpa, chip);
    143                 }
    144 
    145                 psastroVisualPlotRawStars(rawstars, fpa, chip, recipe);
    146 
    147                 if (psTraceGetLevel("psastro.plot") > 0) {
    148                     psastroPlotRawstars (rawstars, fpa, chip, recipe);
    149                 }
     140                // dump or plot the resulting projected positions
     141                if (psTraceGetLevel("psastro.dump") > 0) {
     142                    psastroDumpRawstars (rawstars, fpa, chip);
     143                }
     144
     145                pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe);
     146
     147                if (psTraceGetLevel("psastro.plot") > 0) {
     148                    psastroPlotRawstars (rawstars, fpa, chip, recipe);
     149                }
    150150            }
    151151        }
    … …  
    157157    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
    158158    if (*nStars == 0) {
    159         psLogMsg ("psastro", 2, "no sources available for astrometry\n");
    160         psFree (view);
    161         return true;
    162     }
    163 
    164     psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
    165               DEG_RAD*RAmin, DEG_RAD*DECmin,
    166               DEG_RAD*RAmax, DEG_RAD*DECmax);
     159        psLogMsg ("psastro", 2, "no sources available for astrometry\n");
     160        psFree (view);
     161        return true;
     162    }
     163
     164    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
     165              DEG_RAD*RAmin, DEG_RAD*DECmin,
     166              DEG_RAD*RAmax, DEG_RAD*DECmax);
    167167
    168168    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
    … …  
    203203    pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    204204    if (bilevelAstrometry) {
    205         if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    206             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    207             return false;
    208         }
     205        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
     206            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     207            return false;
     208        }
    209209    } else {
    210         if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
    211             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    212             return false;
    213         }
     210        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
     211            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     212            return false;
     213        }
    214214    }
    215215    return true;
    … …  
    225225    // load mosaic-level astrometry?
    226226    if (phu) {
    227         char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    228         if (ctype) {
    229             *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    230         }
     227        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     228        if (ctype) {
     229            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     230        }
    231231    }
    232232    if (*bilevelAstrometry) {
    233         pmAstromReadBilevelMosaic (fpa, phu->header);
    234     } 
     233        pmAstromReadBilevelMosaic (fpa, phu->header);
     234    }
    235235    psFree (view);
    236236    return true;
    … …  
    245245    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    246246    if (!input) {
    247         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    248         return false;
     247        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     248        return false;
    249249    }
    250250
    … …  
    277277        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    278278
    279         // XXX we are currently inconsistent with marking the good vs the bad data
    280         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    281         // probably wrong since it implies there is no data!
    282 
    283         // skip chips for which the astrometry failed (NASTRO == 0)
    284         if (!chip->cells->n) goto skip_chip;
    285         pmCell *cell = chip->cells->data[0];
    286         if (!cell) goto skip_chip;
    287 
    288         if (!cell->readouts->n) goto skip_chip;
    289         pmReadout *readout = cell->readouts->data[0];
    290         if (!readout) goto skip_chip;
    291 
    292         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    293         if (!updates) goto skip_chip;
    294        
    295         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    296         if (!nAstro) goto skip_chip;
    297 
    298         float astError = psMetadataLookupF32 (&status, updates, "CERROR");
    299         if (fabs(astError) < 1e-6) goto skip_chip;
    300 
    301         psPlane ptCH, ptFP, ptTP;
    302         psSphere ptSky;
    303 
    304         ptCH.x = 0;
    305         ptCH.y = 0;
    306         psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    307         psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    308         psDeproject (&ptSky, &ptTP, fpa->toSky);
    309         psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    310 
    311         // new corner locations based on the calibrated astrometry
    312         psVectorAppend (cornerLn, ptFP.x);
    313         psVectorAppend (cornerMn, ptFP.y);
    314         psVectorAppend (cornerPn, ptTP.x);
    315         psVectorAppend (cornerQn, ptTP.y);
    316         psVectorAppend (cornerRn, ptSky.r);
    317         psVectorAppend (cornerDn, ptSky.d);
    318         psVectorAppend (cornerMK, 0);
    319         continue;
     279        // XXX we are currently inconsistent with marking the good vs the bad data
     280        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     281        // probably wrong since it implies there is no data!
     282
     283        // skip chips for which the astrometry failed (NASTRO == 0)
     284        if (!chip->cells->n) goto skip_chip;
     285        pmCell *cell = chip->cells->data[0];
     286        if (!cell) goto skip_chip;
     287
     288        if (!cell->readouts->n) goto skip_chip;
     289        pmReadout *readout = cell->readouts->data[0];
     290        if (!readout) goto skip_chip;
     291
     292        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     293        if (!updates) goto skip_chip;
     294
     295        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     296        if (!nAstro) goto skip_chip;
     297
     298        float astError = psMetadataLookupF32 (&status, updates, "CERROR");
     299        if (fabs(astError) < 1e-6) goto skip_chip;
     300
     301        psPlane ptCH, ptFP, ptTP;
     302        psSphere ptSky;
     303
     304        ptCH.x = 0;
     305        ptCH.y = 0;
     306        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     307        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     308        psDeproject (&ptSky, &ptTP, fpa->toSky);
     309        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     310
     311        // new corner locations based on the calibrated astrometry
     312        psVectorAppend (cornerLn, ptFP.x);
     313        psVectorAppend (cornerMn, ptFP.y);
     314        psVectorAppend (cornerPn, ptTP.x);
     315        psVectorAppend (cornerQn, ptTP.y);
     316        psVectorAppend (cornerRn, ptSky.r);
     317        psVectorAppend (cornerDn, ptSky.d);
     318        psVectorAppend (cornerMK, 0);
     319        continue;
    320320
    321321    skip_chip:
    322         // new corner locations based on the calibrated astrometry
    323         psVectorAppend (cornerLn, 0.0);
    324         psVectorAppend (cornerMn, 0.0);
    325         psVectorAppend (cornerPn, 0.0);
    326         psVectorAppend (cornerQn, 0.0);
    327         psVectorAppend (cornerRn, 0.0);
    328         psVectorAppend (cornerDn, 0.0);
    329         psVectorAppend (cornerMK, 1);
     322        // new corner locations based on the calibrated astrometry
     323        psVectorAppend (cornerLn, 0.0);
     324        psVectorAppend (cornerMn, 0.0);
     325        psVectorAppend (cornerPn, 0.0);
     326        psVectorAppend (cornerQn, 0.0);
     327        psVectorAppend (cornerRn, 0.0);
     328        psVectorAppend (cornerDn, 0.0);
     329        psVectorAppend (cornerMK, 1);
    330330    }
    331331
    … …  
    336336
    337337    for (int i = 0; i < cornerRo->n; i++) {
    338        
    339         psPlane ptTP;
    340         psSphere ptSky;
    341 
    342         ptSky.r = cornerRo->data.F32[i];
    343         ptSky.d = cornerDo->data.F32[i];
    344 
    345         psProject (&ptTP, &ptSky, fpa->toSky);
    346         psVectorAppend (cornerPs, ptTP.x);
    347         psVectorAppend (cornerQs, ptTP.y);
     338
     339        psPlane ptTP;
     340        psSphere ptSky;
     341
     342        ptSky.r = cornerRo->data.F32[i];
     343        ptSky.d = cornerDo->data.F32[i];
     344
     345        psProject (&ptTP, &ptSky, fpa->toSky);
     346        psVectorAppend (cornerPs, ptTP.x);
     347        psVectorAppend (cornerQs, ptTP.y);
    348348    }
    349349
    … …  
    351351    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    352352    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
    353    
     353
    354354    // fit the valid chips, mask the invalid chips
    355355    psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs);
    356356    psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs);
    357    
     357
    358358    // apply the linear fit...
    359359    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
    … …  
    364364    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
    365365
    366     psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
     366    pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
    367367
    368368    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    … …  
    374374    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
    375375    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
    376    
     376
    377377    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
    378378    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
    … …  
    382382    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
    383383    if (!header) {
    384         header = psMetadataAlloc();
    385         psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
    386         psFree (header);  // drop this reference
     384        header = psMetadataAlloc();
     385        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
     386        psFree (header);  // drop this reference
    387387    }
    388388
    … …  
    395395
    396396    if (DEBUG) {
    397         FILE *f = fopen ("corners.dat", "w");
    398         for (int i = 0; i < cornerRo->n; i++) {
    399             fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
    400                      cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
    401                      cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
    402         }
    403         fclose (f);
     397        FILE *f = fopen ("corners.dat", "w");
     398        for (int i = 0; i < cornerRo->n; i++) {
     399            fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
     400                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
     401                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
     402        }
     403        fclose (f);
    404404    }
    405405
    … …  
    424424    psFree (map);
    425425    psFree (view);
    426    
     426
    427427
    428428    return true;

  • branches/cnb_branch_20090113/psastro/src/psastroCleanup.c

    r20805 r21208  
    44
    55    psFree (config);
     6    pmAstromVisualClose ();
    67
    78    psTimerStop ();
    … …  
    1112    pmConceptsDone ();
    1213    pmConfigDone ();
    13     psastroVisualClose ();
    14     pmAstromVisualClose ();
    1514    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    1615    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");

  • branches/cnb_branch_20090113/psastro/src/psastroFixChips.c

    r20805 r21208  
    1313    bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
    1414    if (!status) {
    15         fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     15        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
    1616    }
    1717    if (!fixChips) return true;
    … …  
    2525    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    2626    if (!input) {
    27         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    28         return false;
     27        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     28        return false;
    2929    }
    3030
    … …  
    4848    // files associated with the science image
    4949    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    50         psError (PS_ERR_IO, false, "Can't load the astrometry model file");
    51         return false;
     50        psError (PS_ERR_IO, false, "Can't load the astrometry model file");
     51        return false;
    5252    }
    5353
    … …  
    6464
    6565    if (DEBUG) {
    66         f = fopen ("corners.raw.dat", "w");
    67         chipName = NULL;
     66        f = fopen ("corners.raw.dat", "w");
     67        chipName = NULL;
    6868    }
    6969
    7070    pmChip *obsChip = NULL;
    7171    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    72         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    73 
    74         // XXX we are currently inconsistent with marking the good vs the bad data
    75         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    76         // probably wrong since it implies there is no data!
    77 
    78         // skip chips for which the astrometry failed (NASTRO == 0)
    79         if (!obsChip->cells->n) continue;
    80         pmCell *cell = obsChip->cells->data[0];
    81         if (!cell) continue;
    82 
    83         if (!cell->readouts->n) continue;
    84         pmReadout *readout = cell->readouts->data[0];
    85         if (!readout) continue;
    86 
    87         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    88         if (!updates) continue;
    89        
    90         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    91         if (!nAstro) continue;
    92 
    93         // set the chip astrometry using the astrom file
    94         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    95 
    96         psRegion *region = pmChipPixels (obsChip);
    97         psPlane ptCP, ptFP;
    98 
    99         ptCP.x = region->x0; ptCP.y = region->y0;
    100         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    101         xObs->data.F32[nPts] = ptFP.x;
    102         yObs->data.F32[nPts] = ptFP.y;
    103         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    104         xRef->data.F32[nPts] = ptFP.x;
    105         yRef->data.F32[nPts] = ptFP.y;
    106 
    107         if (DEBUG) {
    108             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    109             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    110         }
    111         nPts ++;
    112 
    113         ptCP.x = region->x0; ptCP.y = region->y1;
    114         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    115         xObs->data.F32[nPts] = ptFP.x;
    116         yObs->data.F32[nPts] = ptFP.y;
    117         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    118         xRef->data.F32[nPts] = ptFP.x;
    119         yRef->data.F32[nPts] = ptFP.y;
    120 
    121         if (DEBUG) {
    122             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    123             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    124         }
    125         nPts ++;
    126 
    127         ptCP.x = region->x1; ptCP.y = region->y1;
    128         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    129         xObs->data.F32[nPts] = ptFP.x;
    130         yObs->data.F32[nPts] = ptFP.y;
    131         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    132         xRef->data.F32[nPts] = ptFP.x;
    133         yRef->data.F32[nPts] = ptFP.y;
    134 
    135         if (DEBUG) {
    136             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    137             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    138         }
    139         nPts ++;
    140 
    141         ptCP.x = region->x1; ptCP.y = region->y0;
    142         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    143         xObs->data.F32[nPts] = ptFP.x;
    144         yObs->data.F32[nPts] = ptFP.y;
    145         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    146         xRef->data.F32[nPts] = ptFP.x;
    147         yRef->data.F32[nPts] = ptFP.y;
    148 
    149         if (DEBUG) {
    150             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    151             fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
    152         }
    153         nPts ++;
    154 
    155         psFree (region);
     72        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     73
     74        // XXX we are currently inconsistent with marking the good vs the bad data
     75        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     76        // probably wrong since it implies there is no data!
     77
     78        // skip chips for which the astrometry failed (NASTRO == 0)
     79        if (!obsChip->cells->n) continue;
     80        pmCell *cell = obsChip->cells->data[0];
     81        if (!cell) continue;
     82
     83        if (!cell->readouts->n) continue;
     84        pmReadout *readout = cell->readouts->data[0];
     85        if (!readout) continue;
     86
     87        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     88        if (!updates) continue;
     89
     90        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     91        if (!nAstro) continue;
     92
     93        // set the chip astrometry using the astrom file
     94        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     95
     96        psRegion *region = pmChipPixels (obsChip);
     97        psPlane ptCP, ptFP;
     98
     99        ptCP.x = region->x0; ptCP.y = region->y0;
     100        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     101        xObs->data.F32[nPts] = ptFP.x;
     102        yObs->data.F32[nPts] = ptFP.y;
     103        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     104        xRef->data.F32[nPts] = ptFP.x;
     105        yRef->data.F32[nPts] = ptFP.y;
     106
     107        if (DEBUG) {
     108            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     109            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     110        }
     111        nPts ++;
     112
     113        ptCP.x = region->x0; ptCP.y = region->y1;
     114        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     115        xObs->data.F32[nPts] = ptFP.x;
     116        yObs->data.F32[nPts] = ptFP.y;
     117        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     118        xRef->data.F32[nPts] = ptFP.x;
     119        yRef->data.F32[nPts] = ptFP.y;
     120
     121        if (DEBUG) {
     122            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     123            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     124        }
     125        nPts ++;
     126
     127        ptCP.x = region->x1; ptCP.y = region->y1;
     128        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     129        xObs->data.F32[nPts] = ptFP.x;
     130        yObs->data.F32[nPts] = ptFP.y;
     131        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     132        xRef->data.F32[nPts] = ptFP.x;
     133        yRef->data.F32[nPts] = ptFP.y;
     134
     135        if (DEBUG) {
     136            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     137            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     138        }
     139        nPts ++;
     140
     141        ptCP.x = region->x1; ptCP.y = region->y0;
     142        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     143        xObs->data.F32[nPts] = ptFP.x;
     144        yObs->data.F32[nPts] = ptFP.y;
     145        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     146        xRef->data.F32[nPts] = ptFP.x;
     147        yRef->data.F32[nPts] = ptFP.y;
     148
     149        if (DEBUG) {
     150            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     151            fprintf (f, "%s  %f %f  %f %f\n", chipName, xObs->data.F32[nPts], yObs->data.F32[nPts], xRef->data.F32[nPts], yRef->data.F32[nPts]);
     152        }
     153        nPts ++;
     154
     155        psFree (region);
    156156    }
    157157    xObs->n = yObs->n = xRef->n = yRef->n = nPts;
    158158    if (DEBUG) fclose (f);
    159        
     159
    160160    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
    161  
     161
    162162    psVector *mask = psVectorAlloc (nPts, PS_TYPE_U8);
    163163    psVectorInit (mask, 0);
    … …  
    167167
    168168    for (int i = 0; i < 3; i++) {
    169         psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
    170         psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
    171 
    172         psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
    173         psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
     169        psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
     170        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
     171
     172        psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
     173        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
    174174    }
    175175
    176176    // loop over all chips, select the outliers, and replace the measured astrometry with the model
    177     // the measured transformation above must be applied to make the comparison, and also then applied to the 
     177    // the measured transformation above must be applied to make the comparison, and also then applied to the
    178178    // model transformation
    179179
    180180    if (DEBUG) {
    181         f = fopen ("corners.fit.dat", "w");
    182         for (int i = 0; i < xObs->n; i++) {
    183             psPlane obsCoord, refCoord;
    184             refCoord.x = xRef->data.F32[i];
    185             refCoord.y = yRef->data.F32[i];
    186             psPlaneTransformApply (&obsCoord, map, &refCoord);
    187             fprintf (f, "%f %f  %f %f  %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y);
    188         }
    189         fclose (f);
     181        f = fopen ("corners.fit.dat", "w");
     182        for (int i = 0; i < xObs->n; i++) {
     183            psPlane obsCoord, refCoord;
     184            refCoord.x = xRef->data.F32[i];
     185            refCoord.y = yRef->data.F32[i];
     186            psPlaneTransformApply (&obsCoord, map, &refCoord);
     187            fprintf (f, "%f %f  %f %f  %f %f\n", xObs->data.F32[i], yObs->data.F32[i], xRef->data.F32[i], yRef->data.F32[i], obsCoord.x, obsCoord.y);
     188        }
     189        fclose (f);
    190190    }
    191191
    … …  
    199199
    200200    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    201         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    202         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    203 
    204         // set the chip astrometry using the astrom file
    205         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    206 
    207         // bad Astrometry test:  ref pixel or angle outside nominal
    208 
    209         psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
    210         psPlane obsCoord, refCoord, tmpCoord;
    211 
    212         // find location of 0,0 pixel in focal plane coords for this chip
    213         psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
    214 
    215         // find location of 0,0 pixel in focal plane coords for ref chip
    216         // apply the global field rotation and offset before comparing
    217         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
    218         psPlaneTransformApply (&refCoord, map, &tmpCoord);
    219    
    220         psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
    221         psPlane obsOffPt, refOffPt;
    222 
    223         // find location of 0,0 pixel in focal plane coords for this chip
    224         psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
    225 
    226         // find location of 0,0 pixel in focal plane coords for ref chip
    227         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
    228         psPlaneTransformApply (&refOffPt, map, &tmpCoord);
    229    
    230         double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
    231         double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
    232 
    233         bool badAstrom = false;
    234         badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
    235         badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
    236         badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
    237 
    238         fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    239 
    240         // XXX for now, just use first readout
    241         pmCell *cell = obsChip->cells->data[0];
    242         pmReadout *readout = cell->readouts->data[0];
    243 
    244         // update the header (pull or create local view to entry on readout->analysis)
    245         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    246         if (!updates) {
    247             updates = psMetadataAlloc ();
    248             psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    249             psFree (updates);
    250         }
    251 
    252         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
    253         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
    254         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
    255 
    256         // for successful chips, save the measured offsets in the header
    257         if (!badAstrom) continue;
    258 
    259         // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
    260         // catch the failures relative to the model
    261         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
    262         continue;
    263 
    264         psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
    265                   view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    266 
    267         psFree (obsChip->toFPA);
    268         psFree (obsChip->fromFPA);
    269 
    270         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    271         // XXX this only works if toTPA is at most a linear transformation
    272         psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
    273         for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
    274             for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
    275                 double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    276                 double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    277                 toFPA->x->coeff[i][j] = f1 + f2;
    278 
    279                 double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    280                 double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    281                 toFPA->y->coeff[i][j] = g1 + g2;
    282             }
    283         }
    284         toFPA->x->coeff[0][0] += map->x->coeff[0][0];
    285         toFPA->y->coeff[0][0] += map->y->coeff[0][0];
    286 
    287         psRegion *region = pmChipPixels (obsChip);
    288         obsChip->toFPA   = toFPA;
    289         obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
    290         psFree (region);
    291    
    292         // use the new position to re-try the match fit
    293         // select the raw objects for this readout
    294         psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    295         if (rawstars == NULL) { continue; }
    296 
    297         // select the raw objects for this readout
    298         psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    299         if (refstars == NULL) { continue; }
    300 
    301         // the absolute minimum number of stars is 4 (for order = 1)
    302         if ((rawstars->n < 4) || (refstars->n < 4)) {
    303             readout->data_exists = false;
    304             psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
    305                       rawstars->n, refstars->n);
    306             continue;
    307         }
    308 
    309         psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
    310 
    311         // XXX update the header with info to reflect the failure
    312         if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
    313             readout->data_exists = false;
    314             psLogMsg ("psastro", 3, "failed to find a solution\n");
    315             continue;
    316         }
    317 
    318         pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
    319     }
    320 
    321     psastroVisualPlotFixChips (input, xObs, yObs);
     201        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
     202        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     203
     204        // set the chip astrometry using the astrom file
     205        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     206
     207        // bad Astrometry test:  ref pixel or angle outside nominal
     208
     209        psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
     210        psPlane obsCoord, refCoord, tmpCoord;
     211
     212        // find location of 0,0 pixel in focal plane coords for this chip
     213        psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
     214
     215        // find location of 0,0 pixel in focal plane coords for ref chip
     216        // apply the global field rotation and offset before comparing
     217        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
     218        psPlaneTransformApply (&refCoord, map, &tmpCoord);
     219
     220        psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
     221        psPlane obsOffPt, refOffPt;
     222
     223        // find location of 0,0 pixel in focal plane coords for this chip
     224        psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
     225
     226        // find location of 0,0 pixel in focal plane coords for ref chip
     227        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
     228        psPlaneTransformApply (&refOffPt, map, &tmpCoord);
     229
     230        double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
     231        double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
     232
     233        bool badAstrom = false;
     234        badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
     235        badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
     236        badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
     237
     238        fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     239
     240        // XXX for now, just use first readout
     241        pmCell *cell = obsChip->cells->data[0];
     242        pmReadout *readout = cell->readouts->data[0];
     243
     244        // update the header (pull or create local view to entry on readout->analysis)
     245        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     246        if (!updates) {
     247            updates = psMetadataAlloc ();
     248            psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     249            psFree (updates);
     250        }
     251
     252        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
     253        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
     254        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
     255
     256        // for successful chips, save the measured offsets in the header
     257        if (!badAstrom) continue;
     258
     259        // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
     260        // catch the failures relative to the model
     261        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
     262        continue;
     263
     264        psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
     265                  view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     266
     267        psFree (obsChip->toFPA);
     268        psFree (obsChip->fromFPA);
     269
     270        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     271        // XXX this only works if toTPA is at most a linear transformation
     272        psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
     273        for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
     274            for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
     275                double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     276                double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     277                toFPA->x->coeff[i][j] = f1 + f2;
     278
     279                double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     280                double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     281                toFPA->y->coeff[i][j] = g1 + g2;
     282            }
     283        }
     284        toFPA->x->coeff[0][0] += map->x->coeff[0][0];
     285        toFPA->y->coeff[0][0] += map->y->coeff[0][0];
     286
     287        psRegion *region = pmChipPixels (obsChip);
     288        obsChip->toFPA   = toFPA;
     289        obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
     290        psFree (region);
     291
     292        // use the new position to re-try the match fit
     293        // select the raw objects for this readout
     294        psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     295        if (rawstars == NULL) { continue; }
     296
     297        // select the raw objects for this readout
     298        psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     299        if (refstars == NULL) { continue; }
     300
     301        // the absolute minimum number of stars is 4 (for order = 1)
     302        if ((rawstars->n < 4) || (refstars->n < 4)) {
     303            readout->data_exists = false;
     304            psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
     305                      rawstars->n, refstars->n);
     306            continue;
     307        }
     308
     309        psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
     310
     311        // XXX update the header with info to reflect the failure
     312        if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
     313            readout->data_exists = false;
     314            psLogMsg ("psastro", 3, "failed to find a solution\n");
     315            continue;
     316        }
     317
     318        pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
     319    }
     320
     321    pmAstromVisualPlotFixChips (input, xObs, yObs);
    322322    psFree (xObs);
    323323    psFree (yObs);

  • branches/cnb_branch_20090113/psastro/src/psastroLoadRefstars.c

    r20805 r21208  
    138138    }
    139139
    140     psastroVisualPlotRefStars (refstars, recipe);
     140    pmAstromVisualPlotRefStars (refstars, recipe);
    141141
    142142    if (psTraceGetLevel("psastro.plot") > 0) {
    … …  
    239239    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    240240    if (!input) {
    241         psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
    242         photcode = psStringCopy ("NONE");
    243         return photcode;
     241        psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
     242        photcode = psStringCopy ("NONE");
     243        return photcode;
    244244    }
    245245    assert (input->fpa);
    … …  
    247247    *maxRho = psMetadataLookupF32(&status, recipe, "DVO.GETSTAR.MAX.RHO");
    248248    if (!status) {
    249         psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
    250         return NULL;
     249        psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
     250        return NULL;
    251251    }
    252252
    … …  
    266266    if (!status) ESCAPE ("missing DVO.GETSTAR.MIN.MAG.INST");
    267267
    268     // PHOTCODE.DATA is a multi of metadata items 
     268    // PHOTCODE.DATA is a multi of metadata items
    269269    psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false);
    270270
    271271    psMetadataItem *refItem = NULL;
    272272    while ((refItem = psListGetAndIncrement (iter))) {
    273         if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
    274    
    275         char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
    276         if (!status) {
    277             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    278             continue;
    279         }
    280 
    281         // does this entry match the current filter?
    282         if (strcmp (refFilter, filter)) continue;
    283 
    284         psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
    285 
    286         float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
    287         if (!status) {
    288             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    289             continue;
    290         }
    291         photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
    292         if (!status) {
    293             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    294             continue;
    295         }
    296 
    297         // convert the minInst to a calibrated minimum magnitude
    298         *minMag = minInst + 2.5*log10(exptime) + zeropt;
    299 
    300         psFree (iter);
    301         return photcode;
    302     }   
     273        if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
     274
     275        char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
     276        if (!status) {
     277            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     278            continue;
     279        }
     280
     281        // does this entry match the current filter?
     282        if (strcmp (refFilter, filter)) continue;
     283
     284        psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
     285
     286        float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
     287        if (!status) {
     288            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     289            continue;
     290        }
     291        photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
     292        if (!status) {
     293            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     294            continue;
     295        }
     296
     297        // convert the minInst to a calibrated minimum magnitude
     298        *minMag = minInst + 2.5*log10(exptime) + zeropt;
     299
     300        psFree (iter);
     301        return photcode;
     302    }
    303303    psFree (iter);
    304304
    … …  
    306306    photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE");
    307307    PS_ASSERT (photcode, NULL);
    308        
     308
    309309    // give up and use fixed value
    310310    *minMag = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MIN.MAG");

  • branches/cnb_branch_20090113/psastro/src/psastroLuminosityFunction.c

    r20805 r21208  
    129129    lumFunc->sPeak = sPeak;
    130130
    131     psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
     131    pmAstromVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
    132132
    133133    psFree (lnMag);

  • branches/cnb_branch_20090113/psastro/src/psastroMosaicOneChip.c

    r20803 r21208  
    148148
    149149    //plot results
    150     psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
     150    pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
    151151
    152152    psFree (fitStats);

  • branches/cnb_branch_20090113/psastro/src/psastroMosaicSetMatch.c

    r20805 r21208  
    5757                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    5858
    59                 psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
     59                pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
    6060
    6161                // XXX drop the old one

  • branches/cnb_branch_20090113/psastro/src/psastroOneChipFit.c

    r20805 r21208  
    1111
    1212    // default value for match/fit : radius is in pixels
    13     REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
     13    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32);
    1414
    1515    // run the match/fit sequence NITER times
    16     REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
     16    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
    1717
    1818    // correct radius to FP units (physical pixel scale in microns per pixel)
    19     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     19    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
    2020    RADIUS *= pixelScale;
    2121
    … …  
    3232
    3333    for (int iter = 0; iter < nIter; iter++) {
    34        
    35         char name[128];
    3634
    37         sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
    38         float radius = psMetadataLookupF32 (&status, recipe, name);
    39         radius *= pixelScale;
    40         if (!status || (radius == 0.0)) {
    41             radius = RADIUS;
    42         }
     35        char name[128];
     36
     37        sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
     38        float radius = psMetadataLookupF32 (&status, recipe, name);
     39        radius *= pixelScale;
     40        if (!status || (radius == 0.0)) {
     41            radius = RADIUS;
     42        }
    4343
    4444
    45         // use small radius to match stars
    46         match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
    47         if (match == NULL) {
    48             psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
    49             return false;
    50         }
     45        // use small radius to match stars
     46        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
     47        if (match == NULL) {
     48            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
     49            return false;
     50        }
    5151
    52         // modify the order to correspond to the actual number of matched stars:
    53         int Ndof_min = 3;
    54         int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
    55         order = PS_MIN (order, order_max);
     52        // modify the order to correspond to the actual number of matched stars:
     53        int Ndof_min = 3;
     54        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
     55        order = PS_MIN (order, order_max);
    5656
    57         // if ((match->n < 11) && (order >= 3)) order = 2;
    58         // if ((match->n <  7) && (order >= 2)) order = 1;
    59         // if ((match->n <  4) && (order >= 1)) order = 0;
     57        // if ((match->n < 11) && (order >= 3)) order = 2;
     58        // if ((match->n <  7) && (order >= 2)) order = 1;
     59        // if ((match->n <  4) && (order >= 1)) order = 0;
    6060
    61         if (order < 1) {
    62             psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
    63             psFree (match);
    64             return false;
    65         }
     61        if (order < 1) {
     62            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
     63            psFree (match);
     64            return false;
     65        }
    6666
    67         // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    68         psFree (chip->toFPA);
    69         chip->toFPA = psPlaneTransformAlloc (order, order);
    70         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    71             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    72                 if (i + j > order) {
    73                     chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    74                     chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    75                 }
    76             }
    77         }
     67        // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     68        psFree (chip->toFPA);
     69        chip->toFPA = psPlaneTransformAlloc (order, order);
     70        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     71            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     72                if (i + j > order) {
     73                    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     74                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     75                }
     76            }
     77        }
    7878
    79         // XXX allow statitic to be set by the user
    80         // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    81         fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    82         fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
    83         fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
     79        // XXX allow statitic to be set by the user
     80        // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     81        fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     82        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
     83        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
    8484
    85         // improved fit for astrometric terms
    86         results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    87         if (!results) {
    88             psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
    89             psFree (match);
    90             psFree (fitStats);
    91             return false;
    92         }
    93    
    94         // determine fromFPA transformation and apply new transformation to raw & ref stars
    95         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    96    
    97         // toSky converts from FPA & TPA units (microns) to sky units (radians)
    98         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    99         // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    100         float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
    101         int astNstar = results->yStats->clippedNvalues;
    102         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     85        // improved fit for astrometric terms
     86        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     87        if (!results) {
     88            psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
     89            psFree (match);
     90            psFree (fitStats);
     91            return false;
     92        }
    10393
    104         if (iter < nIter - 1) {
    105             psFree (fitStats);
    106             psFree (results);
    107             psFree (match);
    108         }
     94        // determine fromFPA transformation and apply new transformation to raw & ref stars
     95        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     96
     97        // toSky converts from FPA & TPA units (microns) to sky units (radians)
     98        float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
     99        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     100        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     101        int astNstar = results->yStats->clippedNvalues;
     102        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     103
     104        if (iter < nIter - 1) {
     105            psFree (fitStats);
     106            psFree (results);
     107            psFree (match);
     108        }
    109109    }
    110110
    … …  
    127127    if (astError > maxError) {
    128128        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    129         validSolution = false;
     129        validSolution = false;
    130130    }
    131131    if (astNstar < minNstar) {
    132132        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    133         validSolution = false;
     133        validSolution = false;
    134134    }
    135135
    … …  
    138138    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    139139    if (validSolution) {
    140         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    141         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     140        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     141        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    142142    } else {
    143         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    144         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     143        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     144        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    145145    }
    146146    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    … …  
    148148    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
    149149    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    150    
     150
    151151    // XXX check if we correctly applied the new transformation:
    152152    if (psTraceGetLevel("psastro.dump") > 0) {
    153         psastroDumpRawstars (rawstars, fpa, chip);
    154         psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
    155         psastroDumpStars (refstars, "refstars.cal.dat");
     153        psastroDumpRawstars (rawstars, fpa, chip);
     154        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
     155        psastroDumpStars (refstars, "refstars.cal.dat");
    156156    }
    157157
    158     psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);
     158    pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe);
    159159
    160160    if (psTraceGetLevel("psastro.plot") > 0) {
    161         psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     161        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
    162162    }
    163163

  • branches/cnb_branch_20090113/psastro/src/psastroRemoveClumps.c

    r20805 r21208  
    5555    psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit);
    5656
    57     psastroVisualPlotRemoveClumps (input, count, scale, limit);
     57    pmAstromVisualPlotRemoveClumps (input, count, scale, limit);
    5858
    5959    // find and exclude objects in bad pixels

  • branches/cnb_branch_20090113/psastro/src/psastroUtils.c

    r20804 r21208  
    4242            if (!chip->toFPA) { continue; }
    4343
    44             if (chip->cells->n == 0) { continue; }
    45             pmCell *cell = chip->cells->data[0];
     44            if (chip->cells->n == 0) { continue; }
     45            pmCell *cell = chip->cells->data[0];
    4646            if (!cell->process || !cell->file_exists) { continue; }
    4747
    48             if (cell->readouts->n == 0) { continue; }
    49             pmReadout *readout = cell->readouts->data[0];
    50             if (! readout->data_exists) { continue; }
     48            if (cell->readouts->n == 0) { continue; }
     49            pmReadout *readout = cell->readouts->data[0];
     50            if (! readout->data_exists) { continue; }
    5151
    5252            pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
    … …  
    8888    psastroMosaicSetAstrom (fpa);
    8989    if (!useExternal) {
    90         psastroVisualPlotCommonScale (fpa, oldScale);
     90        pmAstromVisualPlotCommonScale (fpa, oldScale);
    9191    }
    9292    psFree (oldScale);
Note: See TracChangeset for help on using the changeset viewer.