IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 22, 2007, 7:53:15 AM (18 years ago)
Author:
eugene
Message:

extensive work to finish

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15884 r15897  
    1010 *
    1111 *  @author EAM, IfA
    12  *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2007-12-19 18:57:05 $
     12 *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2007-12-22 17:53:15 $
    1414 *
    1515 *  Copyright 2007 Institute for Astronomy, University of Hawaii
    … …  
    120120    }
    121121
     122    if (!pmAstromWriteChips (file)) {
     123        psError(PS_ERR_IO, false, "Failed to write Astrometry for chips");
     124        return false;
     125    }
     126
     127    if (!pmAstromWriteFP (file)) {
     128        psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table");
     129        return false;
     130    }
     131
     132    if (!pmAstromWriteTP (file)) {
     133        psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table");
     134        return false;
     135    }
     136
    122137    if (!pmAstromWriteSky (file)) {
    123138        psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table");
    … …  
    125140    }
    126141
    127     if (!pmAstromWriteTP (file)) {
    128         psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table");
    129         return false;
    130     }
    131 
    132     if (!pmAstromWriteFP (file)) {
    133         psError(PS_ERR_IO, false, "Failed to write Sky for Astrometry table");
    134         return false;
    135     }
    136 
    137     if (!pmAstromWriteChips (file)) {
    138         psError(PS_ERR_IO, false, "Failed to write Astrometry for chips");
    139         return false;
    140     }
    141142    return true;
    142143}
    … …  
    174175    psFree (outhead);
    175176
    176     return true;
    177 }
    178 
    179 // first layer is the sky
    180 bool pmAstromWriteSky (pmFPAfile *file) {
    181 
    182     psMetadata *header = psMetadataAlloc();
    183     psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "SKY");
    184     psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
    185     psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "NONE");
    186     psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "NONE");
    187 
    188     psArray *table = psArrayAllocEmpty (1);
    189     psMetadata *row = psMetadataAlloc ();
    190     psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "SKY");
    191     psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
    192    
    193     psArrayAdd (table, 100, row);
    194     psFree (row);
    195 
    196     if (!psFitsWriteTable (file->fits, header, table, "SKY")) {
    197         psError(PS_ERR_IO, false, "writing sky data\n");
    198         psFree(table);
    199         return false;
    200     }
    201 
    202     psFree (table);
    203     psFree (header);
    204     return (true);
    205 }
    206 
    207 // second layer is the tangent plane
    208 bool pmAstromWriteTP (pmFPAfile *file) {
    209 
    210     psMetadata *header = psMetadataAlloc();
    211     psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "TANGENT_PLANE");
    212     psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
    213     psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
    214     psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "PROJECTION");
    215 
    216     psArray *table = psArrayAllocEmpty (1);
    217     psMetadata *row = psMetadataAlloc ();
    218     psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "TANGENT_PLANE");
    219     psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
    220 
    221     psRegion *region = pmAstromFPAExtent (file->fpa);
    222     psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
    223     psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
    224     psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
    225     psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
    226    
    227     psMetadataAddF32(row,    PS_LIST_TAIL, "XSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Xs * PM_DEG_RAD);
    228     psMetadataAddF32(row,    PS_LIST_TAIL, "YSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Ys * PM_DEG_RAD);
    229     psMetadataAddF32(row,    PS_LIST_TAIL, "XREF",     PS_META_REPLACE, "", file->fpa->toSky->R  * PM_DEG_RAD);
    230     psMetadataAddF32(row,    PS_LIST_TAIL, "YREF",     PS_META_REPLACE, "", file->fpa->toSky->D  * PM_DEG_RAD);
    231 
    232     psArrayAdd (table, 100, row);
    233     psFree (row);
    234 
    235     if (!psFitsWriteTable (file->fits, header, table, "TP")) {
    236         psError(PS_ERR_IO, false, "writing sky data\n");
    237         psFree (region);
    238         psFree (table);
    239         psFree (header);
    240         return false;
    241     }
    242 
    243     psFree (region);
    244     psFree (table);
    245     psFree (header);
    246     return (true);
    247 }
    248 
    249 // third layer is the focal plane
    250 bool pmAstromWriteFP (pmFPAfile *file) {
    251 
    252     psMetadata *header = psMetadataAlloc();
    253     psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "FOCAL_PLANE");
    254     psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
    255     psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
    256     psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "POLYNOMIAL");
    257 
    258     psArray *table = psArrayAllocEmpty (1);
    259 
    260     // XXX is this or the tpa region correct?
    261     psRegion *region = pmAstromFPAExtent (file->fpa);
    262 
    263     for (int i = 0; i <= file->fpa->toTPA->x->nX; i++) {
    264         for (int j = 0; j <= file->fpa->toTPA->x->nY; j++) {
    265             psMetadata *row = psMetadataAlloc ();
    266             psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "FOCAL_PLANE");
    267             psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
    268             psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
    269             psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
    270             psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
    271             psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
    272    
    273             psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", file->fpa->toTPA->x->nX);
    274             psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", file->fpa->toTPA->x->nY);
    275             psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", file->fpa->toTPA->x->coeff[i][j]);
    276             psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", file->fpa->toTPA->y->coeff[i][j]);
    277             psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_X",  PS_META_REPLACE, "", file->fpa->toTPA->x->coeffErr[i][j]);
    278             psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_Y",  PS_META_REPLACE, "", file->fpa->toTPA->y->coeffErr[i][j]);
    279 
    280             psArrayAdd (table, 100, row);
    281             psFree (row);
    282         }
    283     }
    284 
    285     if (!psFitsWriteTable (file->fits, header, table, "FP")) {
    286         psError(PS_ERR_IO, false, "writing sky data\n");
    287         psFree(table);
    288         return false;
    289     }
    290 
    291     psFree (table);
    292     psFree (header);
    293     psFree (region);
    294177    return true;
    295178}
    … …  
    332215                psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
    333216   
    334                 psMetadataAddS32(row,    PS_LIST_TAIL, "NX",       PS_META_REPLACE, "", i);
    335                 psMetadataAddS32(row,    PS_LIST_TAIL, "NY",       PS_META_REPLACE, "", j);
     217                psMetadataAddS32(row,    PS_LIST_TAIL, "XORDER",   PS_META_REPLACE, "", i);
     218                psMetadataAddS32(row,    PS_LIST_TAIL, "YORDER",   PS_META_REPLACE, "", j);
     219                psMetadataAddS32(row,    PS_LIST_TAIL, "NXORDER",  PS_META_REPLACE, "", chip->toFPA->x->nX);
     220                psMetadataAddS32(row,    PS_LIST_TAIL, "NYORDER",  PS_META_REPLACE, "", chip->toFPA->x->nY);
    336221                psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", chip->toFPA->x->coeff[i][j]);
    337222                psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", chip->toFPA->y->coeff[i][j]);
    338223                psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_X",  PS_META_REPLACE, "", chip->toFPA->x->coeffErr[i][j]);
    339224                psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_Y",  PS_META_REPLACE, "", chip->toFPA->y->coeffErr[i][j]);
     225                psMetadataAddU8 (row,    PS_LIST_TAIL, "MASK_X",   PS_META_REPLACE, "", chip->toFPA->x->coeffMask[i][j]);
     226                psMetadataAddU8 (row,    PS_LIST_TAIL, "MASK_Y",   PS_META_REPLACE, "", chip->toFPA->y->coeffMask[i][j]);
    340227                psArrayAdd (table, 100, row);
    341228                psFree (row);
    … …  
    357244    psFree (header);
    358245    return true;
     246}
     247
     248// third layer is the focal plane
     249bool pmAstromWriteFP (pmFPAfile *file) {
     250
     251    bool status;
     252
     253    psMetadata *header = psMetadataAlloc();
     254    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "FOCAL_PLANE");
     255    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     256    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
     257    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "POLYNOMIAL");
     258
     259    psArray *table = psArrayAllocEmpty (1);
     260
     261    // XXX is this or the tpa region correct?
     262    psRegion *region = pmAstromFPAExtent (file->fpa);
     263
     264    // rotate the toTPA to have 0.0 posangle
     265    float posangle = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.POSANGLE");
     266
     267    // the to/from TPA transform currently has rotation of posangle; remove it & create the toTPA version
     268    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, -posangle);
     269    psPlaneTransform *toTPA   = psPlaneTransformInvert(NULL, fromTPA, *region, 50);
     270
     271    for (int i = 0; i <= toTPA->x->nX; i++) {
     272        for (int j = 0; j <= toTPA->x->nY; j++) {
     273            psMetadata *row = psMetadataAlloc ();
     274            psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "FOCAL_PLANE");
     275            psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "TANGENT_PLANE");
     276            psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
     277            psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
     278            psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
     279            psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     280   
     281            psMetadataAddS32(row,    PS_LIST_TAIL, "XORDER",   PS_META_REPLACE, "", i);
     282            psMetadataAddS32(row,    PS_LIST_TAIL, "YORDER",   PS_META_REPLACE, "", j);
     283            psMetadataAddS32(row,    PS_LIST_TAIL, "NXORDER",  PS_META_REPLACE, "", toTPA->x->nX);
     284            psMetadataAddS32(row,    PS_LIST_TAIL, "NYORDER",  PS_META_REPLACE, "", toTPA->x->nY);
     285            psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_X",   PS_META_REPLACE, "", toTPA->x->coeff[i][j]);
     286            psMetadataAddF32(row,    PS_LIST_TAIL, "POLY_Y",   PS_META_REPLACE, "", toTPA->y->coeff[i][j]);
     287            psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_X",  PS_META_REPLACE, "", toTPA->x->coeffErr[i][j]);
     288            psMetadataAddF32(row,    PS_LIST_TAIL, "ERROR_Y",  PS_META_REPLACE, "", toTPA->y->coeffErr[i][j]);
     289            psMetadataAddU8 (row,    PS_LIST_TAIL, "MASK_X",   PS_META_REPLACE, "", toTPA->x->coeffMask[i][j]);
     290            psMetadataAddU8 (row,    PS_LIST_TAIL, "MASK_Y",   PS_META_REPLACE, "", toTPA->y->coeffMask[i][j]);
     291
     292            psArrayAdd (table, 100, row);
     293            psFree (row);
     294        }
     295    }
     296
     297    if (!psFitsWriteTable (file->fits, header, table, "FP")) {
     298        psError(PS_ERR_IO, false, "writing sky data\n");
     299        psFree(table);
     300        psFree (header);
     301        psFree (region);
     302        return false;
     303    }
     304
     305    psFree (table);
     306    psFree (header);
     307    psFree (region);
     308    return true;
     309}
     310
     311// second layer is the tangent plane
     312bool pmAstromWriteTP (pmFPAfile *file) {
     313
     314    bool status;
     315
     316    // get the boresite model parameters
     317    float Xo = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.X0");
     318    float Yo = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.Y0");
     319    float RX = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.RX");
     320    float RY = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.RY");
     321    float To = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.T0");
     322    float Po = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.BORE.P0");
     323    float PosZero = psMetadataLookupF32 (&status, file->fpa->concepts, "FPA.POS_ZERO");  /// XXX be consistent with degrees v radians
     324    char *refChip = psMetadataLookupStr (&status, file->fpa->concepts, "FPA.REF.CHIP");
     325
     326    psMetadata *header = psMetadataAlloc();
     327    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "TANGENT_PLANE");
     328    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
     329    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "RECTANGLE");
     330    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "PROJECTION");
     331
     332    psArray *table = psArrayAllocEmpty (1);
     333    psMetadata *row = psMetadataAlloc ();
     334    psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "TANGENT_PLANE");
     335    psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "SKY");
     336
     337    psRegion *region = pmAstromFPAExtent (file->fpa);
     338    psMetadataAddF32(row,    PS_LIST_TAIL, "MINX",     PS_META_REPLACE, "range", region->x0);
     339    psMetadataAddF32(row,    PS_LIST_TAIL, "MAXX",     PS_META_REPLACE, "range", region->x1);
     340    psMetadataAddF32(row,    PS_LIST_TAIL, "MINY",     PS_META_REPLACE, "range", region->y0);
     341    psMetadataAddF32(row,    PS_LIST_TAIL, "MAXY",     PS_META_REPLACE, "range", region->y1);
     342   
     343    // psMetadataAddF32(row,    PS_LIST_TAIL, "XREF",     PS_META_REPLACE, "", file->fpa->toSky->R  * PM_DEG_RAD);
     344    // psMetadataAddF32(row,    PS_LIST_TAIL, "YREF",     PS_META_REPLACE, "", file->fpa->toSky->D  * PM_DEG_RAD);
     345
     346    psMetadataAddF32(row,    PS_LIST_TAIL, "XSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Xs * PM_DEG_RAD);
     347    psMetadataAddF32(row,    PS_LIST_TAIL, "YSCALE",   PS_META_REPLACE, "", file->fpa->toSky->Ys * PM_DEG_RAD);
     348    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_X0",  PS_META_REPLACE, "boresite parameter", Xo);
     349    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_Y0",  PS_META_REPLACE, "boresite parameter", Yo);
     350    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_RX",  PS_META_REPLACE, "boresite parameter", RX);
     351    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_RY",  PS_META_REPLACE, "boresite parameter", RY);
     352    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_T0",  PS_META_REPLACE, "boresite parameter", To);
     353    psMetadataAddF32(row,    PS_LIST_TAIL, "BORE_P0",  PS_META_REPLACE, "boresite parameter", Po);
     354    psMetadataAddF32(row,    PS_LIST_TAIL, "POS_ZERO", PS_META_REPLACE, "boresite parameter", PosZero);
     355    psMetadataAddStr(row,    PS_LIST_TAIL, "REF_CHIP", PS_META_REPLACE, "boresite parameter", refChip);
     356
     357    psArrayAdd (table, 100, row);
     358    psFree (row);
     359
     360    if (!psFitsWriteTable (file->fits, header, table, "TP")) {
     361        psError(PS_ERR_IO, false, "writing sky data\n");
     362        psFree (region);
     363        psFree (table);
     364        psFree (header);
     365        return false;
     366    }
     367
     368    psFree (region);
     369    psFree (table);
     370    psFree (header);
     371    return (true);
     372}
     373
     374// first layer is the sky
     375bool pmAstromWriteSky (pmFPAfile *file) {
     376
     377    psMetadata *header = psMetadataAlloc();
     378    psMetadataAddStr(header, PS_LIST_TAIL, "COORD",    PS_META_REPLACE, "name of this layer",   "SKY");
     379    psMetadataAddStr(header, PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
     380    psMetadataAddStr(header, PS_LIST_TAIL, "BOUNDARY", PS_META_REPLACE, "validity region",      "NONE");
     381    psMetadataAddStr(header, PS_LIST_TAIL, "TRANSFRM", PS_META_REPLACE, "mapping to parent",    "NONE");
     382
     383    psArray *table = psArrayAllocEmpty (1);
     384    psMetadata *row = psMetadataAlloc ();
     385    psMetadataAddStr(row,    PS_LIST_TAIL, "SEGMENT",  PS_META_REPLACE, "name of this segment", "SKY");
     386    psMetadataAddStr(row,    PS_LIST_TAIL, "PARENT",   PS_META_REPLACE, "next layer up",        "NONE");
     387   
     388    psArrayAdd (table, 100, row);
     389    psFree (row);
     390
     391    if (!psFitsWriteTable (file->fits, header, table, "SKY")) {
     392        psError(PS_ERR_IO, false, "writing sky data\n");
     393        psFree(table);
     394        return false;
     395    }
     396
     397    psFree (table);
     398    psFree (header);
     399    return (true);
    359400}
    360401
    … …  
    504545        chip->toFPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X");
    505546        chip->toFPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y");
     547        chip->toFPA->x->coeffMask[ix][iy] = psMetadataLookupU8(&status, row, "MASK_X");
     548        chip->toFPA->y->coeffMask[ix][iy] = psMetadataLookupU8(&status, row, "MASK_Y");
    506549    }
    507550
    … …  
    509552    for (int i = 0; i < file->fpa->chips->n; i++) {
    510553        pmChip *chip = file->fpa->chips->data[i];
     554        if (!chip->toFPA) continue;
    511555        psRegion *region = pmChipPixels (chip);
    512556        psFree (chip->fromFPA);
    … …  
    533577    if (!header) psAbort("cannot read table header");
    534578
    535     // there is only one transformation in this table; the order is defined in the header
    536     int nX = psMetadataLookupS32(&status, header, "NXORDER"); REQUIRE (status, "missing NXORDER");
    537     int nY = psMetadataLookupS32(&status, header, "NYORDER"); REQUIRE (status, "missing NYORDER");
    538 
    539     // allocate the new transformation
     579    // free the old
    540580    psFree (file->fpa->toTPA);
    541     file->fpa->toTPA = psPlaneTransformAlloc(nX, nY);
    542 
    543     // free & DON'T allocate the new transformation fromTPA; this is created after the boresite
    544     // is applied in pmAstromReadTP
    545     psFree (file->fpa->fromTPA);
    546     file->fpa->fromTPA = NULL;
     581    file->fpa->toTPA = NULL;
    547582
    548583    // read the complete table data at one shot
    … …  
    553588    for (int i = 0; i < table->n; i++) {
    554589        psMetadata *row = table->data[i];
     590
     591        // there is only one transformation in this table; the order is defined in the header
     592        int nX = psMetadataLookupS32(&status, row, "NXORDER"); REQUIRE (status, "missing NXORDER");
     593        int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER");
     594        if (file->fpa->toTPA == NULL) {
     595            // allocate the new transformation
     596            file->fpa->toTPA = psPlaneTransformAlloc(nX, nY);
     597        } else {
     598            REQUIRE (file->fpa->toTPA->x->nX == nX, "mismatch in chip order");
     599            REQUIRE (file->fpa->toTPA->x->nY == nY, "mismatch in chip order");
     600            REQUIRE (file->fpa->toTPA->y->nX == nX, "mismatch in chip order");
     601            REQUIRE (file->fpa->toTPA->y->nY == nY, "mismatch in chip order");
     602        }
     603
    555604        int ix = psMetadataLookupS32(&status, row, "XORDER"); REQUIRE (status, "missing XORDER");
    556605        int iy = psMetadataLookupS32(&status, row, "YORDER"); REQUIRE (status, "missing YORDER");
    557         file->fpa->toTPA->x->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_X"); REQUIRE (status, "missing POLY_X");
    558         file->fpa->toTPA->y->coeff[ix][iy]    = psMetadataLookupF32(&status, row, "POLY_Y"); REQUIRE (status, "missing POLY_Y");
    559         file->fpa->toTPA->x->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_X"); REQUIRE (status, "missing ERROR_X");
    560         file->fpa->toTPA->y->coeffErr[ix][iy] = psMetadataLookupF32(&status, row, "ERROR_Y"); REQUIRE (status, "missing ERROR_Y");
    561     }
     606        file->fpa->toTPA->x->coeff[ix][iy]     = psMetadataLookupF32(&status, row, "POLY_X");  REQUIRE (status, "missing POLY_X");
     607        file->fpa->toTPA->y->coeff[ix][iy]     = psMetadataLookupF32(&status, row, "POLY_Y");  REQUIRE (status, "missing POLY_Y");
     608        file->fpa->toTPA->x->coeffErr[ix][iy]  = psMetadataLookupF32(&status, row, "ERROR_X"); REQUIRE (status, "missing ERROR_X");
     609        file->fpa->toTPA->y->coeffErr[ix][iy]  = psMetadataLookupF32(&status, row, "ERROR_Y"); REQUIRE (status, "missing ERROR_Y");
     610        file->fpa->toTPA->x->coeffMask[ix][iy] = psMetadataLookupU8 (&status, row, "MASK_X");  REQUIRE (status, "missing MASK_X");
     611        file->fpa->toTPA->y->coeffMask[ix][iy] = psMetadataLookupU8 (&status, row, "MASK_Y");  REQUIRE (status, "missing MASK_Y");
     612    }
     613
     614    psRegion *region = pmAstromFPAExtent (file->fpa);
     615    psFree (file->fpa->fromTPA);
     616    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50);
    562617
    563618    psFree (table);
    564619    psFree (header);
    565     return true;
    566 }
     620    psFree (region);
     621    return true;
     622}
     623
     624# define TRANSFER(TO,FROM,NAME) { \
     625    psMetadataItem *item = psMetadataLookup(FROM,NAME); \
     626    if (!item) psAbort ("cannot find %s", NAME); \
     627    psMetadataItem *newItem = psMetadataItemCopy(item); \
     628    if (!psMetadataAddItem(TO, newItem, PS_LIST_TAIL, PS_META_REPLACE)) { \
     629        psAbort ("cannot copy %s", NAME); \
     630    } \
     631    psFree (newItem); }
    567632
    568633// third layer applies boresite corrections and converts tangent plane to sky
    569634bool pmAstromReadTP (pmFPAfile *file) {
    570635   
    571     bool status;
    572 
    573     // these externally supplied values are used to set the final transformation terms
    574     double RA  = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA");
    575     double DEC = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.DEC"); REQUIRE (status, "missing FPA.DEC");
    576     double POS = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.POSANGLE"); REQUIRE (status, "missing FPA.POSANGLE");
    577     double ROT = psMetadataLookupF64 (&status, file->fpa->concepts, "FPA.ROTANGLE"); REQUIRE (status, "missing FPA.ROTANGLE");
    578 
    579636    if (!psFitsMoveExtName (file->fits, "TP")) {
    580637        psError(PS_ERR_IO, false, "missing TP extension in astrometry table\n");
    … …  
    591648    psMetadata *row = table->data[0];
    592649
    593     // get projection scale; center is supplied
    594     float Xs = psMetadataLookupF32(&status, row, "XSCALE") * PM_RAD_DEG; REQUIRE (status, "missing XSCALE");
    595     float Ys = psMetadataLookupF32(&status, row, "YSCALE") * PM_RAD_DEG; REQUIRE (status, "missing YSCALE");
    596 
    597     // allocate a new toSky projection
    598     psFree (file->fpa->toSky);
    599     file->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_WRP);
    600 
    601     // get boresite correction terms.  L,M,T define an ellipse along which the boresite travels
    602     double R_L = psMetadataLookupF32(&status, row, "ROTATOR_L"); REQUIRE (status, "missing ROTATOR_L");
    603     double R_M = psMetadataLookupF32(&status, row, "ROTATOR_M"); REQUIRE (status, "missing ROTATOR_M");
    604     double R_T = psMetadataLookupF32(&status, row, "ROTATOR_T"); REQUIRE (status, "missing ROTATOR_T");
    605 
    606     // the position of the boresite on the ellipse is the ROTANGLE - ROT_ZERO
    607     double R_0 = psMetadataLookupF32(&status, row, "ROT_ZERO"); REQUIRE (status, "missing ROT_ZERO");
    608 
    609     // find the current position of the boresite for this image
    610     double Lo = R_L*cos(ROT - R_0)*cos(R_T) + R_M*sin(ROT - R_0)*sin(R_T);
    611     double Mo = R_M*sin(ROT - R_0)*cos(R_T) - R_L*cos(ROT - R_0)*sin(R_T);
    612 
    613     // apply new reference pixel to transformation
    614     psPlaneTransform *tmpToTPA = psPlaneTransformSetCenter (NULL, file->fpa->toTPA, Lo, Mo);
    615     psFree (file->fpa->toTPA);
    616    
    617     // apply POSANGLE
    618     double Po = psMetadataLookupF32(&status, row, "POS_ZERO"); REQUIRE (status, "missing POS_ZERO");
    619     file->fpa->toTPA = psPlaneTransformRotate (NULL, tmpToTPA, POS - Po);
    620     psFree (tmpToTPA);
    621 
    622     psRegion *region = pmAstromFPAExtent (file->fpa);
    623     psFree (file->fpa->fromTPA);
    624     file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50);
    625 
    626     psFree (region);
     650    // move needed items to the concepts
     651    TRANSFER (file->fpa->concepts, row, "XSCALE");
     652    TRANSFER (file->fpa->concepts, row, "XSCALE");
     653    TRANSFER (file->fpa->concepts, row, "YSCALE");
     654    TRANSFER (file->fpa->concepts, row, "BORE_X0");
     655    TRANSFER (file->fpa->concepts, row, "BORE_Y0");
     656    TRANSFER (file->fpa->concepts, row, "BORE_RX");
     657    TRANSFER (file->fpa->concepts, row, "BORE_RY");
     658    TRANSFER (file->fpa->concepts, row, "BORE_T0");
     659    TRANSFER (file->fpa->concepts, row, "BORE_P0");
     660    TRANSFER (file->fpa->concepts, row, "POS_ZERO");
     661    TRANSFER (file->fpa->concepts, row, "REF_CHIP");
     662
    627663    psFree (table);
    628664    psFree (header);
    … …  
    652688}
    653689
     690// third layer applies boresite corrections and converts tangent plane to sky
     691bool pmAstromSetTP (pmFPAfile *file, psMetadata *concepts) {
     692   
     693    bool status;
     694
     695    // these externally supplied values are used to set the final transformation terms
     696    double RA  = psMetadataLookupF64 (&status, concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA");
     697    double DEC = psMetadataLookupF64 (&status, concepts, "FPA.DEC"); REQUIRE (status, "missing FPA.DEC");
     698    double POS = PM_RAD_DEG * psMetadataLookupF64 (&status, concepts, "FPA.POSANGLE"); REQUIRE (status, "missing FPA.POSANGLE");
     699    // double ROT = psMetadataLookupF64 (&status, concepts, "FPA.ROTANGLE"); REQUIRE (status, "missing FPA.ROTANGLE");
     700
     701    // get projection scale; center is supplied
     702    float Xs = psMetadataLookupF32(&status, file->fpa->concepts, "XSCALE") * PM_RAD_DEG; REQUIRE (status, "missing XSCALE");
     703    float Ys = psMetadataLookupF32(&status, file->fpa->concepts, "YSCALE") * PM_RAD_DEG; REQUIRE (status, "missing YSCALE");
     704
     705    // allocate a new toSky projection using the reported position
     706    psFree (file->fpa->toSky);
     707    file->fpa->toSky = psProjectionAlloc (RA, DEC, Xs, Ys, PS_PROJ_WRP);
     708
     709    // get boresite correction terms.  RX,RY,To,Po define an ellipse along which the boresite travels
     710    double Xo = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_X0"); REQUIRE (status, "missing ");
     711    double Yo = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_Y0"); REQUIRE (status, "missing ");
     712    double RX = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_RX"); REQUIRE (status, "missing ");
     713    double RY = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_RY"); REQUIRE (status, "missing ");
     714    double To = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_T0"); REQUIRE (status, "missing ");
     715    double Po = psMetadataLookupF32(&status, file->fpa->concepts, "BORE_P0"); REQUIRE (status, "missing ");
     716
     717    // the true rotation of the instrument is POSANGLE - POS_ZERO
     718    double PosZero = PM_RAD_DEG * psMetadataLookupF32(&status, file->fpa->concepts, "POS_ZERO"); REQUIRE (status, "missing POS_ZERO");
     719    char *refChip  = psMetadataLookupStr(&status, file->fpa->concepts, "REF_CHIP"); REQUIRE (status, "missing REF_CHIP");
     720
     721    // apply true posangle = -(POS - POS_ZERO)
     722    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, (POS - PosZero));
     723    psFree (file->fpa->fromTPA);
     724    file->fpa->fromTPA = fromTPA;
     725
     726    psRegion *region = pmAstromFPAExtent (file->fpa);
     727    psFree (file->fpa->toTPA);
     728    file->fpa->toTPA = psPlaneTransformInvert(NULL, file->fpa->fromTPA, *region, 50);
     729
     730    // current position of the nominal boresite in refChip coordinates
     731    double X = Xo + RX*cos(POS - To)*cos(Po) + RY*sin(POS - To)*sin(Po);
     732    double Y = Yo + RY*sin(POS - To)*cos(Po) - RX*cos(POS - To)*sin(Po);
     733
     734    // get reference chip from name
     735    pmChip *chip = pmConceptsChipFromName (file->fpa, refChip);
     736    if (!chip) psAbort ("invalid chip name for reference");
     737
     738    psPlane *boreCH = psPlaneAlloc();
     739    psPlane *boreFP = psPlaneAlloc();
     740    psPlane *boreTP = psPlaneAlloc();
     741    psSphere *boreSky = psSphereAlloc();
     742
     743    // find the true RA,DEC coord of the mirror of the reported boresite
     744    boreCH->x = -X;
     745    boreCH->y = -Y;
     746    psPlaneTransformApply (boreFP, chip->toFPA, boreCH);
     747    psPlaneTransformApply (boreTP, file->fpa->toTPA, boreFP);
     748    psDeproject (boreSky, boreTP, file->fpa->toSky);
     749
     750    // modify the projection to account for offset between true and reported boresite
     751    file->fpa->toSky->R = boreSky->r;
     752    file->fpa->toSky->D = boreSky->d;
     753
     754    psFree (boreCH);
     755    psFree (boreFP);
     756    psFree (boreTP);
     757    psFree (boreSky);
     758
     759    psFree (region);
     760    return (true);
     761}
     762
Note: See TracChangeset for help on using the changeset viewer.