IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39926 for trunk/psModules


Ignore:
Timestamp:
Jan 6, 2017, 11:30:10 AM (10 years ago)
Author:
eugene
Message:

merging changes from czw dev branch (compare with r39924)

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psModules/src/astrom/pmAstrometryDistortion.c

    r23989 r39926  
    335335
    336336    psFree (fpa->fromTPA);
    337     fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50);
     337    psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4);
     338    fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, *region, 50);
     339    psFree (myPT);
    338340    psFree (region);
    339341

  • trunk/psModules/src/astrom/pmAstrometryModel.c

    r36834 r39926  
    320320    // the PosZero is the offset between the reported and actual POSANGLE values
    321321    float PosZero = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.POS_ZERO");  /// XXX be consistent with degrees v radians
     322    // Indicate which direction the position angle goes
     323    int rotParity = psMetadataLookupS32 (&status, file->fpa->concepts, "FPA.ROT_PARITY");
    322324    char *refChip = psMetadataLookupStr (&status, file->fpa->concepts, "FPA.REF.CHIP");
    323325
    … …  
    349351
    350352    psMetadataAddF32(row,    PS_LIST_TAIL, "POS_ZERO", PS_META_REPLACE, "POSANGLE offset (degrees)", PosZero);
     353    psMetadataAddS32(row,    PS_LIST_TAIL, "ROT_PARITY", PS_META_REPLACE, "rotator parity", rotParity);
    351354    psMetadataAddStr(row,    PS_LIST_TAIL, "REF_CHIP", PS_META_REPLACE, "reference chip for model", refChip);
    352355
    … …  
    624627    TRANSFER (file->fpa->concepts, row, "BORE_P0");
    625628    TRANSFER (file->fpa->concepts, row, "POS_ZERO");
     629    TRANSFER (file->fpa->concepts, row, "ROT_PARITY");
    626630    TRANSFER (file->fpa->concepts, row, "REF_CHIP");
    627631
    … …  
    687691    char *refChip  = psMetadataLookupStr(&status, file->fpa->concepts, "REF_CHIP"); REQUIRE (status, "missing REF_CHIP");
    688692
     693    int rotatorParity = psMetadataLookupS32(&status, file->fpa->concepts, "ROT_PARITY"); REQUIRE (status, "missing ROT_PARITY");
     694   
    689695    // XXX we've swapped the sign and parity of POS (add to model)
    690696    // apply true posangle = -(POS - POS_ZERO)
    691697    psLogMsg ("psModules.astrom", 4, "Position Angle: %f, Model Position Angle Zero Point: %f\n", POS, PosZero);
    692     psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, (PosZero - POS));
     698    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, rotatorParity * (PosZero - POS));
    693699    psFree (file->fpa->fromTPA);
    694700    file->fpa->fromTPA = fromTPA;

  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r36834 r39926  
    225225    psArray *ref,
    226226    psArray *match,
    227     psStats *stats)
     227    psStats *stats,
     228    const psMetadata *config)
    228229{
    229230    PS_ASSERT_PTR_NON_NULL(map, NULL);
    … …  
    232233    PS_ASSERT_PTR_NON_NULL(match, NULL);
    233234    PS_ASSERT_PTR_NON_NULL(stats, NULL);
     235    PS_ASSERT_PTR_NON_NULL(config, NULL);
     236
     237    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
     238    bool status;
     239    double photomWindowSigma  = psMetadataLookupF32 (&status, config, "PSASTRO.PHOTOM.WINDOW.SIGMA");
     240    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma) < 0.01));
     241    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
    234242
    235243    // reassign values for clip fit
    … …  
    252260        y->data.F32[i] = refStar->FP->y;
    253261
    254         wt->data.F32[i] = 1.0;
     262        // wt is used as an error (sqrt(variance)) in the fit. the 1.01 prevents the weight from going to 0.0 for perfect matches
     263        wt->data.F32[i] = 1.01 - exp(photomWindowFactor*PS_SQR(refStar->magCal - rawStar->magCal));
    255264    }
    256265
    … …  
    657666    obj->Color= 0;
    658667    obj->dMag = 0;
     668    obj->magCal = 0;
    659669
    660670    return (obj);
    … …  
    685695    obj->Color =  old->Color;
    686696    obj->dMag  =  old->dMag;
     697    obj->magCal =  old->magCal;
    687698
    688699    return(obj);
    … …  
    774785    double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET");
    775786
    776     // sampling scale of the grid
    777     double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
     787     // sampling scale of the grid
     788     double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
     789
     790    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
     791    double photomWindowSigma  = psMetadataLookupF32 (&status, config, "PSASTRO.PHOTOM.WINDOW.SIGMA");
     792    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma) < 0.01));
     793    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
    778794
    779795    // set the static scaling factors
    … …  
    816832            }
    817833
     834            // XXX should I make the scale factor in front a recipe value?
     835            int Npts = 10 * exp(photomWindowFactor*PS_SQR(ob1->magCal - ob2->magCal));
     836
    818837            // accumulate bin stats
    819             NP[iY][iX] ++;
    820             DX[iY][iX] += dX;
    821             DY[iY][iX] += dY;
    822             D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
     838            NP[iY][iX] += Npts;
     839            DX[iY][iX] += dX*Npts;
     840            DY[iY][iX] += dY*Npts;
     841            D2[iY][iX] += PS_SQR(dX*Npts) + PS_SQR(dY*Npts);
    823842        }
    824843    }
    … …  
    10531072    double tweakNsigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.NSIGMA");
    10541073
     1074    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
     1075    double photomWindowSigma  = psMetadataLookupF32 (&status, recipe, "PSASTRO.PHOTOM.WINDOW.SIGMA");
     1076    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma) < 0.01));
     1077    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
     1078
    10551079    nBin = 2*tweakRange / tweakScale;
    10561080    psVector *xHist = psVectorAlloc (nBin, PS_TYPE_F32);
    … …  
    10791103                continue;
    10801104
    1081             xHist->data.F32[xBin] += 1.0;
    1082             yHist->data.F32[yBin] += 1.0;
     1105            // XXX should I make the scale factor in front a recipe value?
     1106            int Npts = 10 * exp(photomWindowFactor*PS_SQR(ob1->magCal - ob2->magCal));  // sigma = 0.22 mag
     1107
     1108            xHist->data.F32[xBin] += Npts;
     1109            yHist->data.F32[yBin] += Npts;
    10831110        }
    10841111    }

  • trunk/psModules/src/astrom/pmAstrometryObjects.h

    r36834 r39926  
    4343    float dMag;                         ///< error on object magnitude
    4444    float SBinst;                       ///< surface brightness, used for Koppenhoefer correction
     45    float magCal;                       ///< object calibrated magnitude in extracted filter
    4546}
    4647pmAstromObj;
    … …  
    347348    psArray *ref,
    348349    psArray *match,
    349     psStats *stats
     350    psStats *stats,
     351    const psMetadata *config
    350352);
    351353

  • trunk/psModules/src/astrom/pmAstrometryWCS.c

    r36492 r39926  
    4343    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
    4444    if (!wcs) {
    45         return false;
     45        return false;
    4646    }
    4747
    … …  
    7070    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
    7171    if (!wcs) {
    72         return false;
     72        return false;
    7373    }
    7474
    … …  
    8686    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
    8787    if (!wcs) {
    88         psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from header");
    89         return false;
     88        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from header");
     89        return false;
    9090    }
    9191
    … …  
    9696
    9797    if (!status1 || !status2) {
    98         Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1");
    99         Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2");
     98        Nx = psMetadataLookupS32 (&status1, header, "IMNAXIS1");
     99        Ny = psMetadataLookupS32 (&status2, header, "IMNAXIS2");
    100100    }
    101101
    102102    if (!status1 || !status2) {
    103         Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
    104         Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
     103        Nx = psMetadataLookupS32 (&status1, header, "ZNAXIS1");
     104        Ny = psMetadataLookupS32 (&status2, header, "ZNAXIS2");
    105105    }
    106106
    107107    if (!status1 || !status2) {
    108         psFree (wcs);
    109         psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
    110         return false;
     108        psFree (wcs);
     109        psError(PS_ERR_UNKNOWN, false, "missing required FPA size in header");
     110        return false;
    111111    }
    112112
    … …  
    123123    pmAstromWCS *wcs = pmAstromWCSBilevelChipFromFPA (chip, tol);
    124124    if (!wcs) {
    125         psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
    126         return false;
     125        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
     126        return false;
    127127    }
    128128
    … …  
    139139    pmAstromWCS *wcs = pmAstromWCSBilevelMosaicFromFPA (fpa, tol);
    140140    if (!wcs) {
    141         psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
    142         return false;
     141        psError(PS_ERR_UNKNOWN, false, "failure to determine WCS terms from fpa");
     142        return false;
    143143    }
    144144
    … …  
    163163
    164164    if (chip == NULL)
    165         return false;
     165        return false;
    166166    if (sky == NULL)
    167         return false;
     167        return false;
    168168    if (wcs == NULL)
    169         return false;
     169        return false;
    170170
    171171    psPlane *Chip = psPlaneAlloc();
    … …  
    188188
    189189    if (chip == NULL)
    190         return false;
     190        return false;
    191191    if (sky == NULL)
    192         return false;
     192        return false;
    193193    if (wcs == NULL)
    194         return false;
     194        return false;
    195195
    196196    psError(PS_ERR_UNKNOWN, true, "not yet implemented: needs to invert the transformation");
    … …  
    223223    char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
    224224    if (!status) {
    225         psLogMsg ("psastro", 5, "warning: no WCS metadata in header\n");
    226         return NULL;
     225        psLogMsg ("psastro", 5, "warning: no WCS metadata in header\n");
     226        return NULL;
    227227    }
    228228
    … …  
    232232    type = psProjectTypeFromString (ctype);
    233233    if (type == PS_PROJ_NTYPE) {
    234         psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
    235         return NULL;
     234        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);
     235        return NULL;
    236236    }
    237237
    … …  
    243243
    244244    if (cdKeys && pcKeys) {
    245         // XXX make this an option
    246         psLogMsg ("psastro", 5, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
     245        // XXX make this an option
     246        psLogMsg ("psastro", 5, "warning: both CDi_j and PC00i00j defined in headers, using PC00i00j terms\n");
    247247    }
    248248    if (!cdKeys && !pcKeys) {
    249         psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms");
    250         // XXX we could default here to RA, DEC, ROTANGLE
    251         return NULL;
     249        psError(PS_ERR_UNKNOWN, true, "missing both CDi_j and PC00i00j WCS terms");
     250        // XXX we could default here to RA, DEC, ROTANGLE
     251        return NULL;
    252252    }
    253253    if (isPoly) {
    254         if (!pcKeys) {
    255             psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
    256             return NULL;
    257         }
    258         if (fitOrder == 0)
    259             fitOrder = 1;
    260         if ((fitOrder > 3) || (fitOrder < 1)) {
    261             psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
    262             return NULL;
    263         }
     254        if (!pcKeys) {
     255            psError(PS_ERR_UNKNOWN, true, "polynomial terms defined, but missing PC00i00j WCS terms");
     256            return NULL;
     257        }
     258        if (fitOrder == 0)
     259            fitOrder = 1;
     260        if ((fitOrder > 3) || (fitOrder < 1)) {
     261            psError(PS_ERR_UNKNOWN, true, "NPLYTERM value undefined: %d", fitOrder);
     262            return NULL;
     263        }
    264264    } else {
    265         fitOrder = 1;
     265        fitOrder = 1;
    266266    }
    267267
    … …  
    277277    wcs->crpix2 = psMetadataLookupF64 (&status, header, "CRPIX2");
    278278    wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
     279
     280    // XXX if type == ZPN, look for PV2_%d elements:
     281    if (type == PS_PROJ_ZPN) {
     282        psVector *maxRadial = psVectorAlloc (21, PS_TYPE_F64);
     283        for (int i = 0; i <= 20; i++) {
     284            char name[64];
     285            snprintf (name, 64, "PV2_%d", i);
     286
     287            maxRadial->data.F64[i] = 0.0;
     288            double value = psMetadataLookupF64 (&status, header, name);
     289
     290            if (status) {
     291                maxRadial->data.F64[i] = value;
     292                maxRadial->n = i;
     293            }
     294
     295            // PV2_1 is implicit if not present
     296            if ((i == 1) && !status) {
     297                maxRadial->data.F64[i] = 1.0;
     298                continue;
     299            }
     300        }
     301        maxRadial->n ++;
     302        wcs->toSky->radial = maxRadial;
     303    }
    279304
    280305    // These aren't needed but having them empty is disconcerting
    … …  
    289314    // test the CDELTi varient
    290315    if (pcKeys) {
    291         wcs->wcsCDkeys = 0;
    292         wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
    293         wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
    294 
    295         // test the CROTAi varient:
    296         // XXX double check lambda..
    297         double rotate = psMetadataLookupF64 (&status, header, "CROTA2");
    298         if (status) {
    299             wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
    300             wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
    301             wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
    302             wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
    303             return wcs;
    304         }
    305 
    306         // FITS WCS PCi,j has units of unity
    307         // wcs->trans has units of degrees/pixel
    308         wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC001001"); // == PC1_1
    309         wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC001002"); // == PC1_2
    310         wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC002001"); // == PC2_1
    311         wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC002002"); // == PC2_2
    312 
    313         if (isPoly) {
    314             // Elixir-style polynomial terms
    315             // XXX currently, Elixir/DVO cannot accept mixed orders
    316             for (int i = 0; i <= fitOrder; i++) {
    317                 for (int j = 0; j <= fitOrder; j++) {
    318                     if (i + j < 2)
    319                         continue;
    320                     if (i + j > fitOrder) {
    321                         wcs->trans->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    322                         wcs->trans->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    323                         continue;
    324                     }
    325                     sprintf (name, "PCA1X%1dY%1d", i, j);
    326                     wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
    327                     sprintf (name, "PCA2X%1dY%1d", i, j);
    328                     wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
    329                 }
    330             }
    331         }
    332         return wcs;
     316        wcs->wcsCDkeys = 0;
     317        wcs->cdelt1 = psMetadataLookupF64 (&status, header, "CDELT1");
     318        wcs->cdelt2 = psMetadataLookupF64 (&status, header, "CDELT2");
     319
     320        // test the CROTAi varient:
     321        // XXX double check lambda..
     322        double rotate = psMetadataLookupF64 (&status, header, "CROTA2");
     323        if (status) {
     324            wcs->trans->x->coeff[1][0] = +wcs->cdelt1 * cos(rotate*PM_RAD_DEG); // == PC1_1
     325            wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
     326            wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
     327            wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
     328            return wcs;
     329        }
     330
     331        // FITS WCS PCi,j has units of unity
     332        // wcs->trans has units of degrees/pixel
     333        wcs->trans->x->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC001001"); // == PC1_1
     334        wcs->trans->x->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC001002"); // == PC1_2
     335        wcs->trans->y->coeff[1][0] = wcs->cdelt1 * psMetadataLookupF64 (&status, header, "PC002001"); // == PC2_1
     336        wcs->trans->y->coeff[0][1] = wcs->cdelt2 * psMetadataLookupF64 (&status, header, "PC002002"); // == PC2_2
     337
     338        if (isPoly) {
     339            // Elixir-style polynomial terms
     340            // XXX currently, Elixir/DVO cannot accept mixed orders
     341            for (int i = 0; i <= fitOrder; i++) {
     342                for (int j = 0; j <= fitOrder; j++) {
     343                    if (i + j < 2)
     344                        continue;
     345                    if (i + j > fitOrder) {
     346                        wcs->trans->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     347                        wcs->trans->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     348                        continue;
     349                    }
     350                    sprintf (name, "PCA1X%1dY%1d", i, j);
     351                    wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
     352                    sprintf (name, "PCA2X%1dY%1d", i, j);
     353                    wcs->trans->y->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF64 (&status, header, name);
     354                }
     355            }
     356        }
     357        return wcs;
    333358    }
    334359
    335360    // test the CDi_j varient
    336361    if (cdKeys) {
    337         wcs->wcsCDkeys = 1;
    338 
    339         wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
    340         wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
    341         wcs->trans->y->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD2_1"); // == PC2_1
    342         wcs->trans->y->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD2_2"); // == PC2_2
    343         wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
    344         wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
    345         return wcs;
     362        wcs->wcsCDkeys = 1;
     363
     364        wcs->trans->x->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD1_1"); // == PC1_1
     365        wcs->trans->x->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD1_2"); // == PC1_2
     366        wcs->trans->y->coeff[1][0] = psMetadataLookupF64 (&status, header, "CD2_1"); // == PC2_1
     367        wcs->trans->y->coeff[0][1] = psMetadataLookupF64 (&status, header, "CD2_2"); // == PC2_2
     368        wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
     369        wcs->cdelt2 = hypot (wcs->trans->y->coeff[1][0], wcs->trans->y->coeff[0][1]);
     370        return wcs;
    346371    }
    347372    psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");
    … …  
    376401    psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", wcs->crpix2);
    377402
     403    if (wcs->toSky->type == PS_PROJ_ZPN) {
     404        psAssert (wcs->toSky->radial, "missing radial vector");
     405        for (int i = 0; i < wcs->toSky->radial->n; i++) {
     406            if (wcs->toSky->radial->data.F64[i] == 0.0) continue;
     407            snprintf (name, 16, "PV2_%d", i);
     408            psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->toSky->radial->data.F64[i]);
     409        }
     410    }
     411
    378412    // XXX make it optional to write out CDi_j terms, or other versions
    379413    // apply CDELT1,2 (degrees / pixel) to yield PCi,j terms of order unity
    380414    if (!wcs->wcsCDkeys) {
    381415
    382       double cdelt1 = wcs->cdelt1;
    383       double cdelt2 = wcs->cdelt2;
    384       psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
    385       psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
    386 
    387       // test the PC00i00j varient:
    388       psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
    389       psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
    390       psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
    391       psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
    392 
    393       // Elixir-style polynomial terms
    394       // XXX currently, Elixir/DVO cannot accept mixed orders
    395       // XXX need to respect the masks
    396       // XXX is wcs->cdelt1,2 always consistent?
    397       int fitOrder = wcs->trans->x->nX;
    398       if (fitOrder > 1) {
    399         for (int i = 0; i <= fitOrder; i++) {
    400           for (int j = 0; j <= fitOrder; j++) {
    401             if (i + j < 2)
    402               continue;
    403             if (i + j > fitOrder)
    404               continue;
    405             sprintf (name, "PCA1X%1dY%1d", i, j);
    406             psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    407             sprintf (name, "PCA2X%1dY%1d", i, j);
    408             psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
    409           }
    410         }
    411         psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
    412       }
    413 
    414       // remove any existing 'CDi_j style' wcs keywords
    415       if (psMetadataLookup(header, "CD1_1")) {
    416         psMetadataRemoveKey(header, "CD1_1");
    417         psMetadataRemoveKey(header, "CD1_2");
    418         psMetadataRemoveKey(header, "CD2_1");
    419         psMetadataRemoveKey(header, "CD2_2");
    420       }
    421 
    422       // Remove 'CDi_jX' WCS keywords
    423       psString cd11 = psStringCopy("CD1_1 ");
    424       psString cd12 = psStringCopy("CD1_2 ");
    425       psString cd21 = psStringCopy("CD2_1 ");
    426       psString cd22 = psStringCopy("CD2_2 ");
    427       for (char extra = 'A'; extra <= 'Z'; extra++) {
    428           cd11[strlen(cd11)-1] = extra;
    429           if (psMetadataLookup(header, cd11)) {
    430               cd12[strlen(cd12)-1] = extra;
    431               cd21[strlen(cd21)-1] = extra;
    432               cd22[strlen(cd22)-1] = extra;
    433               psMetadataRemoveKey(header, cd11);
    434               psMetadataRemoveKey(header, cd12);
    435               psMetadataRemoveKey(header, cd21);
    436               psMetadataRemoveKey(header, cd22);
    437           }
    438       }
    439       psFree(cd11);
    440       psFree(cd12);
    441       psFree(cd21);
    442       psFree(cd22);
     416        double cdelt1 = wcs->cdelt1;
     417        double cdelt2 = wcs->cdelt2;
     418        psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT1", PS_META_REPLACE, "", cdelt1);
     419        psMetadataAddF64 (header, PS_LIST_TAIL, "CDELT2", PS_META_REPLACE, "", cdelt2);
     420
     421        // test the PC00i00j varient:
     422        psMetadataAddF64 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0] / cdelt1); // == PC1_1
     423        psMetadataAddF64 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1] / cdelt2); // == PC1_2
     424        psMetadataAddF64 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0] / cdelt1); // == PC2_1
     425        psMetadataAddF64 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1] / cdelt2); // == PC2_2
     426
     427        // Elixir-style polynomial terms
     428        // XXX currently, Elixir/DVO cannot accept mixed orders
     429        // XXX need to respect the masks
     430        // XXX is wcs->cdelt1,2 always consistent?
     431        int fitOrder = wcs->trans->x->nX;
     432        if (fitOrder > 1) {
     433            for (int i = 0; i <= fitOrder; i++) {
     434                for (int j = 0; j <= fitOrder; j++) {
     435                    if (i + j < 2)
     436                        continue;
     437                    if (i + j > fitOrder)
     438                        continue;
     439                    sprintf (name, "PCA1X%1dY%1d", i, j);
     440                    psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->x->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
     441                    sprintf (name, "PCA2X%1dY%1d", i, j);
     442                    psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->trans->y->coeff[i][j] / pow(cdelt1, i) / pow(cdelt2, j));
     443                }
     444            }
     445            psMetadataAddS32 (header, PS_LIST_TAIL, "NPLYTERM", PS_META_REPLACE, "", fitOrder);
     446        }
     447
     448        // remove any existing 'CDi_j style' wcs keywords
     449        if (psMetadataLookup(header, "CD1_1")) {
     450            psMetadataRemoveKey(header, "CD1_1");
     451            psMetadataRemoveKey(header, "CD1_2");
     452            psMetadataRemoveKey(header, "CD2_1");
     453            psMetadataRemoveKey(header, "CD2_2");
     454        }
     455
     456        // Remove 'CDi_jX' WCS keywords
     457        psString cd11 = psStringCopy("CD1_1 ");
     458        psString cd12 = psStringCopy("CD1_2 ");
     459        psString cd21 = psStringCopy("CD2_1 ");
     460        psString cd22 = psStringCopy("CD2_2 ");
     461        for (char extra = 'A'; extra <= 'Z'; extra++) {
     462            cd11[strlen(cd11)-1] = extra;
     463            if (psMetadataLookup(header, cd11)) {
     464                cd12[strlen(cd12)-1] = extra;
     465                cd21[strlen(cd21)-1] = extra;
     466                cd22[strlen(cd22)-1] = extra;
     467                psMetadataRemoveKey(header, cd11);
     468                psMetadataRemoveKey(header, cd12);
     469                psMetadataRemoveKey(header, cd21);
     470                psMetadataRemoveKey(header, cd22);
     471            }
     472        }
     473        psFree(cd11);
     474        psFree(cd12);
     475        psFree(cd21);
     476        psFree(cd22);
    443477
    444478
    445479    } else {
    446480
    447       psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
    448       psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
    449       psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
    450       psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
    451 
    452       if (psMetadataLookup(header, "PC001001")) {
    453         psMetadataRemoveKey(header, "PC001001");
    454         psMetadataRemoveKey(header, "PC001002");
    455         psMetadataRemoveKey(header, "PC002001");
    456         psMetadataRemoveKey(header, "PC002002");
    457       }
     481        psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_1", PS_META_REPLACE, "", wcs->trans->x->coeff[1][0]);
     482        psMetadataAddF64 (header, PS_LIST_TAIL, "CD1_2", PS_META_REPLACE, "", wcs->trans->x->coeff[0][1]);
     483        psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_1", PS_META_REPLACE, "", wcs->trans->y->coeff[1][0]);
     484        psMetadataAddF64 (header, PS_LIST_TAIL, "CD2_2", PS_META_REPLACE, "", wcs->trans->y->coeff[0][1]);
     485
     486        if (psMetadataLookup(header, "PC001001")) {
     487            psMetadataRemoveKey(header, "PC001001");
     488            psMetadataRemoveKey(header, "PC001002");
     489            psMetadataRemoveKey(header, "PC002001");
     490            psMetadataRemoveKey(header, "PC002002");
     491        }
    458492    }
    459493
    … …  
    473507    // cdelt1,2 has units of degree/pixel
    474508    for (int i = 0; i <= toFPA->x->nX; i++) {
    475         for (int j = 0; j <= toFPA->x->nX; j++) {
    476             toFPA->x->coeff[i][j] *= pixelScale/wcs->cdelt1;
    477             toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
    478         }
     509        for (int j = 0; j <= toFPA->x->nX; j++) {
     510            toFPA->x->coeff[i][j] *= pixelScale/wcs->cdelt1;
     511            toFPA->y->coeff[i][j] *= pixelScale/wcs->cdelt2;
     512        }
    479513    }
    480514
    … …  
    485519    // projection from TPA (linear microns) to SKY (radians)
    486520    psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG*pdelt1, PM_RAD_DEG*pdelt2, wcs->toSky->type);
     521    toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
    487522
    488523    if (fpa->toSky == NULL) {
    489524        psFree(fpa->toTPA);
    490525        psFree(fpa->fromTPA);
    491         fpa->toTPA = psPlaneTransformIdentity (1);
    492         fpa->fromTPA = psPlaneTransformIdentity (1);
    493         fpa->toSky = toSky;
     526        fpa->toTPA = psPlaneTransformIdentity (1);
     527        fpa->fromTPA = psPlaneTransformIdentity (1);
     528        fpa->toSky = toSky;
    494529    } else {
    495530
    496         // this section allows the loaded chip to be included in an fpa structure in which
    497         // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have
    498         // already been defined).  we have to adjust to match the existing transformation.
    499 
    500         if (fpa->toTPA == NULL)
    501             psAbort("projection defined, tangent-plane not defined");
    502         if (fpa->fromTPA == NULL)
    503             psAbort("projection defined, tangent-plane not defined");
    504 
    505         // convert from pixels on this chip to pixels on reference chip
    506         // rX has units of refpixels / pixel
    507         double rX = toSky->Xs / fpa->toSky->Xs;
    508         double rY = toSky->Ys / fpa->toSky->Ys;
    509 
    510         for (int i = 0; i <= toFPA->x->nX; i++) {
    511             for (int j = 0; j <= toFPA->x->nY; j++) {
    512                 toFPA->x->coeff[i][j] *= rX;
    513                 toFPA->y->coeff[i][j] *= rY;
    514             }
    515         }
    516 
    517         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    518         // XXX this only works if toTPA is at most a linear transformation
    519         psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
    520         for (int i = 0; i <= toFPA->x->nX; i++) {
    521           for (int j = 0; j <= toFPA->x->nY; j++) {
    522             double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
    523             double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
    524             toFPAnew->x->coeff[i][j] = f1 + f2;
    525 
    526             double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
    527             double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
    528             toFPAnew->y->coeff[i][j] = g1 + g2;
    529           }
    530         }
    531         toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
    532         toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
    533 
    534         psFree (toFPA);
    535         toFPA = toFPAnew;
    536 
    537         // adjust reference pixel for new toSky reference coordinate
    538         // find the FPA coordinate of 0,0 for this chip.
    539         psPlane *fpOld = psPlaneAlloc();
    540         psPlane *fpNew = psPlaneAlloc();
    541         psPlane *tp = psPlaneAlloc();
    542         psSphere *sky = psSphereAlloc();
    543 
    544         sky->r = toSky->R;
    545         sky->d = toSky->D;
    546         psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    547         psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
    548 
    549         sky->r = fpa->toSky->R;
    550         sky->d = fpa->toSky->D;
    551         psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
    552         psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
    553 
    554         toFPA->x->coeff[0][0] -= fpNew->x - fpOld->x;
    555         toFPA->y->coeff[0][0] -= fpNew->y - fpOld->y;
    556 
    557         psFree (sky);
    558         psFree (tp);
    559         psFree (fpOld);
    560         psFree (fpNew);
    561 
    562         psFree (toSky);
     531        // this section allows the loaded chip to be included in an fpa structure in which
     532        // other chips have already been loaded (ie, the fpa->toTPA, fpa->toSky components have
     533        // already been defined).  we have to adjust to match the existing transformation.
     534
     535        if (fpa->toTPA == NULL)
     536            psAbort("projection defined, tangent-plane not defined");
     537        if (fpa->fromTPA == NULL)
     538            psAbort("projection defined, tangent-plane not defined");
     539
     540        // convert from pixels on this chip to pixels on reference chip
     541        // rX has units of refpixels / pixel
     542        double rX = toSky->Xs / fpa->toSky->Xs;
     543        double rY = toSky->Ys / fpa->toSky->Ys;
     544
     545        for (int i = 0; i <= toFPA->x->nX; i++) {
     546            for (int j = 0; j <= toFPA->x->nY; j++) {
     547                toFPA->x->coeff[i][j] *= rX;
     548                toFPA->y->coeff[i][j] *= rY;
     549            }
     550        }
     551
     552        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
     553        // XXX this only works if toTPA is at most a linear transformation
     554        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
     555        for (int i = 0; i <= toFPA->x->nX; i++) {
     556            for (int j = 0; j <= toFPA->x->nY; j++) {
     557                double f1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[1][0]*toFPA->x->coeff[i][j];
     558                double f2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->x->coeff[0][1]*toFPA->y->coeff[i][j];
     559                toFPAnew->x->coeff[i][j] = f1 + f2;
     560
     561                double g1 = toFPA->x->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[1][0]*toFPA->x->coeff[i][j];
     562                double g2 = toFPA->y->coeffMask[i][j] ? 0.0 : fpa->fromTPA->y->coeff[0][1]*toFPA->y->coeff[i][j];
     563                toFPAnew->y->coeff[i][j] = g1 + g2;
     564            }
     565        }
     566        toFPAnew->x->coeff[0][0] += fpa->fromTPA->x->coeff[0][0];
     567        toFPAnew->y->coeff[0][0] += fpa->fromTPA->y->coeff[0][0];
     568
     569        psFree (toFPA);
     570        toFPA = toFPAnew;
     571
     572        // adjust reference pixel for new toSky reference coordinate
     573        // find the FPA coordinate of 0,0 for this chip.
     574        psPlane *fpOld = psPlaneAlloc();
     575        psPlane *fpNew = psPlaneAlloc();
     576        psPlane *tp = psPlaneAlloc();
     577        psSphere *sky = psSphereAlloc();
     578
     579        sky->r = toSky->R;
     580        sky->d = toSky->D;
     581        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     582        psPlaneTransformApply (fpOld, fpa->fromTPA, tp);
     583
     584        sky->r = fpa->toSky->R;
     585        sky->d = fpa->toSky->D;
     586        psProject (tp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     587        psPlaneTransformApply (fpNew, fpa->fromTPA, tp);
     588
     589        toFPA->x->coeff[0][0] -= fpNew->x - fpOld->x;
     590        toFPA->y->coeff[0][0] -= fpNew->y - fpOld->y;
     591
     592        psFree (sky);
     593        psFree (tp);
     594        psFree (fpOld);
     595        psFree (fpNew);
     596
     597        psFree (toSky);
    563598    }
    564599
    … …  
    578613    // XXX if the inversion fails, we probably do not have a valid transform anyway
    579614    if (!chip->fromFPA) {
    580       psWarning ("failed to find a valid transformation");
    581       psFree (chip->toFPA);
    582       return false;
     615        psWarning ("failed to find a valid transformation");
     616        psFree (chip->toFPA);
     617        return false;
    583618    }
    584619
    585620    // this can take a very long time...
    586621    while (fpa->toSky->R < 0)
    587         fpa->toSky->R += 2.0*M_PI;
     622        fpa->toSky->R += 2.0*M_PI;
    588623    while (fpa->toSky->R > 2.0*M_PI)
    589         fpa->toSky->R -= 2.0*M_PI;
     624        fpa->toSky->R -= 2.0*M_PI;
    590625
    591626    fpa->wcsCDkeys = wcs->wcsCDkeys;
    592627
    593628    psTrace ("psastro", 5, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
    594              chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
    595              chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
    596              chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
     629             chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
     630             chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
     631             chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
    597632
    598633    psTrace ("psastro", 5, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
    599              chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
    600              chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
    601              chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
     634             chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
     635             chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
     636             chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
    602637
    603638    return true;
    … …  
    613648     */
    614649
    615     // create transformation with 0,0 reference pixel and units of microns/pixel
    616650    psFree (chip->toFPA);
    617     chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
     651    chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY);
     652
     653    // copy the toFPA x,y, transformations to the wcs version
     654    chip->toFPA->x = psPolynomial2DCopy (chip->toFPA->x, wcs->trans->x);
     655    chip->toFPA->y = psPolynomial2DCopy (chip->toFPA->y, wcs->trans->y);
     656
     657    // these need to be set based on crval1,2
     658    chip->toFPA->x->coeff[0][0] = wcs->crval1;
     659    chip->toFPA->y->coeff[0][0] = wcs->crval2;
    618660
    619661    // determine the inverse transformation: we need the chip pixels covered by this transform
    … …  
    633675    // cdelt1,2 has units of degrees/micron
    634676    fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->cdelt1*PM_RAD_DEG, wcs->cdelt2*PM_RAD_DEG, wcs->toSky->type);
     677    fpa->toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
    635678
    636679    // create transformation with 0,0 reference pixel
    … …  
    639682    // convert fpa->toTPA to units of unity (microns/micron)
    640683    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
    641         for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
    642             fpa->toTPA->x->coeff[i][j] /= wcs->cdelt1;
    643             fpa->toTPA->y->coeff[i][j] /= wcs->cdelt2;
    644         }
     684        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
     685            fpa->toTPA->x->coeff[i][j] /= wcs->cdelt1;
     686            fpa->toTPA->y->coeff[i][j] /= wcs->cdelt2;
     687        }
    645688    }
    646689
    … …  
    648691    // the region defines the FPA pixels covered by the tranformation
    649692    psFree (fpa->fromTPA);
    650     fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50);
     693    psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4);
     694    fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, region, 50);
     695    psFree (myPT);
    651696    return true;
    652697}
    … …  
    684729
    685730    for (int i = 0; i <= toTPA->x->nX; i++) {
    686         for (int j = 0; j <= toTPA->x->nY; j++) {
    687             double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
    688             double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
    689             toTPA->x->coeff[i][j] = f1 + f2;
    690 
    691             double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
    692             double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
    693             toTPA->y->coeff[i][j] = g1 + g2;
    694         }
     731        for (int j = 0; j <= toTPA->x->nY; j++) {
     732            double f1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[1][0]*chip->toFPA->x->coeff[i][j];
     733            double f2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->x->coeff[0][1]*chip->toFPA->y->coeff[i][j];
     734            toTPA->x->coeff[i][j] = f1 + f2;
     735
     736            double g1 = chip->toFPA->x->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[1][0]*chip->toFPA->x->coeff[i][j];
     737            double g2 = chip->toFPA->y->coeffMask[i][j] ? 0.0 : fpa->toTPA->y->coeff[0][1]*chip->toFPA->y->coeff[i][j];
     738            toTPA->y->coeff[i][j] = g1 + g2;
     739        }
    695740    }
    696741    toTPA->x->coeff[0][0] += fpa->toTPA->x->coeff[0][0];
    … …  
    701746    // convert projection from FPA to SKY into wcs projection (degrees to radians)
    702747    wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
     748    wcs->toSky->radial = psMemIncrRefCounter (fpa->toSky->radial);
     749
    703750    wcs->crval1 = fpa->toSky->R*PS_DEG_RAD;
    704751    wcs->crval2 = fpa->toSky->D*PS_DEG_RAD;
    … …  
    714761    psPlane *center = psPlaneTransformGetCenter (tpa1, tol);
    715762    if (!center) {
    716         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    717         psFree (toTPA);
    718         psFree (tpa1);
    719         psFree (wcs);
    720         return NULL;
     763        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
     764        psFree (toTPA);
     765        psFree (tpa1);
     766        psFree (wcs);
     767        return NULL;
    721768    }
    722769
    … …  
    751798    // convert wcs->trans to a matrix with units of degrees/pixel
    752799    for (int i = 0; i <= wcs->trans->x->nX; i++) {
    753         for (int j = 0; j <= wcs->trans->x->nY; j++) {
    754             wcs->trans->x->coeff[i][j] *= pdelt1;
    755             wcs->trans->y->coeff[i][j] *= pdelt2;
    756         }
     800        for (int j = 0; j <= wcs->trans->x->nY; j++) {
     801            wcs->trans->x->coeff[i][j] *= pdelt1;
     802            wcs->trans->y->coeff[i][j] *= pdelt2;
     803        }
    757804    }
    758805
    … …  
    773820*/
    774821
    775 // convert the chip-level toFPA to a wcs polynomial transformation
     822// convert the chip-level toFPA to a wcs polynomial transformation.  the pmAstromWCS
     823// structure represents a single layer transformation (e.g., RA-TAN, RA-WRP).  Here we are
     824// converting the chip-level to a WRP projection in the structure.  Later, this will be
     825// converted to the WCS keywords
     826
    776827pmAstromWCS *pmAstromWCSBilevelChipFromFPA (const pmChip *chip, double tol)
    777828{
    … …  
    784835    pmAstromWCS *wcs = pmAstromWCSAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
    785836
     837    // copy the toFPA x,y, transformations to the wcs version
     838    wcs->trans->x = psPolynomial2DCopy (wcs->trans->x, chip->toFPA->x);
     839    wcs->trans->y = psPolynomial2DCopy (wcs->trans->y, chip->toFPA->y);
     840
    786841    // Chip to FPA transformation is a Cartesian 'projection'
    787842    // reference pixel for FPA is 0.0, 0.0
    788843    wcs->toSky = psProjectionAlloc (0.0, 0.0, 1.0, 1.0, PS_PROJ_WRP);
    789     wcs->crval1 = 0.0;
    790     wcs->crval2 = 0.0;
    791 
    792     // given transformation, solve for coordinates which yields output coordinates of 0,0
    793     psPlane *center = psPlaneTransformGetCenter (chip->toFPA, tol);
    794     if (!center) {
    795         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    796         psFree (wcs);
    797         return NULL;
    798     }
    799 
    800     // adjust wcs transform to use center as reference coordinate
    801     // resulting transformation has units of microns/pixel
    802     psPlaneTransformSetCenter (wcs->trans, chip->toFPA, center->x, center->y);
    803 
    804     // calculated center is crpix1,2
    805     wcs->crpix1 = center->x;
    806     wcs->crpix2 = center->y;
    807     psFree (center);
    808 
     844
     845    // reference pixel (CRPIX1,2) is (0.0, 0.0):
     846    wcs->crpix1 = 0.0;
     847    wcs->crpix2 = 0.0;
     848
     849    // we need to set CRVAL1,2 for the 0,0 pixel:
     850    wcs->crval1 = psPolynomial2DEval (chip->toFPA->x, 0.0, 0.0);
     851    wcs->crval2 = psPolynomial2DEval (chip->toFPA->y, 0.0, 0.0);
     852
     853    wcs->toSky->R = wcs->crval1*PM_RAD_DEG;
     854    wcs->toSky->D = wcs->crval2*PM_RAD_DEG;
     855
     856    // these need to be set to 0.0 since they have been moved to crpix1,crpix2
     857    wcs->trans->x->coeff[0][0] = 0.0;
     858    wcs->trans->y->coeff[0][0] = 0.0;
     859   
    809860    // output coordinates are in microns : CDELT1,2 has units of microns/pixel
    810861    wcs->cdelt1 = hypot (wcs->trans->x->coeff[1][0], wcs->trans->x->coeff[0][1]);
    … …  
    834885    psPlane *center = psPlaneTransformGetCenter (fpa->toTPA, tol);
    835886    if (!center) {
    836         psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
    837         psFree (wcs);
    838         return NULL;
     887        psError(PS_ERR_UNKNOWN, false, "Unable to solve for TPA center.");
     888        psFree (wcs);
     889        return NULL;
    839890    }
    840891
    … …  
    854905    // convert wcs->trans to units of degree/micron
    855906    for (int i = 0; i <= wcs->trans->x->nX; i++) {
    856         for (int j = 0; j <= wcs->trans->x->nY; j++) {
    857             wcs->trans->x->coeff[i][j] *= pdelt1;
    858             wcs->trans->y->coeff[i][j] *= pdelt2;
    859         }
     907        for (int j = 0; j <= wcs->trans->x->nY; j++) {
     908            wcs->trans->x->coeff[i][j] *= pdelt1;
     909            wcs->trans->y->coeff[i][j] *= pdelt2;
     910        }
    860911    }
    861912
    … …  
    880931    int k=0;
    881932    for (int j=0; j<nSamples; j++) {
    882         double y = bounds->y0 + (j * deltaY / nSamples);
    883         for (int i=0; i<nSamples; i++) {
    884             psPlane *s = psPlaneAlloc();
    885             s->x = bounds->x0 + (i * deltaX / nSamples);
    886             s->y = y;
    887             psArraySet(src, k, s);
    888             psPlane *d = psPlaneTransformApply(NULL, trans, s);
    889             psArraySet(dst, k, d);
    890             psFree(s);  // drop our refs to s and d
    891             psFree(d);
    892             ++k;
    893         }
     933        double y = bounds->y0 + (j * deltaY / nSamples);
     934        for (int i=0; i<nSamples; i++) {
     935            psPlane *s = psPlaneAlloc();
     936            s->x = bounds->x0 + (i * deltaX / nSamples);
     937            s->y = y;
     938            psArraySet(src, k, s);
     939            psPlane *d = psPlaneTransformApply(NULL, trans, s);
     940            psArraySet(dst, k, d);
     941            psFree(s);  // drop our refs to s and d
     942            psFree(d);
     943            ++k;
     944        }
    894945    }
    895946
    … …  
    897948
    898949    if (!psPlaneTransformFit(newTrans, src, dst, 0, 0)) {
    899         psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
    900         return NULL;
     950        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
     951        return NULL;
    901952    }
    902953
    … …  
    907958    printf("   i     chip_x  tpa_x     tpa_x_fit     dx         chip_y    tpa_y     tpa_y_fit     dy     dx > 0.5 || dy > 0.5\n");
    908959    for (int i=0; i<psArrayLength(dst); i++) {
    909         psPlane *d = (psPlane *) psArrayGet(dst, i);
    910         psPlane *s = (psPlane *) psArrayGet(src, i);
    911 
    912         new = psPlaneTransformApply(new, newTrans, s);
    913 
    914         double xerr = new->x - d->x;
    915         double yerr = new->y - d->y;
    916         bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5);
    917         printf("%4d %9.2f %9.2f %9.2f %9.4f     %9.2f %9.2f %9.2f %9.4f   %s\n"
    918         , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
     960        psPlane *d = (psPlane *) psArrayGet(dst, i);
     961        psPlane *s = (psPlane *) psArrayGet(src, i);
     962
     963        new = psPlaneTransformApply(new, newTrans, s);
     964
     965        double xerr = new->x - d->x;
     966        double yerr = new->y - d->y;
     967        bool bigerr = (fabs(xerr) > .5) || (fabs(yerr) > .5);
     968        printf("%4d %9.2f %9.2f %9.2f %9.4f     %9.2f %9.2f %9.2f %9.4f   %s\n"
     969               , i, s->x, new->x, d->x, xerr, s->y, new->y, d->y, yerr, bigerr ? "BIGERR" : "");
    919970    }
    920971    psFree(new);
    … …  
    934985
    935986    if (outFPA == NULL) {
    936         outFPA = inFPA;
     987        outFPA = inFPA;
    937988    }
    938989    if (outChip == NULL) {
    939         outChip = inChip;
     990        outChip = inChip;
    940991    }
    941992    if (outputBounds == NULL) {
    942         outputBounds = pmChipPixels(outChip);
     993        outputBounds = pmChipPixels(outChip);
    943994    }
    944995
    … …  
    946997    psPlaneTransform *chipToTPA = psPlaneTransformCombine(NULL, inChip->toFPA, inFPA->toTPA, *outputBounds, 50);
    947998    if (!chipToTPA) {
    948         psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
    949         return false;
     999        psError(PS_ERR_UNKNOWN, false, "failed to create chipToTPA");
     1000        return false;
    9501001    }
    9511002
    … …  
    9541005    psFree(chipToTPA);
    9551006    if (!chipToFPA) {
    956         psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
    957         return false;
     1007        psError(PS_ERR_UNKNOWN, false, "linear fit of chip to TPA transform failed");
     1008        return false;
    9581009    }
    9591010
    … …  
    9611012    psPlaneTransform *outToFPA;
    9621013    if (offset_x != 0. && offset_y != 0.) {
    963         outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y);
    964         psFree(chipToFPA);
     1014        outToFPA = psPlaneTransformSetCenter(NULL, chipToFPA, offset_x, offset_y);
     1015        psFree(chipToFPA);
    9651016    } else {
    966         outToFPA = chipToFPA;
     1017        outToFPA = chipToFPA;
    9671018    }
    9681019
    9691020    psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50);
    9701021    if (!outFromFPA) {
    971         psFree(outToFPA);
    972         psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
    973         return false;
     1022        psFree(outToFPA);
     1023        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
     1024        return false;
    9741025    }
    9751026
    … …  
    10231074
    10241075    for (int j = 0; j < nSamples; j++) {
    1025         double y = bounds->y0 + (j * deltaY / nSamples);
    1026         for (int i =  0; i < nSamples; i++) {
    1027 
    1028             psSphere srcSky;
    1029             psPlane *srcChip = psPlaneAlloc();
    1030             psPlane *dstTP = psPlaneAlloc();
    1031 
    1032             srcChip->x = bounds->x0 + (i * deltaX / nSamples);
    1033             srcChip->y = y;
    1034 
    1035             psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
    1036             psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
    1037             psDeproject (&srcSky, &srcTP, inFPA->toSky);
    1038 
    1039             // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
    1040 
    1041             psProject (dstTP, &srcSky, outFPA->toSky);
    1042 
    1043             srcChip->x -= bounds->x0;
    1044             srcChip->y -= bounds->y0;
    1045             psArrayAdd (src, 100, srcChip);
    1046             psArrayAdd (dst, 100, dstTP);
    1047 
    1048             psFree(srcChip);  // drop our refs to s and d
    1049             psFree(dstTP);
    1050         }
     1076        double y = bounds->y0 + (j * deltaY / nSamples);
     1077        for (int i =  0; i < nSamples; i++) {
     1078
     1079            psSphere srcSky;
     1080            psPlane *srcChip = psPlaneAlloc();
     1081            psPlane *dstTP = psPlaneAlloc();
     1082
     1083            srcChip->x = bounds->x0 + (i * deltaX / nSamples);
     1084            srcChip->y = y;
     1085
     1086            psPlaneTransformApply (&srcFP, inChip->toFPA, srcChip);
     1087            psPlaneTransformApply (&srcTP, inFPA->toTPA, &srcFP);
     1088            psDeproject (&srcSky, &srcTP, inFPA->toSky);
     1089
     1090            // fprintf (stderr, "%f %f | %f %f | %f %f | %f %f\n", srcChip->x, srcChip->y, srcFP.x, srcFP.y, srcTP.x, srcTP.y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD);
     1091
     1092            psProject (dstTP, &srcSky, outFPA->toSky);
     1093
     1094            srcChip->x -= bounds->x0;
     1095            srcChip->y -= bounds->y0;
     1096            psArrayAdd (src, 100, srcChip);
     1097            psArrayAdd (dst, 100, dstTP);
     1098
     1099            psFree(srcChip);  // drop our refs to s and d
     1100            psFree(dstTP);
     1101        }
    10511102    }
    10521103
    … …  
    10561107
    10571108    if (!psPlaneTransformFit(newToFPA, src, dst, 0, 0)) {
    1058         psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
    1059         psFree(src);
    1060         psFree(dst);
    1061         return NULL;
     1109        psError(PS_ERR_UNKNOWN, false, "linear fit to transform failed");
     1110        psFree(src);
     1111        psFree(dst);
     1112        return NULL;
    10621113    }
    10631114
    … …  
    10651116    for (int i = 0; i < src->n; i++) {
    10661117
    1067         psSphere srcSky, dstSky;
    1068         psPlane *srcChip = src->data[i];
    1069         psPlane *dstTP   = dst->data[i];
    1070 
    1071         psPlaneTransformApply (&srcFP, newToFPA, srcChip);
    1072         psDeproject (&srcSky, &srcFP, outFPA->toSky);
    1073         psDeproject (&dstSky, dstTP, outFPA->toSky);
    1074 
    1075         double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
    1076         double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
    1077         fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
     1118        psSphere srcSky, dstSky;
     1119        psPlane *srcChip = src->data[i];
     1120        psPlane *dstTP   = dst->data[i];
     1121
     1122        psPlaneTransformApply (&srcFP, newToFPA, srcChip);
     1123        psDeproject (&srcSky, &srcFP, outFPA->toSky);
     1124        psDeproject (&dstSky, dstTP, outFPA->toSky);
     1125
     1126        double dX = (srcSky.r*PS_DEG_RAD - dstSky.r*PS_DEG_RAD)*3600.0;
     1127        double dY = (srcSky.d*PS_DEG_RAD - dstSky.d*PS_DEG_RAD)*3600.0;
     1128        fprintf (stderr, "%f %f | %f %f | %f %f | %f %f | %f %f | %f %f\n", dX, dY, srcChip->x, srcChip->y, srcFP.x, srcFP.y, dstTP->x, dstTP->y, srcSky.r*PS_DEG_RAD, srcSky.d*PS_DEG_RAD, dstSky.r*PS_DEG_RAD, dstSky.d*PS_DEG_RAD);
    10781129
    10791130    }
    … …  
    10861137    psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1);
    10871138    if (!newFromFPA) {
    1088         psFree(newToFPA);
    1089         psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
    1090         return false;
     1139        psFree(newToFPA);
     1140        psError(PS_ERR_UNKNOWN, false, "inversion of fit of output chip toFPA failed");
     1141        return false;
    10911142    }
    10921143
    … …  
    11111162
    11121163    if (!wcs)
    1113         return;
     1164        return;
    11141165    psFree (wcs->trans);
    11151166    psFree (wcs->toSky);
    … …  
    11331184/*****
    11341185
    1135 For mosaic astrometry, we need to have a starting set of projection terms in which the
    1136 chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
    1137 microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
    1138 result in substantial errors between the observed and predicted star positions on the focal
    1139 plane: this is the measurement of the optical distortion in the camera.  At the same time,
    1140 we need to carry around the transformations which allow us to make an accurate calculation
    1141 of the position of the stars based on the input (per-chip) astrometry.  These
    1142 transformations will allow us to match the raw and ref stars robustly.  To convert the
    1143 per-chip astrometry (which may have been calculated with a different plate scale for each
    1144 chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
    1145 the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
    1146 each chip (eg, cdelt1).  Thus, we need to carry around both the
     1186      For mosaic astrometry, we need to have a starting set of projection terms in which the
     1187      chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
     1188      microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
     1189      result in substantial errors between the observed and predicted star positions on the focal
     1190      plane: this is the measurement of the optical distortion in the camera.  At the same time,
     1191      we need to carry around the transformations which allow us to make an accurate calculation
     1192      of the position of the stars based on the input (per-chip) astrometry.  These
     1193      transformations will allow us to match the raw and ref stars robustly.  To convert the
     1194      per-chip astrometry (which may have been calculated with a different plate scale for each
     1195      chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
     1196      the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
     1197      each chip (eg, cdelt1).  Thus, we need to carry around both the
    11471198
    11481199*****/

  • trunk/psModules/src/detrend/pmFringeStats.c

    r34085 r39926  
    336336        dfPt[i] = 1.0 / medianSd->sampleStdev;
    337337
    338         psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f : %s\n", (int)region.x0, (int)region.x1,
    339                 (int)region.y0, (int)region.y1, fPt[i], dfPt[i], readout->parent->hdu->extname);
     338        if (readout->parent->hdu) {
     339          psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f : %s\n", (int)region.x0, (int)region.x1,
     340                  (int)region.y0, (int)region.y1, fPt[i], dfPt[i], readout->parent->hdu->extname);
     341        }
     342        else {
     343          psTrace("psModules.detrend", 7, "[%d:%d,%d:%d]: %f %f : THIS_IS_A_SPOOKY_GHOST_CELL\n", (int)region.x0, (int)region.x1,
     344                  (int)region.y0, (int)region.y1, fPt[i], dfPt[i]);
     345        }
    340346    }
    341347    psFree(sky);

  • trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

    r36859 r39926  
    480480# if (TESTCOPY)
    481481    psImageCopy (pcm->modelConvFlux, pcm->modelFlux, pcm->modelFlux->type.type);
    482 # else
     482# else // TESTCOPY
    483483    if (pcm->use1Dgauss) {
    484484
    … …  
    496496        psImageConvolveKernel (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psfFFT);
    497497    }
    498 # endif
     498# endif // TESTCOPY
    499499
    500500    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
    … …  
    505505# if (TESTCOPY)
    506506        psImageCopy (dmodelConv, dmodel, dmodel->type.type);
    507 # else
     507# else // TESTCOPY
    508508        if (pcm->use1Dgauss) {
    509509            if (USE_1D_CACHE) {
    … …  
    520520            psImageConvolveKernel (dmodelConv, dmodel, NULL, 0, pcm->psfFFT);
    521521        }
    522 # endif
    523     }
    524 # else
     522# endif // TESTCOPY
     523    }
     524# else // PRE_CONVOLVE
    525525    // convolve model image and derivative images with psf via FFT
    526526    psImageConvolveFFT (pcm->modelConvFlux, pcm->modelFlux, NULL, 0, pcm->psf);
    … …  
    547547    }
    548548# endif // PRE-CONVOLVE
    549 # else
     549# else // USE_FFT
    550550    // convolve model image and derivative images with psf direct
    551551    psImageConvolveDirect (pcm->modelConvFlux, pcm->modelFlux, pcm->psf);

  • trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c

    r36856 r39926  
    119119                        psMetadataAdd (row, PS_LIST_TAIL, "RA_REF",     PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r);
    120120                        psMetadataAdd (row, PS_LIST_TAIL, "DEC_REF",    PS_DATA_F64, "declination (deg, J2000)",     PM_DEG_RAD*ref->sky->d);
     121                        psMetadataAdd (row, PS_LIST_TAIL, "RA_RAW",     PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*raw->sky->r);
     122                        psMetadataAdd (row, PS_LIST_TAIL, "DEC_RAW",    PS_DATA_F64, "declination (deg, J2000)",     PM_DEG_RAD*raw->sky->d);
    121123                        psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_REF", PS_DATA_F32, "x fitted coord on chip",       ref->chip->x);
    122124                        psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_REF", PS_DATA_F32, "y fitted coord on chip",       ref->chip->y);
Note: See TracChangeset for help on using the changeset viewer.