IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7332


Ignore:
Timestamp:
Jun 4, 2006, 7:02:16 AM (20 years ago)
Author:
eugene
Message:

substantial updates

Location:
trunk/psastro
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/Makefile

    r7154 r7332  
    3939$(SRC)/psastroMosaicSetAstrom.$(ARCH).o     \
    4040$(SRC)/psastroMosaicSetMatch.$(ARCH).o      \
     41$(SRC)/psastroMosaicHeaders.$(ARCH).o       \
     42$(SRC)/psastroMosaicRescaleChips.$(ARCH).o          \
    4143$(SRC)/psastroWCS.$(ARCH).o       
    4244

  • trunk/psastro/src/psastro.c

    r7154 r7332  
    2626    psArray *refs = psastroLoadReferences (config);
    2727
    28     psastroChipAstrom (config, refs);
    29     psastroMosaicAstrom (config, refs);
     28    psastroMosaicGetRefstars (config, refs);
    3029
    31     # if 0
    32     // perform the desired astrometry analysis
    33     switch (mode) {
    34       case stack:
    35         psastroStackAstrom (config, refs);
    36         break
    37       case chip:
     30    char *mosastro = psMetadataLookupStr (NULL, config->arguments, "MOSASTRO");
     31    if (mosastro == NULL) {
    3832        psastroChipAstrom (config, refs);
    39         break
    40       case mosaic:
     33    } else {
    4134        psastroMosaicAstrom (config, refs);
    42         break;
    43       default:
    44         break;
    4535    }
    46     # endif
     36
     37    // XXX how do we specify stack astrometry?
     38    // psastroStackAstrom (config, refs);
    4739
    4840    // write out the results
    … …  
    6052    pmConceptsDone ();
    6153    pmConfigDone ();
    62     fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    63     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
     54    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
     55    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");
    6456
    6557    exit (0);

  • trunk/psastro/src/psastro.h

    r7154 r7332  
    2323bool              psastroAstromGuess (pmConfig *config);
    2424
    25 bool              pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale);
     25// bool              pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale);
     26bool              pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic);
    2627bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale);
    2728psPlaneDistort   *psPlaneDistortIdentity ();
    … …  
    4142psArray *psastroMosaicGetGrads (pmFPA *fpa, psMetadata *recipe);
    4243bool psastroMosaicAstrom (pmConfig *config, psArray *refs);
    43 bool psastroMosaicGetRefstars (pmFPA *fpa, psArray *refs, psMetadata *recipe);
     44bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs);
    4445bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe);
    4546bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe);
    4647bool psastroMosaicSetAstrom (pmFPA *fpa);
     48bool psastroMosaicHeaders (pmConfig *config);
     49bool psastroMosaicRescaleChips (pmFPA *fpa);
     50
     51bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale);
     52psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale);

  • trunk/psastro/src/psastroArguments.c

    r7014 r7332  
    4545    }
    4646
     47    // mosastro mode also specifies output header file
     48    if ((N = psArgumentGet (*argc, argv, "-mosastro"))) {
     49        psArgumentRemove (N, argc, argv);
     50        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MOSASTRO", PS_DATA_STRING, "", psStringCopy(argv[N]));
     51        psArgumentRemove (N, argc, argv);
     52    }
     53
    4754    status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list");
    4855    if (!status) { usage ();}

  • trunk/psastro/src/psastroAstromGuess.c

    r7014 r7332  
    33// XXX this function assumes the initial astrometry arrives in the form of
    44// XXX WCS keywords in the headers corresponding to the chips.
     5// XXX this function is both converting the header WCS astrometry terms to the fpa terms
     6//     and applying the astrometry to the detected objects. 
    57bool psastroAstromGuess (pmConfig *config) {
    68
    79    bool newFPA = true;
    810    bool status = false;
     11    bool isMosaic = false;
    912    double RAmin, RAmax, RAminSky, RAmaxSky;
    1013    double DECmin, DECmax;
    … …  
    1316    pmCell *cell = NULL;
    1417    pmReadout *readout = NULL;
     18
     19    // is this a mosaic astrometry analysis?
     20    psMetadataLookupStr (&isMosaic, config->arguments, "MOSASTRO");
    1521
    1622    // select the current recipe
    … …  
    3036    double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE");
    3137    if (!status) plateScale = 1.0;
     38    plateScale = 1.0;
    3239
    3340    pmFPAview *view = pmFPAviewAlloc (0);
    … …  
    4148        pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    4249
    43         pmAstromReadWCS (&fpa->projection, &fpa->toTangentPlane, &chip->toFPA, hdu->header, plateScale);
     50        pmAstromReadWCS (fpa, chip, hdu->header, plateScale, isMosaic);
    4451        if (newFPA) {
    4552            newFPA = false;
    … …  
    4956            DECmin = DECmax = fpa->projection->D;
    5057        }
    51 
    52         // build projection inversions
    53         // XXX region should be set from image data
    54         psRegion region = {0, 0, 2000, 4500};
    55         fpa->fromTangentPlane = psPlaneDistortIdentity ();
    56         chip->fromFPA = psPlaneTransformInvert (NULL, chip->toFPA, region, 20);
    5758
    5859        // apply the new WCS guess data to all of the data in the readouts
    … …  
    6667                if (! readout->data_exists) { continue; }
    6768
    68                 psArray *objects = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
    69                 if (objects == NULL) { continue; }
     69                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     70                if (rawstars == NULL) { continue; }
    7071
    71                 for (int i = 0; i < objects->n; i++) {
    72                     pmAstromObj *object = objects->data[i];
     72                for (int i = 0; i < rawstars->n; i++) {
     73                    pmAstromObj *raw = rawstars->data[i];
    7374       
    74                     psPlaneTransformApply (object->FP, chip->toFPA, object->chip);
    75                     psPlaneDistortApply (object->TP, fpa->toTangentPlane, object->FP, 0.0, 0.0);
    76                     p_psDeproject (object->sky, object->TP, fpa->projection);
     75                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
     76                    psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
     77                    p_psDeproject (raw->sky, raw->TP, fpa->projection);
     78
     79                    if (i < 0) {
     80                        fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n",
     81                                 raw->chip->x, raw->chip->y,
     82                                 raw->FP->x, raw->FP->y,
     83                                 raw->TP->x, raw->TP->y,
     84                                 raw->sky->r, raw->sky->d);
     85
     86                        psPlane *fp = psPlaneAlloc();
     87                        psPlane *tp = psPlaneAlloc();
     88                        psPlane *ch = psPlaneAlloc();
     89
     90                        p_psProject (tp, raw->sky, fpa->projection);
     91                        psPlaneDistortApply (fp, fpa->fromTangentPlane, tp, 0.0, 0.0);
     92                        psPlaneTransformApply (ch, chip->fromFPA, fp);
     93
     94                        fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n",
     95                                 ch->x, ch->y,
     96                                 fp->x, fp->y,
     97                                 tp->x, tp->y,
     98                                 raw->sky->r, raw->sky->d);
     99                    }
    77100
    78101                    // rationalize ra to sky range centered on boresite
    79                     while (object->sky->r < RAminSky) object->sky->r += M_PI;
    80                     while (object->sky->r > RAmaxSky) object->sky->r -= M_PI;
     102                    while (raw->sky->r < RAminSky) raw->sky->r += M_PI;
     103                    while (raw->sky->r > RAmaxSky) raw->sky->r -= M_PI;
    81104
    82                     RAmin = PS_MIN (object->sky->r, RAmin);
    83                     RAmax = PS_MAX (object->sky->r, RAmax);
     105                    RAmin = PS_MIN (raw->sky->r, RAmin);
     106                    RAmax = PS_MAX (raw->sky->r, RAmax);
    84107
    85                     DECmin = PS_MIN (object->sky->d, DECmin);
    86                     DECmax = PS_MAX (object->sky->d, DECmax);
     108                    DECmin = PS_MIN (raw->sky->d, DECmin);
     109                    DECmax = PS_MAX (raw->sky->d, DECmax);
    87110                }
    88111            }

  • trunk/psastro/src/psastroChipAstrom.c

    r7014 r7332  
    4141
    4242                // select the raw objects for this readout
    43                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     43                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    4444                if (rawstars == NULL) { continue; }
    4545
    46                 // the refstars is a subset within range of this chip
    47                 psArray *refstars = psArrayAlloc (100);
     46                // select the raw objects for this readout
     47                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     48                if (refstars == NULL) { continue; }
    4849
    49                 // select the reference objects within range of this readout
    50                 // project the reference objects to this chip
    51                 for (int i = 0; i < refs->n; i++) {
    52                     pmAstromObj *ref = refs->data[i];
    53        
    54                     p_psProject (ref->TP, ref->sky, fpa->projection);
    55                     psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0);
    56                     psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    57 
    58                     // XXX what is the X,Y range of the selected chip?
    59                     // XXX for the moment, force these to be good for Megacam
    60                     if (ref->chip->x < -500) continue;
    61                     if (ref->chip->x > 2500) continue;
    62                     if (ref->chip->y < -500) continue;
    63                     if (ref->chip->y > 5000) continue;
    64                    
    65                     psArrayAdd (refstars, 100, ref);
    66                 }
    6750                psastroOneChip (fpa, chip, refstars, rawstars, recipe);
    68                 psFree (refstars);
    6951
    7052                // read WCS data from the corresponding header

  • trunk/psastro/src/psastroConvert.c

    r7014 r7332  
    99        return false;
    1010
    11     psArray *objects = pmSourceToAstromObj (sources);
     11    psArray *rawstars = pmSourceToAstromObj (sources);
    1212
    13     psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.OBJECTS", PS_DATA_ARRAY, "astrometry objects", objects);
    14     psFree (objects);
     13    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars);
     14    psFree (rawstars);
    1515 
    1616    return true;

  • trunk/psastro/src/psastroLoadReferences.c

    r7014 r7332  
    6262
    6363    // convert the Average table to the pmAstromObj entries
    64     psArray *objects = psArrayAlloc (table->n);
     64    psArray *refs = psArrayAlloc (table->n);
    6565    for (int i = 0; i < table->n; i++) {
    66         pmAstromObj *object = pmAstromObjAlloc ();
     66        pmAstromObj *ref = pmAstromObjAlloc ();
    6767
    6868        psMetadata *row = table->data[i];
    … …  
    7070        // DVO tables are stored in degrees
    7171        # if ELIXIR_MODE
    72         object->sky->r   = RAD_DEG*psMetadataLookupF32 (&status, row, "RA");
    73         object->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
    74         object->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // XXX ELIXIR uses millimags, PANSTARRS uses mags
     72        ref->sky->r   = RAD_DEG*psMetadataLookupF32 (&status, row, "RA");
     73        ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
     74        ref->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // XXX ELIXIR uses millimags, PANSTARRS uses mags
    7575        # else
    76         object->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
    77         object->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
    78         object->Mag      = psMetadataLookupF32 (&status, row, "MAG");
     76        ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
     77        ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
     78        ref->Mag      = psMetadataLookupF32 (&status, row, "MAG");
    7979        # endif
    8080
    81         psArrayAdd (objects, 100, object);
    82         psFree (object);
     81        psArrayAdd (refs, 100, ref);
     82        psFree (ref);
    8383    }
    8484    psFree (header);
    … …  
    8686
    8787    psTrace (__func__, 3, "loaded %d reference stars from (%10.6f,%10.6f) - (%10.6f,%10.6f)\n",
    88              objects->n, RAmin, DECmin, RAmax, DECmax);
    89     return objects;
     88             refs->n, RAmin, DECmin, RAmax, DECmax);
     89    return refs;
    9090}
    9191

  • trunk/psastro/src/psastroMosaicAstrom.c

    r7154 r7332  
    2020    }
    2121
    22     double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE");
    23     if (!status) plateScale = 1.0;
     22    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.DISTORTION.ORDER");
     23    if (!status) order = 3.0;
    2424
    2525    pmFPA *fpa = input->fpa;
    2626
    27     psastroMosaicGetRefstars (fpa, refs, recipe);
     27    // XXX test point
     28    if (0) {
     29        pmChip *chip = fpa->chips->data[10];
     30        pmCell *cell = chip->cells->data[0];
     31        pmReadout *readout = cell->readouts->data[0];
     32        psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     33        pmAstromObj *ref = refstars->data[0];
     34        fprintf (stderr, "test 1: %f %f  %f %f  %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d);
     35    }
     36
    2837    psastroMosaicSetMatch (fpa, recipe);
     38
     39    // XXX no input distortion model yet; existing fpa distortion is identity: replace it
     40    psFree (fpa->toTangentPlane);
     41    psFree (fpa->fromTangentPlane);
     42    fpa->toTangentPlane   = psPlaneDistortIdentity (order);
     43    fpa->fromTangentPlane = psPlaneDistortIdentity (order);
    2944
    3045    grads = psastroMosaicGetGrads (fpa, recipe);
    … …  
    3247    // fit the measured gradients with the telescope distortion model (3rd order polynomial)
    3348    pmAstromFitDistortion (fpa, grads, recipe);
     49    psastroMosaicRescaleChips (fpa);
    3450
    3551    // apply the new distortion terms up and down
    … …  
    5066    // re-fit the chips with higher-order fits
    5167    psastroMosaicChipAstrom (fpa, recipe);
    52     // psastroMosaicHeaders (fpa);
     68    psastroMosaicHeaders (config);
    5369
    5470    return true;

  • trunk/psastro/src/psastroMosaicChipAstrom.c

    r7154 r7332  
    2323
    2424                // select the raw objects for this readout
    25                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     25                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    2626                if (rawstars == NULL) { continue; }
    2727

  • trunk/psastro/src/psastroMosaicGetGrads.c

    r7154 r7332  
    2525
    2626                // select the raw objects for this readout
    27                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     27                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    2828                if (rawstars == NULL) { continue; }
    2929

  • trunk/psastro/src/psastroMosaicGetRefstars.c

    r7154 r7332  
    11# include "psastro.h"
    22
    3 bool psastroMosaicGetRefstars (pmFPA *fpa, psArray *refs, psMetadata *recipe) {
     3bool psastroMosaicGetRefstars (pmConfig *config, psArray *refs) {
    44
     5    bool status;
    56    pmChip *chip = NULL;
    67    pmCell *cell = NULL;
    78    pmReadout *readout = NULL;
     9
     10    // select the current recipe
     11    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
     12    if (!recipe) {
     13        psErrorStackPrint(stderr, "Can't find PSASTRO recipe!\n");
     14        exit(EXIT_FAILURE);
     15    }
     16
     17    // select the input data sources
     18    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     19    if (!input) {
     20        psErrorStackPrint(stderr, "Can't find input data!\n");
     21        exit(EXIT_FAILURE);
     22    }
     23
     24    float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA");
     25    if (!status) fieldExtra = 0.0;
     26
    827    pmFPAview *view = pmFPAviewAlloc (0);
     28    pmFPA *fpa = input->fpa;
    929
    1030    // this loop selects the matched stars for all chips
    … …  
    2242                if (! readout->data_exists) { continue; }
    2343
    24                 // select the raw objects for this readout
    25                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
    26                 if (rawstars == NULL) { continue; }
     44                // read WCS data from the corresponding header
     45                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    2746
     47                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
     48                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");
     49
     50                float minX = -fieldExtra*Nx;
     51                float maxX = (1+fieldExtra)*Nx;
     52                float minY = -fieldExtra*Ny;
     53                float maxY = (1+fieldExtra)*Ny;
     54               
    2855                // the refstars is a subset within range of this chip
    2956                psArray *refstars = psArrayAlloc (100);
     57
     58                int Nst = 0;
    3059
    3160                // select the reference objects within range of this readout
    3261                // project the reference objects to this chip
    3362                for (int i = 0; i < refs->n; i++) {
    34                     pmAstromObj *ref = refs->data[i];
     63                    pmAstromObj *ref = pmAstromObjCopy(refs->data[i]);
    3564       
    3665                    p_psProject (ref->TP, ref->sky, fpa->projection);
    … …  
    3867                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    3968
    40                     // XXX what is the X,Y range of the selected chip?
    41                     // XXX for the moment, force these to be good for Megacam
    42                     if (ref->chip->x <  -10) continue;
    43                     if (ref->chip->x > 2060) continue;
    44                     if (ref->chip->y <  -10) continue;
    45                     if (ref->chip->y > 4610) continue;
     69                    // limit the X,Y range of the refs to the selected chip
     70                    if (ref->chip->x < minX) continue;
     71                    if (ref->chip->x > maxX) continue;
     72                    if (ref->chip->y < minY) continue;
     73                    if (ref->chip->y > maxY) continue;
    4674                   
     75                    if (Nst < 0) {
     76                        fprintf (stderr, "dn: %f,%f <- %f,%f <- %f,%f <- %f,%f\n",
     77                                 ref->chip->x, ref->chip->y,
     78                                 ref->FP->x, ref->FP->y,
     79                                 ref->TP->x, ref->TP->y,
     80                                 ref->sky->r, ref->sky->d);
     81
     82                        psSphere *sk = psSphereAlloc();
     83                        psPlane *fp = psPlaneAlloc();
     84                        psPlane *tp = psPlaneAlloc();
     85
     86                        psPlaneTransformApply (fp, chip->toFPA, ref->chip);
     87                        psPlaneDistortApply (tp, fpa->toTangentPlane, fp, 0.0, 0.0);
     88                        p_psDeproject (sk, tp, fpa->projection);
     89
     90                        fprintf (stderr, "up: %f,%f -> %f,%f -> %f,%f -> %f,%f\n",
     91                                 ref->chip->x, ref->chip->y,
     92                                 fp->x, fp->y,
     93                                 tp->x, tp->y,
     94                                 sk->r, sk->d);
     95                    }
     96                    Nst ++;
     97
    4798                    psArrayAdd (refstars, 100, ref);
    4899                }
     100                psTrace (__func__, 4, "Added %d refstars\n", refstars->n);
     101
    49102                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
     103
     104                // XXX test point
     105                if (0) {
     106                    pmChip *chip1 = fpa->chips->data[10];
     107                    pmCell *cell1 = chip1->cells->data[0];
     108                    pmReadout *readout1 = cell1->readouts->data[0];
     109                    psArray *refstars1 = psMetadataLookupPtr (NULL, readout1->analysis, "PSASTRO.REFSTARS");
     110                    pmAstromObj *ref1 = refstars1->data[0];
     111                    fprintf (stderr, "test 3: %f %f  %f %f  %f %f\n", ref1->chip->x, ref1->chip->y, ref1->FP->x, ref1->FP->y, ref1->sky->r, ref1->sky->d);
     112                }
     113
     114                psFree (refstars);
    50115            }
    51116        }
    52117    }
     118
     119    psFree (view);
    53120    return true;
    54121}

  • trunk/psastro/src/psastroMosaicSetAstrom.c

    r7154 r7332  
    2323
    2424                // select the raw objects for this readout
    25                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     25                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    2626                if (rawstars == NULL) { continue; }
    2727

  • trunk/psastro/src/psastroMosaicSetMatch.c

    r7154 r7332  
    77    pmReadout *readout = NULL;
    88    pmFPAview *view = pmFPAviewAlloc (0);
     9
     10    FILE *f;
     11    char name[128];
     12
     13    FILE *g1 = fopen ("raw.ps.dat", "w");
     14    FILE *g2 = fopen ("ref.ps.dat", "w");
    915
    1016    // this loop selects the matched stars for all chips
    … …  
    2329
    2430                // select the raw objects for this readout
    25                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.OBJECTS");
     31                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    2632                if (rawstars == NULL) { continue; }
    2733
    … …  
    2935                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    3036                if (refstars == NULL) { continue; }
     37                psTrace (__func__, 4, "Trying %d refstars\n", refstars->n);
    3138
    3239                // use small radius to match stars (assume starting astrometry is good)
    3340                // XXX should this take a (double radius)?
    34                 psArray *matches = pmAstromRadiusMatch (rawstars, refstars, recipe);
     41                psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, recipe);
     42                psTrace (__func__, 4, "Matched %d refstars\n", matches->n);
     43
     44                sprintf (name, "raw.%02d.dat", view->chip);
     45                f = fopen (name, "w");
     46                for (int i = 0; i < rawstars->n; i++) {
     47                    pmAstromObj *raw = rawstars->data[i];
     48                    fprintf (f, "%f %f  %f %f  %f %f\n", raw->chip->x, raw->chip->y, raw->FP->x, raw->FP->y, raw->sky->r, raw->sky->d);
     49                }
     50                fclose (f);
     51
     52                sprintf (name, "ref.%02d.dat", view->chip);
     53                f = fopen (name, "w");
     54                for (int i = 0; i < refstars->n; i++) {
     55                    pmAstromObj *ref = refstars->data[i];
     56                    fprintf (f, "%f %f  %f %f  %f %f\n", ref->chip->x, ref->chip->y, ref->FP->x, ref->FP->y, ref->sky->r, ref->sky->d);
     57                }
     58                fclose (f);
     59
     60                for (int i = 0; i < matches->n; i++) {
     61                    pmAstromMatch *match = matches->data[i];
     62
     63                    pmAstromObj *raw = rawstars->data[match->raw];
     64                    fprintf (g1, "%d %f %f  %f %f  %f %f  %f %f\n", i,
     65                             DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d,
     66                             raw->TP->x, raw->TP->y,
     67                             raw->FP->x, raw->FP->y,
     68                             raw->chip->x, raw->chip->y);
     69
     70                    pmAstromObj *ref = refstars->data[match->ref];
     71                    fprintf (g2, "%d %f %f  %f %f  %f %f  %f %f\n", i,
     72                             DEG_RAD*ref->sky->r, DEG_RAD*ref->sky->d,
     73                             ref->TP->x, ref->TP->y,
     74                             ref->FP->x, ref->FP->y,
     75                             ref->chip->x, ref->chip->y);
     76                }
    3577
    3678                // XXX drop the old one
    37                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY, "astrometry matches", matches);
     79                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches);
    3880            }
    3981        }
    4082    }
     83    fclose (g1);
     84    fclose (g2);
    4185    return true;
    4286}

  • trunk/psastro/src/psastroWCS.c

    r7014 r7332  
    33// interpret header WCS (only handles traditional WCS for the moment)
    44// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
    5 bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) {
    6 
    7     psPlaneTransform *toFPA;
    8     psPlaneDistort *toTPA;
    9     psProjection *toSky;
     5bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic) {
     6
    107    psProjectionType type;
    118    bool status;
    … …  
    1310    float pc1_1, pc1_2, pc2_1, pc2_2;
    1411
    15     // *** interpret header data, convert to crval(i), etc
     12    // interpret header data, convert to crval(i), etc
    1613    char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
    1714    if (!status) {
    … …  
    8885got_matrix:
    8986
    90     // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
    91     //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
    92     //           scale = scale(i)/scale(0) (i == chip #)
    93     //           project crval1(0),crval2(0 using projection
    94     toFPA = psPlaneTransformAlloc (1, 1);
     87    /*****
     88
     89    For mosaic astrometry, we need to have a starting set of projection terms in which the
     90    chip-to-FPA terms result in a fixed physical unit on the focal plane (eg, pixels or
     91    microns).  This set of projections, coupled with an identity toTPA (ie, no distortion) will
     92    result in substantial errors between the observed and predicted star positions on the focal
     93    plane: this is the measurement of the optical distortion in the camera.  At the same time,
     94    we need to carry around the transformations which allow us to make an accurate calculation
     95    of the position of the stars based on the input (per-chip) astrometry.  These
     96    transformations will allow us to match the raw and ref stars robustly.  To convert the
     97    per-chip astrometry (which may have been calculated with a different plate scale for each
     98    chip) to a collection of astrometry terms for chips in a single mosaic, we need to adjust
     99    the chip-to-FPA scaling (eg, pc11) to match the variations in the effective plate scale for
     100    each chip (eg, cdelt1).  Thus, we need to carry around both the
     101
     102    *****/
     103
     104    {
     105        double rX = 1.0;
     106        double rY = 1.0;
     107
     108        // XXX free an existing toFPA?
     109        psPlaneTransform *toFPA = psPlaneTransformAlloc (1, 1);
    95110   
    96     double cdelt = hypot (cdelt1, cdelt2) / plateScale;  // degrees / micron (eg, in fact, whatever unit user chooses for focal plane)
    97     cdelt1 /= cdelt;
    98     cdelt2 /= cdelt;
    99 
    100     toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2);
    101     toFPA->x->coeff[1][0] = +(pc1_1*cdelt1);
    102     toFPA->x->coeff[0][1] = +(pc1_2*cdelt2);
    103     toFPA->x->mask[1][1]  = 1;
    104 
    105     toFPA->y->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2);
    106     toFPA->y->coeff[1][0] = +(pc2_1*cdelt1);
    107     toFPA->y->coeff[0][1] = +(pc2_2*cdelt2);
    108     toFPA->y->mask[1][1]  = 1;
    109     *toFPAout = toFPA;
    110 
    111     if (*toSkyOut != NULL) {
    112         if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
    113 
    114         psSphere sky;
    115         psPlane tpa;
    116 
    117         sky.r = crval1*RAD_DEG;
    118         sky.d = crval2*RAD_DEG;
    119         p_psProject (&tpa, &sky, *toSkyOut);
     111        // basic transformation from chip to FPA
     112        toFPA->x->coeff[0][0] = -(pc1_1*crpix1 + pc1_2*crpix2);
     113        toFPA->x->coeff[1][0] = pc1_1;
     114        toFPA->x->coeff[0][1] = pc1_2;
     115        toFPA->x->mask[1][1]  = 1;
     116
     117        toFPA->y->coeff[0][0] = -(pc2_1*crpix1 + pc2_2*crpix2);
     118        toFPA->y->coeff[1][0] = pc2_1;
     119        toFPA->y->coeff[0][1] = pc2_2;
     120        toFPA->y->mask[1][1]  = 1;
     121
     122        // projection from TPA to SKY
     123        psProjection *toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     124
     125        if (fpa->toSky == NULL) {
     126            // XXX for now, use the identity for TPA <--> FPA
     127            fpa->toTangentPlane = psPlaneDistortIdentity (1);
     128            fpa->fromTangentPlane = psPlaneDistortIdentity (1);
     129            fpa->toSky = toSky;
     130        } else {
     131            if (fpa->toTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
     132            if (fpa->fromTangentPlane == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
     133
     134            // adjust for common toSky for mosaic:
     135            // ignore the TPA since toTPA is identity
     136            // find the FPA coordinate of 0,0 for this chip.
     137            psPlane *fp = psPlaneAlloc();
     138            psPlane *chip = psPlaneAlloc();
     139            psSphere *sky = psSphereAlloc();
     140            chip->x = chip->y = 0;
    120141       
    121         // XXX for the moment, assume toTPA is the identity transformation
    122         toFPA->x->coeff[0][0] = tpa.x;
    123         toFPA->y->coeff[0][0] = tpa.y;
    124     } else {
    125         // XXX for now, use the identity for TPA <--> FPA
    126         toTPA = psPlaneDistortIdentity ();
    127 
    128         // center of projection is (0,0) coordinate of TPA
    129         toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type);
    130 
    131         *toTPAout = toTPA;
    132         *toSkyOut = toSky;
     142            psPlaneTransformApply (fp, toFPA, chip); // find the focal-plane coordinate of this chip's 0,0 coordinate
     143            p_psDeproject (sky, fp, toSky); // find the RA,DEC coord of the focal-plane coordinate
     144            p_psProject (fp, sky, fpa->toSky); // find the focal-plane coord of this RA,DEC coord using the ref chip projection
     145
     146            toFPA->x->coeff[0][0] = fp->x;
     147            toFPA->y->coeff[0][0] = fp->y;
     148
     149            rX = toSky->Xs/fpa->toSky->Xs;
     150            rY = toSky->Ys/fpa->toSky->Ys;
     151
     152            // correct to common plate scale
     153            toFPA->x->coeff[1][0] *= rX;
     154            toFPA->x->coeff[0][1] *= rX;
     155            toFPA->y->coeff[1][0] *= rY;
     156            toFPA->y->coeff[0][1] *= rY;
     157
     158            psFree (fp);
     159            psFree (sky);
     160            psFree (chip);
     161            psFree (toSky);
     162        }
     163
     164        chip->toFPA = toFPA;
     165        chip->fromFPA = p_psPlaneTransformLinearInvert(toFPA);
     166
     167        // remove the correction to the common plate scale
     168        if (isMosaic) {
     169            chip->toFPA->x->coeff[1][0] /= rX;
     170            chip->toFPA->x->coeff[0][1] /= rX;
     171            chip->toFPA->y->coeff[1][0] /= rY;
     172            chip->toFPA->y->coeff[0][1] /= rY;
     173        }
     174
     175        fprintf (stderr, "toFPA: %f %f  (%f,%f),(%f,%f)\n",
     176                 chip->toFPA->x->coeff[0][0], chip->toFPA->y->coeff[0][0],
     177                 chip->toFPA->x->coeff[1][0], chip->toFPA->x->coeff[0][1],
     178                 chip->toFPA->y->coeff[1][0], chip->toFPA->y->coeff[0][1]);
     179
     180        fprintf (stderr, "frFPA: %f %f  (%f,%f),(%f,%f)\n",
     181                 chip->fromFPA->x->coeff[0][0], chip->fromFPA->y->coeff[0][0],
     182                 chip->fromFPA->x->coeff[1][0], chip->fromFPA->x->coeff[0][1],
     183                 chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
     184
     185        fprintf (stderr, "field: %f,%f  %f,%f\n",
     186                 DEG_RAD*fpa->toSky->R, DEG_RAD*fpa->toSky->D,
     187                 3600*DEG_RAD*fpa->toSky->Xs, 3600*DEG_RAD*fpa->toSky->Ys);
     188
    133189    }
    134190    return true;
    135191}
    136 
    137192
    138193// convert toFPA / toSky components to traditional WCS
    … …  
    169224    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
    170225    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
    171     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
    172     psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
    173 
    174     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
    175     psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
    176     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
    177     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
     226    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
     227    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     228
     229    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
     230    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
     231    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
     232    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
    178233
    179234    // alternative representations use
    … …  
    184239}
    185240
    186 psPlaneDistort *psPlaneDistortIdentity () {
     241// convert toFPA / toSky components to traditional WCS
     242// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     243bool pmAstromWriteBilevelChip (psPlaneTransform *toFPA, psMetadata *header, double plateScale) {
     244
     245    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---WRP");
     246    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--WRP");
     247
     248    // XXX not really right: needs to deal with non-identity toTP coeffs & plateScale
     249    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     250    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
     251    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", 0.0);
     252    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", 0.0);
     253    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", 1.0);
     254    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", 1.0);
     255
     256    if (toFPA->x->nX != toFPA->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
     257    if (toFPA->x->nX != toFPA->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
     258    if (toFPA->x->nX != toFPA->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
     259
     260    switch (toFPA->x->nX) {
     261      case 3:
     262        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toFPA->x->coeff[3][0]);
     263        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toFPA->x->coeff[2][1]);
     264        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toFPA->x->coeff[1][2]);
     265        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toFPA->x->coeff[0][3]);
     266
     267        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toFPA->y->coeff[3][0]);
     268        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toFPA->y->coeff[2][1]);
     269        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toFPA->y->coeff[1][2]);
     270        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toFPA->y->coeff[0][3]);
     271
     272      case 2:
     273        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toFPA->x->coeff[2][0]);
     274        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toFPA->x->coeff[1][1]);
     275        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toFPA->x->coeff[0][2]);
     276
     277        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toFPA->y->coeff[2][0]);
     278        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toFPA->y->coeff[1][1]);
     279        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toFPA->y->coeff[0][2]);
     280
     281      case 1:
     282        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]/plateScale);
     283        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]/plateScale);
     284        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
     285        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
     286        break;
     287
     288      case 0:
     289        psAbort ("psastro", "invalid tangent plane order");
     290    }
     291    return true;
     292}
     293
     294// convert toFPA / toSky components to traditional WCS
     295// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
     296psMetadata *pmAstromWriteBilevelMosaic (psProjection *toSky, psPlaneDistort *toTP, double plateScale) {
     297
     298    psMetadata *header = psMetadataAlloc ();
     299
     300    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---DIS");
     301    psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--DIS");
     302
     303    // XXX need to handle the plateScale??
     304    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toSky->R*DEG_RAD);
     305    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toSky->D*DEG_RAD);
     306    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", 0.0);
     307    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", 0.0);
     308    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD*plateScale);
     309    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD*plateScale);
     310
     311    if (toTP->x->nX != toTP->x->nY) psAbort ("psastro", "mis-matched tangent plane orders (1)");
     312    if (toTP->x->nX != toTP->y->nX) psAbort ("psastro", "mis-matched tangent plane orders (2)");
     313    if (toTP->x->nX != toTP->y->nY) psAbort ("psastro", "mis-matched tangent plane orders (3)");
     314
     315    switch (toTP->x->nX) {
     316      case 3:
     317        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X3Y0", PS_META_REPLACE, "", toTP->x->coeff[3][0][0][0]);
     318        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y1", PS_META_REPLACE, "", toTP->x->coeff[2][1][0][0]);
     319        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y2", PS_META_REPLACE, "", toTP->x->coeff[1][2][0][0]);
     320        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y3", PS_META_REPLACE, "", toTP->x->coeff[0][3][0][0]);
     321
     322        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X3Y0", PS_META_REPLACE, "", toTP->y->coeff[3][0][0][0]);
     323        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y1", PS_META_REPLACE, "", toTP->y->coeff[2][1][0][0]);
     324        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y2", PS_META_REPLACE, "", toTP->y->coeff[1][2][0][0]);
     325        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y3", PS_META_REPLACE, "", toTP->y->coeff[0][3][0][0]);
     326
     327      case 2:
     328        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X2Y0", PS_META_REPLACE, "", toTP->x->coeff[2][0][0][0]);
     329        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X1Y1", PS_META_REPLACE, "", toTP->x->coeff[1][1][0][0]);
     330        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA1X0Y2", PS_META_REPLACE, "", toTP->x->coeff[0][2][0][0]);
     331
     332        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X2Y0", PS_META_REPLACE, "", toTP->y->coeff[2][0][0][0]);
     333        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X1Y1", PS_META_REPLACE, "", toTP->y->coeff[1][1][0][0]);
     334        psMetadataAddF32 (header, PS_LIST_TAIL, "PCA2X0Y2", PS_META_REPLACE, "", toTP->y->coeff[0][2][0][0]);
     335
     336      case 1:
     337        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toTP->x->coeff[1][0][0][0]);
     338        psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toTP->x->coeff[0][1][0][0]);
     339        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toTP->y->coeff[1][0][0][0]);
     340        psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toTP->y->coeff[0][1][0][0]);
     341        break;
     342
     343      case 0:
     344        psAbort ("psastro", "invalid tangent plane order");
     345    }
     346    return header;
     347}
     348
     349// XXX for the moment, we mimic the limited Elixir mosastro orders
     350psPlaneDistort *psPlaneDistortIdentity (int order) {
    187351
    188352    psPlaneDistort *distort;
    189353
    190     distort = psPlaneDistortAlloc (1, 1, 0, 0);
    191     distort->x->coeff[0][0][0][0] = 0;
     354    if (order < 1) psAbort ("psastro", "invalid order");
     355    if (order > 3) psAbort ("psastro", "invalid order");
     356   
     357    // all coeffs and masks initially set to 0
     358    distort = psPlaneDistortAlloc (order, order, 0, 0);
     359
     360    for (int i = 0; i <= order; i++) {
     361        for (int j = 0; j <= order; j++) {
     362            if (i + j > order) {
     363                distort->x->mask [1][1][0][0] = 1;
     364                distort->y->mask [1][1][0][0] = 1;
     365            }
     366        }
     367    }
    192368    distort->x->coeff[1][0][0][0] = 1;
    193     distort->x->coeff[0][1][0][0] = 0;
    194     distort->x->mask [1][1][0][0] = 1;
    195 
    196     distort->y->coeff[0][0][0][0] = 0;
    197     distort->y->coeff[1][0][0][0] = 0;
    198369    distort->y->coeff[0][1][0][0] = 1;
    199     distort->y->mask [1][1][0][0] = 1;
    200370
    201371    return distort;
Note: See TracChangeset for help on using the changeset viewer.