IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 10880


Ignore:
Timestamp:
Jan 2, 2007, 3:26:52 PM (19 years ago)
Author:
eugene
Message:

fixed mosaic analysis; cleaned up test code (some in traces, some removed), general output cleanups

Location:
trunk/psastro/src
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastro.h

    r10864 r10880  
    5656
    5757// mosaic fitting functions
    58 psArray *psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe);
     58bool psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe);
    5959bool psastroMosaicCommonScale (pmFPA *fpa, psMetadata *recipe);
    6060bool psastroMosaicAstrom (pmConfig *config, psArray *refs);
    61 bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration);
     61bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, bool nonlinear);
    6262bool psastroMosaicSetMatch (pmFPA *fpa, psMetadata *recipe, int iteration);
    6363bool psastroMosaicSetAstrom (pmFPA *fpa);
    6464bool psastroMosaicHeaders (pmConfig *config);
    6565bool psastroMosaicRescaleChips (pmFPA *fpa);
    66 bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration);
     66bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, bool nonlinear);
    6767
    6868// Return version strings.
    … …  
    7777bool psastroDumpRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip);
    7878bool psastroDumpRefstars (psArray *refstars);
     79bool psastroDumpMatches (pmFPA *fpa, char *filename);
    7980
    8081bool psastroMosaicSetAstrom_tmp (pmFPA *fpa);

  • trunk/psastro/src/psastroChipAstrom.c

    r10855 r10880  
    44bool psastroChipAstrom (pmConfig *config, psArray *refs) {
    55
    6     bool status;
    76    pmChip *chip = NULL;
    87    pmCell *cell = NULL;
    … …  
    2221        return false;
    2322    }
    24 
    25     double plateScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLATE.SCALE");
    26     if (!status) plateScale = 1.0;
    2723
    2824    pmFPAview *view = pmFPAviewAlloc (0);

  • trunk/psastro/src/psastroChooseRefstars.c

    r10855 r10880  
    5656                psArray *refstars = psArrayAllocEmpty (100);
    5757
    58                 FILE *f = fopen ("refstars.dat", "w");
    59 
    6058                // select the reference objects within range of this readout
    6159                // project the reference objects to this chip
    … …  
    6664                    psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
    6765                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    68 
    69                     fprintf (f, "%d  %f %f  %f %f  %f %f  %f %f\n", i,
    70                              ref->sky->r, ref->sky->d,
    71                              ref->TP->x, ref->TP->y,
    72                              ref->FP->x, ref->FP->y,
    73                              ref->chip->x, ref->chip->y);
    7466
    7567                    // limit the X,Y range of the refs to the selected chip
    … …  
    8375                    psFree (ref);
    8476                }
    85                 fclose (f);
    8677                psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
    8778

  • trunk/psastro/src/psastroDemoDump.c

    r10864 r10880  
    5858    return true;
    5959}
     60
     61bool psastroDumpMatches (pmFPA *fpa, char *filename) {
     62
     63    pmChip *chip = NULL;
     64    pmCell *cell = NULL;
     65    pmReadout *readout = NULL;
     66    pmFPAview *view = pmFPAviewAlloc (0);
     67
     68    FILE *f = fopen (filename, "w");
     69
     70    // this loop selects the matched stars for all chips
     71    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
     72        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     73        if (!chip->process || !chip->file_exists) continue;
     74       
     75        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     76            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     77            if (!cell->process || !cell->file_exists) continue;
     78
     79            // process each of the readouts
     80            // XXX there can only be one readout per chip, right?
     81            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     82                if (! readout->data_exists) continue;
     83
     84                // select the raw objects for this readout
     85                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     86                if (rawstars == NULL) continue;
     87
     88                // select the raw objects for this readout
     89                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     90                if (refstars == NULL) continue;
     91                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
     92
     93                psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH");
     94                if (matches == NULL) continue;
     95
     96                for (int i = 0; i < matches->n; i++) {
     97                    pmAstromMatch *match = matches->data[i];
     98
     99                    pmAstromObj *raw = rawstars->data[match->raw];
     100                    fprintf (f, "%f %f  %f %f  %f %f  %f %f  %f   |   ", 
     101                             DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d,
     102                             raw->TP->x, raw->TP->y,
     103                             raw->FP->x, raw->FP->y,
     104                             raw->chip->x, raw->chip->y, raw->Mag);
     105
     106                    pmAstromObj *ref = refstars->data[match->ref];
     107                    fprintf (f, "%f %f  %f %f  %f %f  %f %f  %f\n",
     108                             DEG_RAD*ref->sky->r, DEG_RAD*ref->sky->d,
     109                             ref->TP->x, ref->TP->y,
     110                             ref->FP->x, ref->FP->y,
     111                             ref->chip->x, ref->chip->y, ref->Mag);
     112                }
     113            }
     114        }
     115    }
     116    fclose (f);
     117    psFree (view);
     118    return true;
     119}

  • trunk/psastro/src/psastroLuminosityFunction.c

    r10855 r10880  
    6868  // exclude bins with no stars
    6969  // exclude points after the peak with N < 0.8*peak
    70   FILE *f = fopen ("psastro.dat", "w");
    7170  int n = 0;
    7271  for (int i = 0; i < nMags->n; i++) {
    … …  
    7574    lnMag->data.F32[n] = log10 (nMags->data.F32[i]);
    7675    Mag->data.F32[n] = mMin + (i + 0.5)*dMag;
    77     fprintf (f, "%d  %f %.0f %f\n", n, Mag->data.F32[n], nMags->data.F32[i], lnMag->data.F32[n]);
    7876    n++;
    7977  }
    80   fclose (f);
    8178  lnMag->n = n;
    8279  Mag->n = n;

  • trunk/psastro/src/psastroMosaicAstrom.c

    r10864 r10880  
    44// XXX require this fpa to have multiple chip extensions and a PHU?
    55bool psastroMosaicAstrom (pmConfig *config, psArray *refs) {
    6 
    7     bool status;
    8     psArray *gradients = NULL;
    96
    107    // select the current recipe
    … …  
    2926
    3027    // given the existing per-chip astrometry, determine matches between raw and ref stars
     28    // is this needed? yes, if we didn't do SingleChip astrometry first
    3129    psastroMosaicSetMatch (fpa, recipe, 0);
     30    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.0.dat"); }
    3231
    3332    // fitted chips will follow the local plate-scale, hiding the distortion
    … …  
    3534    // then recalculate raw and ref positions
    3635    psastroMosaicCommonScale (fpa, recipe);
    37     psastroMosaicSetAstrom (fpa);
     36    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.1.dat"); }
    3837
    39     gradients = psastroMosaicGradients (fpa, recipe);
     38    // fit the distortion by fitting its gradient
     39    // apply the new distortion terms up and down
     40    // refit the per-chip terms with linear fits only
     41    psastroMosaicGradients (fpa, recipe);
     42    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.2.dat"); }
    4043
    41     // allocate mosaic-level polynomial transformation and set masks needed by DVO
    42     int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
    43     if (!status) {
    44         psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
    45         return false;
    46     }
    47     psFree (fpa->toTPA);
    48     fpa->toTPA = psPlaneTransformAlloc (order, order);
    49     for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
    50         for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
    51             if (i + j > order) {
    52                 fpa->toTPA->x->mask[i][j] = 1;
    53                 fpa->toTPA->y->mask[i][j] = 1;
    54             }
    55         }
    56     }
     44    psastroMosaicChipAstrom (fpa, recipe, false);
     45    if (psTraceGetLevel("psastro.dump") > 0) { psastroDumpMatches (fpa, "match.3.dat"); }
    5746
    58     // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA)
    59     pmAstromFitDistortion (fpa, gradients, recipe);
    60     psFree (gradients);
    61 
    62     // apply the new distortion terms up and down
    6347    // re-perform the match with a slightly tighter circle
    64     // XXX modify match radius
    65     // XXX set chip.order to 1
    66     psastroMosaicSetAstrom_tmp (fpa);
    6748    psastroMosaicSetMatch (fpa, recipe, 1);
    68     psastroMosaicChipAstrom (fpa, recipe, 1);
    6949
    7050    // do a second pass on the distortion with improved chip positions and rotations
    71     // XXX do we need to iterate with this step until the distortions don't change?
    7251    psastroMosaicCommonScale (fpa, recipe);
    73     gradients = psastroMosaicGradients (fpa, recipe);
    74     pmAstromFitDistortion (fpa, gradients, recipe);
    75     psFree (gradients);
     52    psastroMosaicGradients (fpa, recipe);
     53    psastroMosaicChipAstrom (fpa, recipe, false);
     54    psastroMosaicSetMatch (fpa, recipe, 1);
    7655   
    7756    // now fit the chips under the common distortion with higher-order terms
    78     // XXX modify match radius
    79     // XXX set chip.order to NORDER
    80     psastroMosaicSetAstrom (fpa);
    81     psastroMosaicSetMatch (fpa, recipe, 2);
    82     psastroMosaicChipAstrom (fpa, recipe, 2);
     57    // first, re-perform the match with a slightly tighter circle
     58    psastroMosaicChipAstrom (fpa, recipe, true);
    8359
    8460    // save WCS and analysis metadata in update header
    85     // XXX need to add global summary statistics
    8661    psMetadata *updates = psMetadataAlloc();
    8762    pmAstromWriteBilevelMosaic (updates, fpa, NONLIN_TOL);
    88     psMetadataAdd (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
     63    psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    8964    psFree (updates);
    9065
    9166    // update the headers based on the results
     67    // XXX need to add global summary statistics
    9268    // psastroMosaicHeaders (config);
    9369

  • trunk/psastro/src/psastroMosaicChipAstrom.c

    r10830 r10880  
    22# define NONLIN_TOL 0.001 /* tolerance in pixels */
    33
    4 bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, int iteration) {
     4bool psastroMosaicChipAstrom (pmFPA *fpa, psMetadata *recipe, bool nonlinear) {
    55
    66    pmChip *chip = NULL;
    … …  
    2626                psMetadata *updates = psMetadataAlloc();
    2727
    28                 psastroMosaicOneChip (chip, readout, recipe, updates, iteration);
     28                psastroMosaicOneChip (chip, readout, recipe, updates, nonlinear);
    2929
     30                // create the header keywords to descripe the results
    3031                pmAstromWriteBilevelChip (updates, chip, NONLIN_TOL);
    31                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
     32                psMetadataAddMetadata (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    3233                psFree (updates);
    3334            }
    3435        }
    3536    }
     37    psFree (view);
    3638    return true;
    3739}
    38 
    39 /* the iteration value may be 1 or 2.
    40  */

  • trunk/psastro/src/psastroMosaicGradients.c

    r10830 r10880  
    11# include "psastro.h"
    22
    3 psArray *psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) {
     3bool psastroMosaicGradients (pmFPA *fpa, psMetadata *recipe) {
    44
     5    bool status;
    56    pmChip *chip = NULL;
    67    pmCell *cell = NULL;
    … …  
    3637
    3738                // measure the local gradients for this set of stars
    38                 // XXX update the function prototype to accept an incoming gradient structure to which the new elements are added
    3939                gradients = pmAstromMeasureGradients (gradients, rawstars, refstars, match, recipe);
    4040            }
    4141        }
    4242    }
    43     return (gradients);
     43
     44    // allocate mosaic-level polynomial transformation and set masks needed by DVO
     45    int order = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.ORDER");
     46    if (!status) {
     47        psError(PSASTRO_ERR_UNKNOWN, false, "failed to find single-chip fit order\n");
     48        return false;
     49    }
     50    psFree (fpa->toTPA);
     51    fpa->toTPA = psPlaneTransformAlloc (order, order);
     52    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
     53        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {
     54            if (i + j > order) {
     55                fpa->toTPA->x->mask[i][j] = 1;
     56                fpa->toTPA->y->mask[i][j] = 1;
     57            }
     58        }
     59    }
     60
     61    // fit the measured gradients with the telescope distortion model (polynomial order based on toTPA)
     62    pmAstromFitDistortion (fpa, gradients, recipe);
     63    psastroMosaicSetAstrom (fpa);
     64
     65    psFree (gradients);
     66    psFree (view);
     67    return true;
    4468}

  • trunk/psastro/src/psastroMosaicOneChip.c

    r10830 r10880  
    11# include "psastro.h"
    22
    3 bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) {
     3# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
     4  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
     5  if (!status) { \
     6   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
     7   return false; }
     8
     9bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, bool nonlinear) {
    410
    511    bool status;
    … …  
    2026    if (match == NULL) return false;
    2127
    22     if ((iteration < 1) || (iteration > 2)) {
    23         psError(PSASTRO_ERR_UNKNOWN, false, "invalid iteration number %d\n", iteration);
    24         return false;
     28    // correct radius to FP units (physical pixel scale in microns per pixel)
     29    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale");
     30
     31    // allowed limits for valid solutions
     32    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MOSAIC.MAX.ERROR", F32, "failed to find single-chip max allowed error\n");
     33    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MOSAIC.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
     34
     35    // set the order of the per-chip fit (higher order only if nonlinear == true)
     36    int order = 1;
     37    if (nonlinear) {
     38        // select the desired chip order
     39        REQUIRED_RECIPE_VALUE (order, "PSASTRO.MOSAIC.CHIP.ORDER", S32, "failed to find mosaic chip-level fit order\n");
     40
     41        // modify the order to correspond to the actual number of matched stars:
     42        if ((match->n < 11) && (order >= 3)) order = 2;
     43        if ((match->n <  7) && (order >= 2)) order = 1;
     44        if ((match->n <  4) && (order >= 1)) order = 0;
     45        if (order < 1) {
     46            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
     47            return false;
     48        }
    2549    }
    2650
    27     // correct to FP units (physical pixel scale in microns per pixel)
    28     double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    29     if (!status) {
    30         psError(PS_ERR_IO, false, "Failed to lookup pixel scale");
    31         return false;
    32     }
    33 
    34     // set the order of the per-chip fit (iteration 1: 1, iteration 2: from recipe)
    35     int order = 1;
    36     if (iteration == 2) {
    37         order = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.ORDER");
    38         if (!status) {
    39             psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic chip-level fit order\n");
    40             return false;
    41         }
    42     }
     51    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    4352    psFree (chip->toFPA);
    4453    chip->toFPA = psPlaneTransformAlloc (order, order);
    … …  
    6473    }
    6574
    66     // allowed limits for valid solutions
    67     // XXX have these values depend on the iteration?
    68     float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.MAX.ERROR");
    69     int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.MIN.NSTAR");
     75    // toSky converts from FPA & TPA units (microns) to sky units (radians)
     76    pmFPA *fpa = chip->parent;
     77    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    7078
    71     // astError is the average 1D scatter in pixels ('results' are in FPA units = microns)
    72     // XXX currently using supplied pixel scale: use toFPA->x->coeff[1][0], etc?
    73     float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
     79    // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
     80    float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
     81
     82    // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
     83    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    7484    int astNstar = results->yStats->clippedNvalues;
    7585
    … …  
    7787
    7888    // XXX should these result in errors or be handled another way?
    79     psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
     89    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
    8090    if (astError > maxError) {
    8191        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    … …  
    8898
    8999    // DVO expects NASTRO = 0 if we fail to find a solution
     100    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
     101    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    90102    if (validSolution) {
    91         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", astNstar);
    92         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar));
     103        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     104        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    93105    } else {
    94         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", 0);
    95         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0);
     106        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     107        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    96108    }
    97     psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "", astError);
    98109    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
     110
     111    // determine fromFPA transformation and apply new transformation to raw & ref stars
     112    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    99113
    100114    psFree (fitStats);

  • trunk/psastro/src/psastroMosaicSetAstrom.c

    r10864 r10880  
    4747        }
    4848    }
     49    psFree (view);
    4950    return true;
    5051}
    51 
    52 bool psastroMosaicSetAstrom_tmp (pmFPA *fpa) {
    53 
    54     pmChip *chip = NULL;
    55     pmCell *cell = NULL;
    56     pmReadout *readout = NULL;
    57     pmFPAview *view = pmFPAviewAlloc (0);
    58 
    59     FILE *f1 = fopen ("raw.tmp.dat", "w");
    60     FILE *f2 = fopen ("ref.tmp.dat", "w");
    61 
    62     // this loop selects the matched stars for all chips
    63     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
    64         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    65         if (!chip->process || !chip->file_exists) { continue; }
    66        
    67         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    68             psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    69             if (!cell->process || !cell->file_exists) { continue; }
    70 
    71             // process each of the readouts
    72             // XXX there can only be one readout per chip, right?
    73             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    74                 if (! readout->data_exists) { continue; }
    75 
    76                 // select the raw objects for this readout
    77                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
    78                 if (rawstars == NULL) { continue; }
    79 
    80                 for (int i = 0; i < rawstars->n; i++) {
    81                     pmAstromObj *raw = rawstars->data[i];
    82        
    83                     fprintf (f1, "%f %f  %f %f  %f %f  ", raw->sky->r, raw->sky->d, raw->TP->x, raw->TP->y, raw->FP->x, raw->FP->y);
    84 
    85                     psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    86                     psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
    87                     psDeproject (raw->sky, raw->TP, fpa->toSky);
    88 
    89                     fprintf (f1, "|  %f %f  %f %f  %f %f   %f\n", raw->sky->r, raw->sky->d, raw->TP->x, raw->TP->y, raw->FP->x, raw->FP->y, raw->Mag);
    90                 }
    91 
    92                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
    93                 if (refstars == NULL) { continue; }
    94 
    95                 for (int i = 0; i < refstars->n; i++) {
    96                     pmAstromObj *ref = refstars->data[i];
    97        
    98                     fprintf (f2, "%f %f  %f %f  %f %f  ", ref->sky->r, ref->sky->d, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y);
    99 
    100                     psProject (ref->TP, ref->sky, fpa->toSky);
    101                     psPlaneTransformApply (ref->FP, fpa->fromTPA, ref->TP);
    102                     psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    103 
    104                     fprintf (f2, "|  %f %f  %f %f  %f %f   %f\n", ref->sky->r, ref->sky->d, ref->TP->x, ref->TP->y, ref->FP->x, ref->FP->y, ref->Mag);
    105                 }
    106             }
    107         }
    108     }
    109     fclose (f1);
    110     fclose (f2);
    111     return true;
    112 }

  • trunk/psastro/src/psastroMosaicSetMatch.c

    r10830 r10880  
    88    pmFPAview *view = pmFPAviewAlloc (0);
    99    char radiusWord[64];
    10 
    11     FILE *g1 = fopen ("raw.ps.dat", "w");
    12     FILE *g2 = fopen ("ref.ps.dat", "w");
    1310
    1411    // use small radius to match stars (assume starting astrometry is good)
    … …  
    4744                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
    4845
    49                 for (int i = 0; i < matches->n; i++) {
    50                     pmAstromMatch *match = matches->data[i];
    51 
    52                     pmAstromObj *raw = rawstars->data[match->raw];
    53                     fprintf (g1, "%d %f %f  %f %f  %f %f  %f %f\n", i,
    54                              DEG_RAD*raw->sky->r, DEG_RAD*raw->sky->d,
    55                              raw->TP->x, raw->TP->y,
    56                              raw->FP->x, raw->FP->y,
    57                              raw->chip->x, raw->chip->y);
    58 
    59                     pmAstromObj *ref = refstars->data[match->ref];
    60                     fprintf (g2, "%d %f %f  %f %f  %f %f  %f %f\n", i,
    61                              DEG_RAD*ref->sky->r, DEG_RAD*ref->sky->d,
    62                              ref->TP->x, ref->TP->y,
    63                              ref->FP->x, ref->FP->y,
    64                              ref->chip->x, ref->chip->y);
    65                 }
    66 
    6746                // XXX drop the old one
    6847                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.MATCH", PS_DATA_ARRAY | PS_META_REPLACE, "astrometry matches", matches);
     48                psFree (matches);
    6949            }
    7050        }
    7151    }
    72     fclose (g1);
    73     fclose (g2);
     52    psFree (view);
    7453    return true;
    7554}

  • trunk/psastro/src/psastroOneChip.c

    r10864 r10880  
    2323    // allowed limits for valid solutions
    2424    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32, "failed to find single-chip max allowed error\n");
    25 
    2625    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
    2726
    … …  
    5352    }
    5453
    55     // create the output fit model, modify the order to correspond to the actual number of
    56     // matched stars:
     54    // modify the order to correspond to the actual number of matched stars:
    5755    if ((match->n < 11) && (order >= 3)) order = 2;
    5856    if ((match->n <  7) && (order >= 2)) order = 1;
    … …  
    6664    }
    6765
    68     // set masks appropriate to the Elixir DVO astrometry format
     66    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    6967    psFree (chip->toFPA);
    7068    chip->toFPA = psPlaneTransformAlloc (order, order);
    … …  
    125123    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    126124
    127     // write results
     125    // determine fromFPA transformation and apply new transformation to raw & ref stars
    128126    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    129127   

  • trunk/psastro/src/psastroUtils.c

    r10855 r10880  
    7373
    7474    }
     75    psastroMosaicSetAstrom (fpa);
    7576    return true;
    7677}
Note: See TracChangeset for help on using the changeset viewer.