IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21422 for trunk/psastro/src


Ignore:
Timestamp:
Feb 9, 2009, 11:25:34 AM (17 years ago)
Author:
Paul Price
Message:

Merging cnb_branch_20090113. No major conflicts. Compiles, but not tested.

Location:
trunk/psastro/src
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/Makefile.am

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

  • trunk/psastro/src/psastro.h

    r21409 r21422  
    11/** @file psastro.h
    22 *
    3  *  @brief This file defines the library functions available to external 
    4  *  programs. 
     3 *  @brief This file defines the library functions available to external
     4 *  programs.
    55 *
    6  *  It must be included by programs which are compiled against 
     6 *  It must be included by programs which are compiled against
    77 *  psphot functions.
    88 *
    … …  
    1010 *
    1111 *  @author IfA
    12  *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2009-02-07 02:03:34 $
     12 *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2009-02-09 21:25:34 $
    1414 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1515 */
    … …  
    3030# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
    3131
     32#if 0
    3233/**
    3334 * this structure represents a fit to the logN / logS curve for a set of stars
    … …  
    4344    int sPeak;                          ///< sum of stars to peak bin
    4445} pmLumFunc;
     46#endif
    4547
    4648bool              psastroDataSave (pmConfig *config);
    … …  
    9597psString          psastroVersionLong(void);
    9698
    97 // psastroVisual functions
    98 bool psastroSetVisual (bool mode);
    99 bool psastroVisualClose();
    100 bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
    101 bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
    102 bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe);
    103 bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit);
    104 bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld);
    105 bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    106 bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd);
    107 bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
    108 bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale);
    109 bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe);
    110 
    11199// demo plots
    112100bool              psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);

  • trunk/psastro/src/psastroArguments.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    100100    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    101101        psArgumentRemove (N, &argc, argv);
    102         psastroSetVisual (true);
    103102        pmAstromSetVisual (true);
    104103    }

  • trunk/psastro/src/psastroAstromGuess.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    4242    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    4343    if (!recipe) {
    44         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
    45         return false;
     44        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
     45        return false;
    4646    }
    4747
    … …  
    4949    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
    5050    if (!status) {
    51         useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     51        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
    5252    }
    5353
    … …  
    5555    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    5656    if (!input) {
    57         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    58         return false;
     57        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     58        return false;
    5959    }
    6060
    … …  
    6262    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    6363    if (!status) {
    64         psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
    65         return false;
    66     } 
     64        psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
     65        return false;
     66    }
    6767
    6868    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    8080    bool bilevelAstrometry = false;
    8181    if (!useModel) {
    82         psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
     82        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    8383    }
    8484
    … …  
    8888        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    8989
    90         if (!useModel) {
    91             if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
    92         }
     90        if (!useModel) {
     91            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
     92        }
    9393
    9494        if (newFPA) {
    9595            newFPA = false;
    96             while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
    97             while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
     96            while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
     97            while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
    9898            RAminSky = fpa->toSky->R - M_PI;
    9999            RAmaxSky = fpa->toSky->R + M_PI;
    100100        }
    101101
    102         // report and save the current best guess for the chip 0,0 pixel coordinates
    103         {
    104             psPlane ptCH, ptFP, ptTP;
    105             psSphere ptSky;
    106 
    107             ptCH.x = 0;
    108             ptCH.y = 0;
    109             psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    110             psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    111             psDeproject (&ptSky, &ptTP, fpa->toSky);
    112             psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    113 
    114             psVectorAppend (cornerL, ptFP.x);
    115             psVectorAppend (cornerM, ptFP.y);
    116             psVectorAppend (cornerP, ptTP.x);
    117             psVectorAppend (cornerQ, ptTP.y);
    118             psVectorAppend (cornerR, ptSky.r);
    119             psVectorAppend (cornerD, ptSky.d);
    120         }
     102        // report and save the current best guess for the chip 0,0 pixel coordinates
     103        {
     104            psPlane ptCH, ptFP, ptTP;
     105            psSphere ptSky;
     106
     107            ptCH.x = 0;
     108            ptCH.y = 0;
     109            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     110            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     111            psDeproject (&ptSky, &ptTP, fpa->toSky);
     112            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     113
     114            psVectorAppend (cornerL, ptFP.x);
     115            psVectorAppend (cornerM, ptFP.y);
     116            psVectorAppend (cornerP, ptTP.x);
     117            psVectorAppend (cornerQ, ptTP.y);
     118            psVectorAppend (cornerR, ptSky.r);
     119            psVectorAppend (cornerD, ptSky.d);
     120        }
    121121
    122122        // apply the new WCS guess data to all of the data in the readouts
    … …  
    132132                if (rawstars == NULL) { continue; }
    133133
    134                 *nStars += rawstars->n;
     134                *nStars += rawstars->n;
    135135                for (int i = 0; i < rawstars->n; i++) {
    136136                    pmAstromObj *raw = rawstars->data[i];
    … …  
    151151                }
    152152
    153                 // dump or plot the resulting projected positions
    154                 if (psTraceGetLevel("psastro.dump") > 0) {
    155                     psastroDumpRawstars (rawstars, fpa, chip);
    156                 }
    157 
    158                 psastroVisualPlotRawStars(rawstars, fpa, chip, recipe);
    159 
    160                 if (psTraceGetLevel("psastro.plot") > 0) {
    161                     psastroPlotRawstars (rawstars, fpa, chip, recipe);
    162                 }
     153                // dump or plot the resulting projected positions
     154                if (psTraceGetLevel("psastro.dump") > 0) {
     155                    psastroDumpRawstars (rawstars, fpa, chip);
     156                }
     157
     158                pmAstromVisualPlotRawStars(rawstars, fpa, chip, recipe);
     159
     160                if (psTraceGetLevel("psastro.plot") > 0) {
     161                    psastroPlotRawstars (rawstars, fpa, chip, recipe);
     162                }
    163163            }
    164164        }
    … …  
    170170    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
    171171    if (*nStars == 0) {
    172         psLogMsg ("psastro", 2, "no sources available for astrometry\n");
    173         psFree (view);
    174         return true;
    175     }
    176 
    177     psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
    178               DEG_RAD*RAmin, DEG_RAD*DECmin,
    179               DEG_RAD*RAmax, DEG_RAD*DECmax);
     172        psLogMsg ("psastro", 2, "no sources available for astrometry\n");
     173        psFree (view);
     174        return true;
     175    }
     176
     177    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
     178              DEG_RAD*RAmin, DEG_RAD*DECmin,
     179              DEG_RAD*RAmax, DEG_RAD*DECmax);
    180180
    181181    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
    … …  
    216216    pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    217217    if (bilevelAstrometry) {
    218         if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    219             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    220             return false;
    221         }
     218        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
     219            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     220            return false;
     221        }
    222222    } else {
    223         if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
    224             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    225             return false;
    226         }
     223        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
     224            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     225            return false;
     226        }
    227227    }
    228228    return true;
    … …  
    238238    // load mosaic-level astrometry?
    239239    if (phu) {
    240         char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    241         if (ctype) {
    242             *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    243         }
     240        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     241        if (ctype) {
     242            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     243        }
    244244    }
    245245    if (*bilevelAstrometry) {
    246         pmAstromReadBilevelMosaic (fpa, phu->header);
    247     } 
     246        pmAstromReadBilevelMosaic (fpa, phu->header);
     247    }
    248248    psFree (view);
    249249    return true;
    … …  
    258258    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    259259    if (!input) {
    260         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    261         return false;
     260        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     261        return false;
    262262    }
    263263
    … …  
    290290        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    291291
    292         // XXX we are currently inconsistent with marking the good vs the bad data
    293         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    294         // probably wrong since it implies there is no data!
    295 
    296         // skip chips for which the astrometry failed (NASTRO == 0)
    297         if (!chip->cells->n) goto skip_chip;
    298         pmCell *cell = chip->cells->data[0];
    299         if (!cell) goto skip_chip;
    300 
    301         if (!cell->readouts->n) goto skip_chip;
    302         pmReadout *readout = cell->readouts->data[0];
    303         if (!readout) goto skip_chip;
    304 
    305         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    306         if (!updates) goto skip_chip;
    307        
    308         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    309         if (!nAstro) goto skip_chip;
    310 
    311         float astError = psMetadataLookupF32 (&status, updates, "CERROR");
    312         if (fabs(astError) < 1e-6) goto skip_chip;
    313 
    314         psPlane ptCH, ptFP, ptTP;
    315         psSphere ptSky;
    316 
    317         ptCH.x = 0;
    318         ptCH.y = 0;
    319         psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    320         psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    321         psDeproject (&ptSky, &ptTP, fpa->toSky);
    322         psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    323 
    324         // new corner locations based on the calibrated astrometry
    325         psVectorAppend (cornerLn, ptFP.x);
    326         psVectorAppend (cornerMn, ptFP.y);
    327         psVectorAppend (cornerPn, ptTP.x);
    328         psVectorAppend (cornerQn, ptTP.y);
    329         psVectorAppend (cornerRn, ptSky.r);
    330         psVectorAppend (cornerDn, ptSky.d);
    331         psVectorAppend (cornerMK, 0);
    332         continue;
     292        // XXX we are currently inconsistent with marking the good vs the bad data
     293        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     294        // probably wrong since it implies there is no data!
     295
     296        // skip chips for which the astrometry failed (NASTRO == 0)
     297        if (!chip->cells->n) goto skip_chip;
     298        pmCell *cell = chip->cells->data[0];
     299        if (!cell) goto skip_chip;
     300
     301        if (!cell->readouts->n) goto skip_chip;
     302        pmReadout *readout = cell->readouts->data[0];
     303        if (!readout) goto skip_chip;
     304
     305        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     306        if (!updates) goto skip_chip;
     307
     308        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     309        if (!nAstro) goto skip_chip;
     310
     311        float astError = psMetadataLookupF32 (&status, updates, "CERROR");
     312        if (fabs(astError) < 1e-6) goto skip_chip;
     313
     314        psPlane ptCH, ptFP, ptTP;
     315        psSphere ptSky;
     316
     317        ptCH.x = 0;
     318        ptCH.y = 0;
     319        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     320        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     321        psDeproject (&ptSky, &ptTP, fpa->toSky);
     322        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     323
     324        // new corner locations based on the calibrated astrometry
     325        psVectorAppend (cornerLn, ptFP.x);
     326        psVectorAppend (cornerMn, ptFP.y);
     327        psVectorAppend (cornerPn, ptTP.x);
     328        psVectorAppend (cornerQn, ptTP.y);
     329        psVectorAppend (cornerRn, ptSky.r);
     330        psVectorAppend (cornerDn, ptSky.d);
     331        psVectorAppend (cornerMK, 0);
     332        continue;
    333333
    334334    skip_chip:
    335         // new corner locations based on the calibrated astrometry
    336         psVectorAppend (cornerLn, 0.0);
    337         psVectorAppend (cornerMn, 0.0);
    338         psVectorAppend (cornerPn, 0.0);
    339         psVectorAppend (cornerQn, 0.0);
    340         psVectorAppend (cornerRn, 0.0);
    341         psVectorAppend (cornerDn, 0.0);
    342         psVectorAppend (cornerMK, 1);
     335        // new corner locations based on the calibrated astrometry
     336        psVectorAppend (cornerLn, 0.0);
     337        psVectorAppend (cornerMn, 0.0);
     338        psVectorAppend (cornerPn, 0.0);
     339        psVectorAppend (cornerQn, 0.0);
     340        psVectorAppend (cornerRn, 0.0);
     341        psVectorAppend (cornerDn, 0.0);
     342        psVectorAppend (cornerMK, 1);
    343343    }
    344344
    … …  
    349349
    350350    for (int i = 0; i < cornerRo->n; i++) {
    351        
    352         psPlane ptTP;
    353         psSphere ptSky;
    354 
    355         ptSky.r = cornerRo->data.F32[i];
    356         ptSky.d = cornerDo->data.F32[i];
    357 
    358         psProject (&ptTP, &ptSky, fpa->toSky);
    359         psVectorAppend (cornerPs, ptTP.x);
    360         psVectorAppend (cornerQs, ptTP.y);
     351
     352        psPlane ptTP;
     353        psSphere ptSky;
     354
     355        ptSky.r = cornerRo->data.F32[i];
     356        ptSky.d = cornerDo->data.F32[i];
     357
     358        psProject (&ptTP, &ptSky, fpa->toSky);
     359        psVectorAppend (cornerPs, ptTP.x);
     360        psVectorAppend (cornerQs, ptTP.y);
    361361    }
    362362
    … …  
    364364    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    365365    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
    366    
     366
    367367    // fit the valid chips, mask the invalid chips
    368368    psVectorFitPolynomial2D (map->x, cornerMK, 1, cornerPn, NULL, cornerPs, cornerQs);
    369369    psVectorFitPolynomial2D (map->y, cornerMK, 1, cornerQn, NULL, cornerPs, cornerQs);
    370    
     370
    371371    // apply the linear fit...
    372372    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
    … …  
    377377    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
    378378
    379     psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
     379    pmAstromVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
    380380
    381381    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    … …  
    387387    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
    388388    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
    389    
     389
    390390    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
    391391    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
    … …  
    395395    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
    396396    if (!header) {
    397         header = psMetadataAlloc();
    398         psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
    399         psFree (header);  // drop this reference
     397        header = psMetadataAlloc();
     398        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
     399        psFree (header);  // drop this reference
    400400    }
    401401
    … …  
    408408
    409409    if (DEBUG) {
    410         FILE *f = fopen ("corners.dat", "w");
    411         for (int i = 0; i < cornerRo->n; i++) {
    412             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",
    413                      cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
    414                      cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
    415         }
    416         fclose (f);
     410        FILE *f = fopen ("corners.dat", "w");
     411        for (int i = 0; i < cornerRo->n; i++) {
     412            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",
     413                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
     414                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
     415        }
     416        fclose (f);
    417417    }
    418418
    … …  
    437437    psFree (map);
    438438    psFree (view);
    439    
     439
    440440
    441441    return true;

  • trunk/psastro/src/psastroCleanup.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    1616
    1717    psFree (config);
     18    pmAstromVisualClose ();
    1819
    1920    psTimerStop ();
    … …  
    2324    pmConceptsDone ();
    2425    pmConfigDone ();
    25     psastroVisualClose ();
    26     pmAstromVisualClose ();
    2726    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    2827    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");

  • trunk/psastro/src/psastroFixChips.c

    r21409 r21422  
    11/** @file psastroFixChips.c
    22 *
    3  *  @brief 
     3 *  @brief
    44 *
    55 *  @ingroup libpsastro
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    2525    bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
    2626    if (!status) {
    27         fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     27        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
    2828    }
    2929    if (!fixChips) return true;
    … …  
    3737    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    3838    if (!input) {
    39         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    40         return false;
     39        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     40        return false;
    4141    }
    4242
    … …  
    6060    // files associated with the science image
    6161    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    62         psError (PS_ERR_IO, false, "Can't load the astrometry model file");
    63         return false;
     62        psError (PS_ERR_IO, false, "Can't load the astrometry model file");
     63        return false;
    6464    }
    6565
    … …  
    7676
    7777    if (DEBUG) {
    78         f = fopen ("corners.raw.dat", "w");
    79         chipName = NULL;
     78        f = fopen ("corners.raw.dat", "w");
     79        chipName = NULL;
    8080    }
    8181
    8282    pmChip *obsChip = NULL;
    8383    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    84         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    85 
    86         // XXX we are currently inconsistent with marking the good vs the bad data
    87         // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
    88         // probably wrong since it implies there is no data!
    89 
    90         // skip chips for which the astrometry failed (NASTRO == 0)
    91         if (!obsChip->cells->n) continue;
    92         pmCell *cell = obsChip->cells->data[0];
    93         if (!cell) continue;
    94 
    95         if (!cell->readouts->n) continue;
    96         pmReadout *readout = cell->readouts->data[0];
    97         if (!readout) continue;
    98 
    99         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    100         if (!updates) continue;
    101        
    102         int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    103         if (!nAstro) continue;
    104 
    105         // set the chip astrometry using the astrom file
    106         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    107 
    108         psRegion *region = pmChipPixels (obsChip);
    109         psPlane ptCP, ptFP;
    110 
    111         ptCP.x = region->x0; ptCP.y = region->y0;
    112         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    113         xObs->data.F32[nPts] = ptFP.x;
    114         yObs->data.F32[nPts] = ptFP.y;
    115         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    116         xRef->data.F32[nPts] = ptFP.x;
    117         yRef->data.F32[nPts] = ptFP.y;
    118 
    119         if (DEBUG) {
    120             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    121             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]);
    122         }
    123         nPts ++;
    124 
    125         ptCP.x = region->x0; ptCP.y = region->y1;
    126         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    127         xObs->data.F32[nPts] = ptFP.x;
    128         yObs->data.F32[nPts] = ptFP.y;
    129         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    130         xRef->data.F32[nPts] = ptFP.x;
    131         yRef->data.F32[nPts] = ptFP.y;
    132 
    133         if (DEBUG) {
    134             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    135             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]);
    136         }
    137         nPts ++;
    138 
    139         ptCP.x = region->x1; ptCP.y = region->y1;
    140         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    141         xObs->data.F32[nPts] = ptFP.x;
    142         yObs->data.F32[nPts] = ptFP.y;
    143         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    144         xRef->data.F32[nPts] = ptFP.x;
    145         yRef->data.F32[nPts] = ptFP.y;
    146 
    147         if (DEBUG) {
    148             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    149             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]);
    150         }
    151         nPts ++;
    152 
    153         ptCP.x = region->x1; ptCP.y = region->y0;
    154         psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
    155         xObs->data.F32[nPts] = ptFP.x;
    156         yObs->data.F32[nPts] = ptFP.y;
    157         psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
    158         xRef->data.F32[nPts] = ptFP.x;
    159         yRef->data.F32[nPts] = ptFP.y;
    160 
    161         if (DEBUG) {
    162             chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
    163             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]);
    164         }
    165         nPts ++;
    166 
    167         psFree (region);
     84        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     85
     86        // XXX we are currently inconsistent with marking the good vs the bad data
     87        // psastroChipAstrom sets data_exists to false if the fit is bad.  this is
     88        // probably wrong since it implies there is no data!
     89
     90        // skip chips for which the astrometry failed (NASTRO == 0)
     91        if (!obsChip->cells->n) continue;
     92        pmCell *cell = obsChip->cells->data[0];
     93        if (!cell) continue;
     94
     95        if (!cell->readouts->n) continue;
     96        pmReadout *readout = cell->readouts->data[0];
     97        if (!readout) continue;
     98
     99        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     100        if (!updates) continue;
     101
     102        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
     103        if (!nAstro) continue;
     104
     105        // set the chip astrometry using the astrom file
     106        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     107
     108        psRegion *region = pmChipPixels (obsChip);
     109        psPlane ptCP, ptFP;
     110
     111        ptCP.x = region->x0; ptCP.y = region->y0;
     112        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     113        xObs->data.F32[nPts] = ptFP.x;
     114        yObs->data.F32[nPts] = ptFP.y;
     115        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     116        xRef->data.F32[nPts] = ptFP.x;
     117        yRef->data.F32[nPts] = ptFP.y;
     118
     119        if (DEBUG) {
     120            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     121            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]);
     122        }
     123        nPts ++;
     124
     125        ptCP.x = region->x0; ptCP.y = region->y1;
     126        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     127        xObs->data.F32[nPts] = ptFP.x;
     128        yObs->data.F32[nPts] = ptFP.y;
     129        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     130        xRef->data.F32[nPts] = ptFP.x;
     131        yRef->data.F32[nPts] = ptFP.y;
     132
     133        if (DEBUG) {
     134            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     135            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]);
     136        }
     137        nPts ++;
     138
     139        ptCP.x = region->x1; ptCP.y = region->y1;
     140        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     141        xObs->data.F32[nPts] = ptFP.x;
     142        yObs->data.F32[nPts] = ptFP.y;
     143        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     144        xRef->data.F32[nPts] = ptFP.x;
     145        yRef->data.F32[nPts] = ptFP.y;
     146
     147        if (DEBUG) {
     148            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     149            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]);
     150        }
     151        nPts ++;
     152
     153        ptCP.x = region->x1; ptCP.y = region->y0;
     154        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
     155        xObs->data.F32[nPts] = ptFP.x;
     156        yObs->data.F32[nPts] = ptFP.y;
     157        psPlaneTransformApply (&ptFP, refChip->toFPA, &ptCP);
     158        xRef->data.F32[nPts] = ptFP.x;
     159        yRef->data.F32[nPts] = ptFP.y;
     160
     161        if (DEBUG) {
     162            chipName = psMetadataLookupStr(NULL, obsChip->concepts, "CHIP.NAME");
     163            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]);
     164        }
     165        nPts ++;
     166
     167        psFree (region);
    168168    }
    169169    xObs->n = yObs->n = xRef->n = yRef->n = nPts;
    170170    if (DEBUG) fclose (f);
    171        
     171
    172172    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
    173  
     173
    174174    psVector *mask = psVectorAlloc (nPts, PS_TYPE_VECTOR_MASK);
    175175    psVectorInit (mask, 0);
    … …  
    179179
    180180    for (int i = 0; i < 3; i++) {
    181         psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
    182         psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
    183 
    184         psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
    185         psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
     181        psVectorClipFitPolynomial2D (map->x, stats, mask, 0xff, xObs, NULL, xRef, yRef);
     182        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, xObs->n);
     183
     184        psVectorClipFitPolynomial2D (map->y, stats, mask, 0xff, yObs, NULL, xRef, yRef);
     185        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues, yObs->n);
    186186    }
    187187
    188188    // loop over all chips, select the outliers, and replace the measured astrometry with the model
    189     // the measured transformation above must be applied to make the comparison, and also then applied to the 
     189    // the measured transformation above must be applied to make the comparison, and also then applied to the
    190190    // model transformation
    191191
    192192    if (DEBUG) {
    193         f = fopen ("corners.fit.dat", "w");
    194         for (int i = 0; i < xObs->n; i++) {
    195             psPlane obsCoord, refCoord;
    196             refCoord.x = xRef->data.F32[i];
    197             refCoord.y = yRef->data.F32[i];
    198             psPlaneTransformApply (&obsCoord, map, &refCoord);
    199             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);
    200         }
    201         fclose (f);
     193        f = fopen ("corners.fit.dat", "w");
     194        for (int i = 0; i < xObs->n; i++) {
     195            psPlane obsCoord, refCoord;
     196            refCoord.x = xRef->data.F32[i];
     197            refCoord.y = yRef->data.F32[i];
     198            psPlaneTransformApply (&obsCoord, map, &refCoord);
     199            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);
     200        }
     201        fclose (f);
    202202    }
    203203
    … …  
    211211
    212212    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    213         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    214         if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    215 
    216         // set the chip astrometry using the astrom file
    217         pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    218 
    219         // bad Astrometry test:  ref pixel or angle outside nominal
    220 
    221         psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
    222         psPlane obsCoord, refCoord, tmpCoord;
    223 
    224         // find location of 0,0 pixel in focal plane coords for this chip
    225         psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
    226 
    227         // find location of 0,0 pixel in focal plane coords for ref chip
    228         // apply the global field rotation and offset before comparing
    229         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
    230         psPlaneTransformApply (&refCoord, map, &tmpCoord);
    231    
    232         psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
    233         psPlane obsOffPt, refOffPt;
    234 
    235         // find location of 0,0 pixel in focal plane coords for this chip
    236         psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
    237 
    238         // find location of 0,0 pixel in focal plane coords for ref chip
    239         psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
    240         psPlaneTransformApply (&refOffPt, map, &tmpCoord);
    241    
    242         double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
    243         double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
    244 
    245         bool badAstrom = false;
    246         badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
    247         badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
    248         badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
    249 
    250         fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    251 
    252         // XXX for now, just use first readout
    253         pmCell *cell = obsChip->cells->data[0];
    254         pmReadout *readout = cell->readouts->data[0];
    255 
    256         // update the header (pull or create local view to entry on readout->analysis)
    257         psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    258         if (!updates) {
    259             updates = psMetadataAlloc ();
    260             psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    261             psFree (updates);
    262         }
    263 
    264         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
    265         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
    266         psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
    267 
    268         // for successful chips, save the measured offsets in the header
    269         if (!badAstrom) continue;
    270 
    271         // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
    272         // catch the failures relative to the model
    273         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
    274         continue;
    275 
    276         psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
    277                   view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
    278 
    279         psFree (obsChip->toFPA);
    280         psFree (obsChip->fromFPA);
    281 
    282         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    283         // XXX this only works if toTPA is at most a linear transformation
    284         psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
    285         for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
    286             for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
    287                 double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    288                 double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    289                 toFPA->x->coeff[i][j] = f1 + f2;
    290 
    291                 double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
    292                 double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
    293                 toFPA->y->coeff[i][j] = g1 + g2;
    294             }
    295         }
    296         toFPA->x->coeff[0][0] += map->x->coeff[0][0];
    297         toFPA->y->coeff[0][0] += map->y->coeff[0][0];
    298 
    299         psRegion *region = pmChipPixels (obsChip);
    300         obsChip->toFPA   = toFPA;
    301         obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
    302         psFree (region);
    303    
    304         // use the new position to re-try the match fit
    305         // select the raw objects for this readout
    306         psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    307         if (rawstars == NULL) { continue; }
    308 
    309         // select the raw objects for this readout
    310         psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    311         if (refstars == NULL) { continue; }
    312 
    313         // the absolute minimum number of stars is 4 (for order = 1)
    314         if ((rawstars->n < 4) || (refstars->n < 4)) {
    315             readout->data_exists = false;
    316             psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
    317                       rawstars->n, refstars->n);
    318             continue;
    319         }
    320 
    321         psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
    322 
    323         // XXX update the header with info to reflect the failure
    324         if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
    325             readout->data_exists = false;
    326             psLogMsg ("psastro", 3, "failed to find a solution\n");
    327             continue;
    328         }
    329 
    330         pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
    331     }
    332 
    333     psastroVisualPlotFixChips (input, xObs, yObs);
     213        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
     214        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     215
     216        // set the chip astrometry using the astrom file
     217        pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     218
     219        // bad Astrometry test:  ref pixel or angle outside nominal
     220
     221        psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
     222        psPlane obsCoord, refCoord, tmpCoord;
     223
     224        // find location of 0,0 pixel in focal plane coords for this chip
     225        psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);
     226
     227        // find location of 0,0 pixel in focal plane coords for ref chip
     228        // apply the global field rotation and offset before comparing
     229        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &refPixel);
     230        psPlaneTransformApply (&refCoord, map, &tmpCoord);
     231
     232        psPlane offPixel = {100.0, 0.0, 100.0, 0.0};
     233        psPlane obsOffPt, refOffPt;
     234
     235        // find location of 0,0 pixel in focal plane coords for this chip
     236        psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);
     237
     238        // find location of 0,0 pixel in focal plane coords for ref chip
     239        psPlaneTransformApply (&tmpCoord, refChip->toFPA, &offPixel);
     240        psPlaneTransformApply (&refOffPt, map, &tmpCoord);
     241
     242        double obsAngle = PM_DEG_RAD*atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
     243        double refAngle = PM_DEG_RAD*atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);
     244
     245        bool badAstrom = false;
     246        badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
     247        badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
     248        badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
     249
     250        fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     251
     252        // XXX for now, just use first readout
     253        pmCell *cell = obsChip->cells->data[0];
     254        pmReadout *readout = cell->readouts->data[0];
     255
     256        // update the header (pull or create local view to entry on readout->analysis)
     257        psMetadata *updates = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
     258        if (!updates) {
     259            updates = psMetadataAlloc ();
     260            psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     261            psFree (updates);
     262        }
     263
     264        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DX", PS_META_REPLACE, "chip x offset wrt model", obsCoord.x - refCoord.x);
     265        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
     266        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
     267
     268        // for successful chips, save the measured offsets in the header
     269        if (!badAstrom) continue;
     270
     271        // XXX for now, let's just fail on the bad chips.  In the future, let's try to recover, but we still need to
     272        // catch the failures relative to the model
     273        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO", PS_META_REPLACE, "number of astrometry stars", 0);
     274        continue;
     275
     276        psLogMsg ("psastro", PS_LOG_INFO, "fixing chip %d, angle: %f, pixel: %f,%f\n",
     277                  view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     278
     279        psFree (obsChip->toFPA);
     280        psFree (obsChip->fromFPA);
     281
     282        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     283        // XXX this only works if toTPA is at most a linear transformation
     284        psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
     285        for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
     286            for (int j = 0; j <= refChip->toFPA->x->nY; j++) {
     287                double f1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->x->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     288                double f2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->x->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     289                toFPA->x->coeff[i][j] = f1 + f2;
     290
     291                double g1 = refChip->toFPA->x->coeffMask[i][j] ? 0.0 : map->y->coeff[1][0]*refChip->toFPA->x->coeff[i][j];
     292                double g2 = refChip->toFPA->y->coeffMask[i][j] ? 0.0 : map->y->coeff[0][1]*refChip->toFPA->y->coeff[i][j];
     293                toFPA->y->coeff[i][j] = g1 + g2;
     294            }
     295        }
     296        toFPA->x->coeff[0][0] += map->x->coeff[0][0];
     297        toFPA->y->coeff[0][0] += map->y->coeff[0][0];
     298
     299        psRegion *region = pmChipPixels (obsChip);
     300        obsChip->toFPA   = toFPA;
     301        obsChip->fromFPA = psPlaneTransformInvert(NULL, obsChip->toFPA, *region, 50);
     302        psFree (region);
     303
     304        // use the new position to re-try the match fit
     305        // select the raw objects for this readout
     306        psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     307        if (rawstars == NULL) { continue; }
     308
     309        // select the raw objects for this readout
     310        psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     311        if (refstars == NULL) { continue; }
     312
     313        // the absolute minimum number of stars is 4 (for order = 1)
     314        if ((rawstars->n < 4) || (refstars->n < 4)) {
     315            readout->data_exists = false;
     316            psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
     317                      rawstars->n, refstars->n);
     318            continue;
     319        }
     320
     321        psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
     322
     323        // XXX update the header with info to reflect the failure
     324        if (!psastroOneChipFit (input->fpa, obsChip, refstars, rawstars, recipe, updates)) {
     325            readout->data_exists = false;
     326            psLogMsg ("psastro", 3, "failed to find a solution\n");
     327            continue;
     328        }
     329
     330        pmAstromWriteWCS (updates, input->fpa, obsChip, NONLIN_TOL);
     331    }
     332
     333    pmAstromVisualPlotFixChips (input, xObs, yObs);
    334334    psFree (xObs);
    335335    psFree (yObs);

  • trunk/psastro/src/psastroLoadRefstars.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    150150    }
    151151
    152     psastroVisualPlotRefStars (refstars, recipe);
     152    pmAstromVisualPlotRefStars (refstars, recipe);
    153153
    154154    if (psTraceGetLevel("psastro.plot") > 0) {
    … …  
    251251    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    252252    if (!input) {
    253         psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
    254         photcode = psStringCopy ("NONE");
    255         return photcode;
     253        psLogMsg ("psastro", PS_LOG_DETAIL, "no supplied reference header data");
     254        photcode = psStringCopy ("NONE");
     255        return photcode;
    256256    }
    257257    assert (input->fpa);
    … …  
    259259    *maxRho = psMetadataLookupF32(&status, recipe, "DVO.GETSTAR.MAX.RHO");
    260260    if (!status) {
    261         psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
    262         return NULL;
     261        psError(PSASTRO_ERR_CONFIG, false, "DVO.GETSTAR.MAX.RHO missing from recipe");
     262        return NULL;
    263263    }
    264264
    … …  
    278278    if (!status) ESCAPE ("missing DVO.GETSTAR.MIN.MAG.INST");
    279279
    280     // PHOTCODE.DATA is a multi of metadata items 
     280    // PHOTCODE.DATA is a multi of metadata items
    281281    psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false);
    282282
    283283    psMetadataItem *refItem = NULL;
    284284    while ((refItem = psListGetAndIncrement (iter))) {
    285         if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
    286    
    287         char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
    288         if (!status) {
    289             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    290             continue;
    291         }
    292 
    293         // does this entry match the current filter?
    294         if (strcmp (refFilter, filter)) continue;
    295 
    296         psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
    297 
    298         float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
    299         if (!status) {
    300             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    301             continue;
    302         }
    303         photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
    304         if (!status) {
    305             psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
    306             continue;
    307         }
    308 
    309         // convert the minInst to a calibrated minimum magnitude
    310         *minMag = minInst + 2.5*log10(exptime) + zeropt;
    311 
    312         psFree (iter);
    313         return photcode;
    314     }   
     285        if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
     286
     287        char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
     288        if (!status) {
     289            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     290            continue;
     291        }
     292
     293        // does this entry match the current filter?
     294        if (strcmp (refFilter, filter)) continue;
     295
     296        psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);
     297
     298        float zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
     299        if (!status) {
     300            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     301            continue;
     302        }
     303        photcode = psMetadataLookupStr (&status, refItem->data.md, "PHOTCODE");
     304        if (!status) {
     305            psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
     306            continue;
     307        }
     308
     309        // convert the minInst to a calibrated minimum magnitude
     310        *minMag = minInst + 2.5*log10(exptime) + zeropt;
     311
     312        psFree (iter);
     313        return photcode;
     314    }
    315315    psFree (iter);
    316316
    … …  
    318318    photcode = psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR.PHOTCODE");
    319319    PS_ASSERT (photcode, NULL);
    320        
     320
    321321    // give up and use fixed value
    322322    *minMag = psMetadataLookupF32(NULL, recipe, "DVO.GETSTAR.MIN.MAG");

  • trunk/psastro/src/psastroLuminosityFunction.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    141141    lumFunc->sPeak = sPeak;
    142142
    143     psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
     143    pmAstromVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
    144144
    145145    psFree (lnMag);

  • trunk/psastro/src/psastroMosaicOneChip.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    160160
    161161    //plot results
    162     psastroVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
     162    pmAstromVisualPlotMosaicOneChip(rawstars, refstars, match, recipe);
    163163
    164164    psFree (fitStats);

  • trunk/psastro/src/psastroMosaicSetMatch.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    6969                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    7070
    71                 psastroVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
     71                pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);
    7272
    7373                // XXX drop the old one

  • trunk/psastro/src/psastroOneChipFit.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    2323
    2424    // default value for match/fit : radius is in pixels
    25     REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
     25    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32);
    2626
    2727    // run the match/fit sequence NITER times
    28     REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
     28    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
    2929
    3030    // correct radius to FP units (physical pixel scale in microns per pixel)
    31     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     31    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
    3232    RADIUS *= pixelScale;
    3333
    … …  
    4444
    4545    for (int iter = 0; iter < nIter; iter++) {
    46        
    47         char name[128];
    4846
    49         sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
    50         float radius = psMetadataLookupF32 (&status, recipe, name);
    51         radius *= pixelScale;
    52         if (!status || (radius == 0.0)) {
    53             radius = RADIUS;
    54         }
     47        char name[128];
     48
     49        sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
     50        float radius = psMetadataLookupF32 (&status, recipe, name);
     51        radius *= pixelScale;
     52        if (!status || (radius == 0.0)) {
     53            radius = RADIUS;
     54        }
    5555
    5656
    57         // use small radius to match stars
    58         match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
    59         if (match == NULL) {
    60             psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
    61             return false;
    62         }
     57        // use small radius to match stars
     58        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
     59        if (match == NULL) {
     60            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
     61            return false;
     62        }
    6363
    64         // modify the order to correspond to the actual number of matched stars:
    65         int Ndof_min = 3;
    66         int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
    67         order = PS_MIN (order, order_max);
     64        // modify the order to correspond to the actual number of matched stars:
     65        int Ndof_min = 3;
     66        int order_max = 0.5*(3 + sqrt(4*match->n - 4*Ndof_min + 1));
     67        order = PS_MIN (order, order_max);
    6868
    69         // if ((match->n < 11) && (order >= 3)) order = 2;
    70         // if ((match->n <  7) && (order >= 2)) order = 1;
    71         // if ((match->n <  4) && (order >= 1)) order = 0;
     69        // if ((match->n < 11) && (order >= 3)) order = 2;
     70        // if ((match->n <  7) && (order >= 2)) order = 1;
     71        // if ((match->n <  4) && (order >= 1)) order = 0;
    7272
    73         if (order < 1) {
    74             psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
    75             psFree (match);
    76             return false;
    77         }
     73        if (order < 1) {
     74            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
     75            psFree (match);
     76            return false;
     77        }
    7878
    79         // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    80         psFree (chip->toFPA);
    81         chip->toFPA = psPlaneTransformAlloc (order, order);
    82         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    83             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    84                 if (i + j > order) {
    85                     chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    86                     chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    87                 }
    88             }
    89         }
     79        // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     80        psFree (chip->toFPA);
     81        chip->toFPA = psPlaneTransformAlloc (order, order);
     82        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     83            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     84                if (i + j > order) {
     85                    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     86                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     87                }
     88            }
     89        }
    9090
    91         // XXX allow statitic to be set by the user
    92         // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    93         fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    94         fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
    95         fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
     91        // XXX allow statitic to be set by the user
     92        // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     93        fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     94        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
     95        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
    9696
    97         // improved fit for astrometric terms
    98         results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    99         if (!results) {
    100             psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
    101             psFree (match);
    102             psFree (fitStats);
    103             return false;
    104         }
    105    
    106         // determine fromFPA transformation and apply new transformation to raw & ref stars
    107         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    108    
    109         // toSky converts from FPA & TPA units (microns) to sky units (radians)
    110         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    111         // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    112         float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
    113         int astNstar = results->yStats->clippedNvalues;
    114         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     97        // improved fit for astrometric terms
     98        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     99        if (!results) {
     100            psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
     101            psFree (match);
     102            psFree (fitStats);
     103            return false;
     104        }
    115105
    116         if (iter < nIter - 1) {
    117             psFree (fitStats);
    118             psFree (results);
    119             psFree (match);
    120         }
     106        // determine fromFPA transformation and apply new transformation to raw & ref stars
     107        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     108
     109        // toSky converts from FPA & TPA units (microns) to sky units (radians)
     110        float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
     111        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     112        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     113        int astNstar = results->yStats->clippedNvalues;
     114        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     115
     116        if (iter < nIter - 1) {
     117            psFree (fitStats);
     118            psFree (results);
     119            psFree (match);
     120        }
    121121    }
    122122
    … …  
    139139    if (astError > maxError) {
    140140        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    141         validSolution = false;
     141        validSolution = false;
    142142    }
    143143    if (astNstar < minNstar) {
    144144        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    145         validSolution = false;
     145        validSolution = false;
    146146    }
    147147
    … …  
    150150    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    151151    if (validSolution) {
    152         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    153         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     152        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     153        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    154154    } else {
    155         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    156         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     155        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     156        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    157157    }
    158158    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
    … …  
    160160    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
    161161    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    162    
     162
    163163    // XXX check if we correctly applied the new transformation:
    164164    if (psTraceGetLevel("psastro.dump") > 0) {
    165         psastroDumpRawstars (rawstars, fpa, chip);
    166         psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
    167         psastroDumpStars (refstars, "refstars.cal.dat");
     165        psastroDumpRawstars (rawstars, fpa, chip);
     166        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
     167        psastroDumpStars (refstars, "refstars.cal.dat");
    168168    }
    169169
    170     psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);
     170    pmAstromVisualPlotOneChipFit (rawstars, refstars, match, recipe);
    171171
    172172    if (psTraceGetLevel("psastro.plot") > 0) {
    173         psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     173        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
    174174    }
    175175

  • trunk/psastro/src/psastroRemoveClumps.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    6969    psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit);
    7070
    71     psastroVisualPlotRemoveClumps (input, count, scale, limit);
     71    pmAstromVisualPlotRemoveClumps (input, count, scale, limit);
    7272
    7373    // find and exclude objects in bad pixels

  • trunk/psastro/src/psastroUtils.c

    r21409 r21422  
    66 *
    77 *  @author IfA
    8  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2009-02-07 02:03:34 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2009-02-09 21:25:34 $
    1010 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    5454            if (!chip->toFPA) { continue; }
    5555
    56             if (chip->cells->n == 0) { continue; }
    57             pmCell *cell = chip->cells->data[0];
     56            if (chip->cells->n == 0) { continue; }
     57            pmCell *cell = chip->cells->data[0];
    5858            if (!cell->process || !cell->file_exists) { continue; }
    5959
    60             if (cell->readouts->n == 0) { continue; }
    61             pmReadout *readout = cell->readouts->data[0];
    62             if (! readout->data_exists) { continue; }
     60            if (cell->readouts->n == 0) { continue; }
     61            pmReadout *readout = cell->readouts->data[0];
     62            if (! readout->data_exists) { continue; }
    6363
    6464            pixelScale1 = hypot (chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1]);
    … …  
    100100    psastroMosaicSetAstrom (fpa);
    101101    if (!useExternal) {
    102         psastroVisualPlotCommonScale (fpa, oldScale);
     102        pmAstromVisualPlotCommonScale (fpa, oldScale);
    103103    }
    104104    psFree (oldScale);
Note: See TracChangeset for help on using the changeset viewer.