IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 15891


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

finished psastroModel for boresite, created gpcModel for rough theoretical model, added useModel and fixChips options

Location:
trunk/psastro/src
Files:
5 added
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/Makefile.am

    r15887 r15891  
    33libpsastro_la_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
    44
    5 bin_PROGRAMS = psastro psastroModel
     5bin_PROGRAMS = psastro psastroModel gpcModel
    66
    77psastro_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
    … …  
    1212psastroModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
    1313psastroModel_LDADD = libpsastro.la
     14
     15gpcModel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSASTRO_CFLAGS)
     16gpcModel_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSASTRO_LIBS)
     17gpcModel_LDADD = libpsastro.la
    1418
    1519psastro_SOURCES = \
    … …  
    2933        psastroModelBoresite.c      \
    3034        psastroModelFitBoresite.c   \
     35        psastroModelAdjust.c        \
     36        psastroModelDataSave.c      \
    3137        psastroCleanup.c
     38
     39gpcModel_SOURCES = gpcModel.c
    3240
    3341## move DataSave to psastro_SOURCES?
    … …  
    4856        psastroLuminosityFunction.c \
    4957        psastroRefstarSubset.c      \
    50         psastroEnforceChips.c       \
     58        psastroFixChips.c           \
     59        psastroUseModel.c           \
    5160        psastroMosaicAstrom.c       \
    5261        psastroMosaicGradients.c    \

  • trunk/psastro/src/psastro.h

    r15562 r15891  
    8080int               psastroSortByMag (const void *a, const void *b);
    8181
    82 bool              psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe);
     82bool              psastroFixChips (pmConfig *config, psMetadata *recipe);
     83bool              psastroUseModel (pmConfig *config, psMetadata *recipe);
    8384
    84 psArray *psastroReadGetstarCatalog (psFits *fits);
    85 psArray *psastroReadGetstar_PS1_DEV_0 (psFits *fits);
     85psArray          *psastroReadGetstarCatalog (psFits *fits);
     86psArray          *psastroReadGetstar_PS1_DEV_0 (psFits *fits);
     87
     88bool              psastroAstromGuessSetChip (pmFPA *fpa, pmChip *chip, pmFPAview *view, double pixelScale, bool bilevelAstrometry);
     89bool              psastroAstromGuessSetFPA (pmFPA *fpa, bool *bilevelAstrometry);
    8690
    8791# endif /* PSASTRO_H */

  • trunk/psastro/src/psastroAnalysis.c

    r15562 r15891  
    55    bool status;
    66    int nStars;
     7
     8    // select the current recipe
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
     10    if (!recipe) {
     11        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe");
     12        return false;
     13    }
     14
     15    if (!psastroUseModel (config, recipe)) {
     16        psError (PSASTRO_ERR_UNKNOWN, false, "failed to set model astrometry\n");
     17        return false;
     18    }
    719
    820    // interpret the available initial astrometric information
    … …  
    2638    if (!psastroChooseRefstars (config, refs)) {
    2739        psError (PSASTRO_ERR_UNKNOWN, false, "failed to select reference data for chips\n");
    28         return false;
    29     }
    30 
    31     // select the current recipe
    32     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
    33     if (!recipe) {
    34         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe");
    3540        return false;
    3641    }

  • trunk/psastro/src/psastroArguments.c

    r15562 r15891  
    4343        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true);
    4444    }
     45    // no valid header WCS: supply from astrom model
     46    if ((N = psArgumentGet (argc, argv, "-use-astrommodel"))) {
     47        psArgumentRemove (N, &argc, argv);
     48        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.USE.MODEL", PS_META_REPLACE, "", true);
     49    }
    4550    // define the reference astrometry file
    46     status = pmConfigFileSetsMD (config->arguments, &argc, argv, "REF.ASTROM", "-refastrom", "-refastromlist");
     51    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist");
    4752    if (status) {
    4853        // if supplied, assume -fixchips is desired
    4954        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true);
    5055    }
     56
    5157
    5258    // apply mosastro mode?

  • trunk/psastro/src/psastroAstromGuess.c

    r15200 r15891  
    3232    }
    3333
     34    // have we already supplied the astrometry from the model?
     35    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
     36    if (!status) {
     37        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     38    }
     39
    3440    // select the input data sources
    3541    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    … …  
    4652    }
    4753
    48     pmFPAview *view = pmFPAviewAlloc (0);
    4954    pmFPA *fpa = input->fpa;
    5055
    5156    // load mosaic-level astrometry?
    5257    bool bilevelAstrometry = false;
    53     pmHDU *phu = pmFPAviewThisPHU (view, fpa);
    54     if (phu) {
    55       char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    56       if (ctype) {
    57         bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    58       }
     58    if (!useModel) {
     59        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    5960    }
    60     if (bilevelAstrometry) {
    61       pmAstromReadBilevelMosaic (fpa, phu->header);
    62     }
    6361
     62    pmFPAview *view = pmFPAviewAlloc (0);
    6463    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
    6564        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    6665        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    6766
    68         // read WCS data from the corresponding header
    69         pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    70         if (bilevelAstrometry) {
    71             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    72                 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    73                 continue;
    74             }
    75         } else {
    76             if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
    77                 psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    78                 continue;
    79             }
     67        if (!useModel) {
     68            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
    8069        }
    8170
    … …  
    160149   sky (ra, dec)
    161150*/
     151
     152bool psastroAstromGuessSetChip (pmFPA *fpa, pmChip *chip, pmFPAview *view, double pixelScale, bool bilevelAstrometry) {
     153
     154    // read WCS data from the corresponding header
     155    pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
     156    if (bilevelAstrometry) {
     157        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
     158            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     159            return false;
     160        }
     161    } else {
     162        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
     163            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     164            return false;
     165        }
     166    }
     167    return true;
     168}
     169
     170bool psastroAstromGuessSetFPA (pmFPA *fpa, bool *bilevelAstrometry) {
     171
     172    pmFPAview *view = pmFPAviewAlloc (0);
     173    pmHDU *phu = pmFPAviewThisPHU (view, fpa);
     174
     175    *bilevelAstrometry = false;
     176
     177    // load mosaic-level astrometry?
     178    if (phu) {
     179        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     180        if (ctype) {
     181            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     182        }
     183    }
     184    if (*bilevelAstrometry) {
     185        pmAstromReadBilevelMosaic (fpa, phu->header);
     186    }
     187    psFree (view);
     188    return true;
     189}

  • trunk/psastro/src/psastroDefineFiles.c

    r15669 r15891  
    2020    output->save = true;
    2121
    22     bool fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
    23     if (fixChips) {
    24         if (!psastroDefineFile (config, input->fpa, "PSASTRO.REF.ASTROM", "REF.ASTROM", PM_FPA_FILE_ASTROM, PM_DETREND_TYPE_ASTROM)) {
    25             psError (PS_ERR_IO, false, "Can't find a reference astrometry file");
     22    bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
     23    if (!status) {
     24        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     25    }
     26    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
     27    if (!status) {
     28        fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     29    }
     30    if (fixChips || useModel) {
     31        if (!psastroDefineFile (config, input->fpa, "PSASTRO.MODEL", "ASTROM.MODEL", PM_FPA_FILE_ASTROM, PM_DETREND_TYPE_ASTROM)) {
     32            psError (PS_ERR_IO, false, "Can't find an astrometry model file");
    2633            return NULL;
    2734        }
    28     }
    29 
    30     bool saveAstrom = psMetadataLookupBool (&status, recipe, "SAVE.REF.ASTROM");
    31     if (saveAstrom) {
    32         pmFPAfile *outAstrom = pmFPAfileDefineOutputFromFile (config, input, "PSASTRO.OUT.ASTROM");
    33         if (!outAstrom) {
    34             psError (PS_ERR_IO, false, "Can't find a reference astrometry file");
    35             return NULL;
    36         }
    37         outAstrom->save = true;
    3835    }
    3936

  • trunk/psastro/src/psastroEnforceChips.c

    r15670 r15891  
    11# include "psastroInternal.h"
     2# define NONLIN_TOL 0.001 /* tolerance in pixels */
    23
    3 bool psastroEnforceChips (pmConfig *config, pmFPA *fpa, psMetadata *recipe) {
     4XXX add a function to supply for all entries, as opposed to fixing the poor astrometry ones.
     5bool psastroEnforceChips (pmConfig *config, psMetadata *recipe) {
    46
    57  bool status;
    68
     9  bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
     10  if (!status) {
     11      fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
     12  }
     13  if (!fixChips) return true;
     14
    715  // identify reference astrometry table
    816  // if not defined, correction was not requested; skip step
    9   pmFPAfile *refAstrom = psMetadataLookupPtr (NULL, config->files, "REF.ASTROM");
    10   if (!refAstrom) return true;
     17  pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
     18  if (!astrom) return true;
     19
     20  // select the input data sources
     21  pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
     22  if (!input) {
     23      psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     24      return false;
     25  }
     26
     27  // make sure the astrometry model is loaded
     28  // de-activate all files except PSASTRO.MODEL
     29  // XXX do I need to supply the input->fpa->concepts to the astrom->fpa->concepts?
     30  psFree (astrom->fpa->concepts);
     31  astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts);
     32  pmFPAfileActivate (config->files, false, NULL);
     33  pmFPAfileActivate (config->files, true, "PSASTRO.MODEL");
    1134
    1235  pmFPAview *view = pmFPAviewAlloc (0);
     36
     37  // files associated with the science image
     38  if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     39      psError (PS_ERR_IO, false, "Can't load the astrometry model file");
     40      return false;
     41  }
     42
     43  // loop over all chips, replace input astrometry elements with those from astrom
     44  pmChip *obsChip = NULL;
     45  while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     46    psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
     47    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
     48
     49    // set the chip astrometry using the astrom file
     50    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
     51
     52    psFree (obsChip->toFPA);
     53    psFree (obsChip->fromFPA);
     54
     55    // XXX can we do a straight copy?  I think so, but we may need to apply a scaling factor based
     56    // on the toSky or toTPA values
     57    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
     58    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
     59
     60    // XXX test to write out adjusted header
     61    pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
     62  }
     63
     64  psFree (input->fpa->toSky);
     65  psFree (input->fpa->toTPA);
     66  psFree (input->fpa->fromTPA);
     67  input->fpa->toSky   = psMemIncrRefCounter (astrom->fpa->toSky);
     68  input->fpa->toTPA   = psMemIncrRefCounter (astrom->fpa->toTPA);
     69  input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
     70
     71  psMetadata *updates = psMetadataAlloc();
     72  pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL);
     73  psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     74  psFree (updates);
     75
     76  psFree (view);
     77  return true;
     78}
     79
     80# if (0)
    1381
    1482  // physical pixel scale in microns per pixel
    … …  
    33101    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }
    34102
    35     // set the chip astrometry using the refAstrom file
    36     pmChip *refChip = pmFPAviewThisChip (view, refAstrom->fpa);
     103    // set the chip astrometry using the astrom file
     104    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);
    37105
    38106    // bad Astrometry test:  ref pixel or angle outside nominal
    … …  
    80148// XXX for this function to work, I need to:
    81149// 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?)
    82 // 2) load the refAstrom fpa in the pmFPAfileIO stage
     150// 2) load the astrom fpa in the pmFPAfileIO stage
    83151// 3) allow for this operation to be optional
     152# endif

  • trunk/psastro/src/psastroModel.c

    r15880 r15891  
    3131    // run the full astrometry analysis (chip and/or mosaic)
    3232    if (!psastroModelAnalysis (config)) {
    33         psErrorStackPrint(stderr, "failure in psastro analysis\n");
     33        psErrorStackPrint(stderr, "failure in psastro model analysis\n");
    3434        exit (1);
    3535    }
    3636   
     37    // run the full astrometry analysis (chip and/or mosaic)
     38    if (!psastroModelAdjust (config)) {
     39        psErrorStackPrint(stderr, "failure in psastro model adjust\n");
     40        exit (1);
     41    }
     42   
     43    // save the model
     44    if (!psastroModelDataSave (config)) {
     45        psErrorStackPrint(stderr, "error saving output data\n");
     46        exit (1);
     47    }
     48
    3749    psLogMsg ("psastro", 3, "complete psastro run: %f sec\n", psTimerMark ("complete"));
    3850

  • trunk/psastro/src/psastroModelAnalysis.c

    r15887 r15891  
    1212    }
    1313
     14    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
     15    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
     16
    1417    // physical pixel scale in microns per pixel
    1518    char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP");
    … …  
    1922    }
    2023
    21     /* model analysis steps:
    22 
    23     1) determine POS_ZERO:
    24     * choose a reference chip
    25     * what is posangle of chip
    26     * posangle = FPA.POSANGLE - POS_ZERO
    27     * POS_ZERO = FPA.POSANGLE - posangle
    28 
    29     ** average over input images
    30 
    31     2) determine boresite model
    32 
    33     * find CRPIX1,2 for a series of images
    34     * fit ellipse
    35 
    36     */
    37 
    38     /***** find the POSANGLE offset *****/
     24    /* model analysis:
     25     *
     26     * determine POS_ZERO via comparison of measured and reported posangles
     27     * POS_ZERO = FPA.POSANGLE - posangle
     28     *
     29     * determine boresite model:
     30     * X = Xo + R_X cos(FPA.POSANGLE - T_0) cos(P_0) + R_Y sin(FPA.POSANGLE - T_0) sin(P_0)
     31     * Y = Yo + R_Y sin(FPA.POSANGLE - T_0) cos(P_0) - R_X cos(FPA.POSANGLE - T_0) sin(P_0)
     32     * position of reported boresite in reference chip pixels
     33     * Xo, Yo : true coordinate of boresite (rotator center) in reference chip pixels
     34     * R_X, R_Y : amplitude of boresite offset
     35     * T_0 : reference angle for rotator
     36     * P_0 : orientation of boresite ellipse
     37     */
    3938
    4039    // select the input pmFPAfile pointers
    … …  
    4443    psArray *files = psListToArray (item->data.list);
    4544
    46     // temp data vectors
     45    // data storage vectors for measurements
    4746    psVector *posZero  = psVectorAlloc (files->n, PS_TYPE_F32);
    4847    psVector *Po       = psVectorAlloc (files->n, PS_TYPE_F32);
    49     psVector *Lo       = psVectorAlloc (files->n, PS_TYPE_F32);
    50     psVector *Mo       = psVectorAlloc (files->n, PS_TYPE_F32);
     48    psVector *Xo       = psVectorAlloc (files->n, PS_TYPE_F32);
     49    psVector *Yo       = psVectorAlloc (files->n, PS_TYPE_F32);
    5150
     51    // counter for accepted measured values
    5252    int n = 0;
    5353
    54     // convert the headers for the input file into fpa astrometry terms
     54    // Re-select the output fpa.  The output fpa needs to have valid astrometry for the
     55    // refChip.  output->fpa is a copy of the pointer to one of the input->fpa, but the choice
     56    // is arbitrary.  select a new one that has an existing ref chip
     57    psFree (output->fpa);
     58    output->fpa = NULL;
     59
     60    // extract the relevant measured and reported values from the reference chip
    5561    for (int i = 0; i < files->n; i++) {
    5662        psMetadataItem *file = files->data[i];
    5763        pmFPAfile *input = file->data.V;
    5864
     65        // reported rotator position angle (this should perhaps be ROT, not POS)
    5966        double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE");
    6067        if (!status) psAbort ("missing FPA.POSANGLE");
    6168
    62         // get chip from name
     69        // get reference chip from name
    6370        pmChip *chip = pmConceptsChipFromName (input->fpa, refChip);
    6471        if (!chip) psAbort ("invalid chip name for reference");
    6572
     73        fprintf (stderr, "input %d : %zx : %zx : %zx\n", i, (size_t) input, (size_t) chip, (size_t) chip->toFPA);
    6674        if (!chip->toFPA) continue;
     75
     76        if (!output->fpa) {
     77            // this one matches
     78            output->fpa = psMemIncrRefCounter(input->fpa);
     79        }
    6780
    6881        // we have two measurements of the posangle (may be parity flipped to different quadrants)
    … …  
    7285        while (posZero->data.F32[n] <   0.0) posZero->data.F32[n] += 360.0;
    7386
    74         Po->data.F32[n] = POSANGLE * PM_RAD_DEG;
    75         Lo->data.F32[n] = chip->toFPA->x->coeff[0][0];
    76         Mo->data.F32[n] = chip->toFPA->y->coeff[0][0];
    77         fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Lo->data.F32[n], Mo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
     87        Po->data.F32[n] = POSANGLE * PM_RAD_DEG; // reported position angle
     88        Xo->data.F32[n] = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
     89        Yo->data.F32[n] = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
     90        fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
    7891        n ++;
    7992    }
    8093     
    81     Lo->n = n;
    82     Mo->n = n;
     94    Xo->n = n;
     95    Yo->n = n;
    8396    Po->n = n;
    8497    posZero->n = n;
    … …  
    88101
    89102    fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     103    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POS_ZERO", PS_META_REPLACE, "offset between obs and meas posangle", stats->sampleMedian);
    90104
    91     psastroModelFitBoresite (Lo, Mo, Po);
     105    psVector *params = psastroModelFitBoresite (Xo, Yo, Po);
     106    if (params->n != 6) psAbort ("error");
     107
     108    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_X0]);
     109    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_Y0]);
     110    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RX]);
     111    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RY]);
     112    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_T0]);
     113    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_P0]);
     114    psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip);
    92115
    93116    return true;

  • trunk/psastro/src/psastroModelBoresite.c

    r15887 r15891  
    11# include "psastroStandAlone.h"
    22
    3 // chisq = sum ((L_obs - L_fit(t))^2 + (M_obs - M_fit(t))^2)
    4 // coord[0] = measured L or measured M
     3// the full chisq is built of two associated sums over coordinates:
     4// chisq = sum ((X_obs - X_fit(t))^2 + (Y_obs - Y_fit(t))^2)
     5// we use split this into a 2x long vector and use coord[1] to distinguish the X and Y terms:
     6// coord[0] = measured X or measured Y
    57// coord[1] =          0 or          1
    68psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord) {
    … …  
    1517    float sntht = sin(dtheta);
    1618
    17     // value is L
     19    // value is X
    1820    if (coord->data.F32[1] == 0) {
    1921
    20         float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi + PAR[PAR_RM]*sntht*snphi;
     22        float value = PAR[PAR_X0] + PAR[PAR_RX]*cstht*csphi + PAR[PAR_RY]*sntht*snphi;
    2123
    2224        if (deriv) {
    2325            psF32 *dPAR = deriv->data.F32;
    24             dPAR[PAR_L0] = 1.0;
    25             dPAR[PAR_M0] = 0.0;
    26             dPAR[PAR_RL] = +cstht*csphi;
    27             dPAR[PAR_RM] = +sntht*snphi;
    28             dPAR[PAR_P0] = -PAR[PAR_RL]*cstht*snphi + PAR[PAR_RM]*sntht*csphi;
    29             dPAR[PAR_T0] =  PAR[PAR_RL]*sntht*csphi - PAR[PAR_RM]*cstht*snphi;
     26            dPAR[PAR_X0] = 1.0;
     27            dPAR[PAR_Y0] = 0.0;
     28            dPAR[PAR_RX] = +cstht*csphi;
     29            dPAR[PAR_RY] = +sntht*snphi;
     30            dPAR[PAR_P0] = -PAR[PAR_RX]*cstht*snphi + PAR[PAR_RY]*sntht*csphi;
     31            dPAR[PAR_T0] =  PAR[PAR_RX]*sntht*csphi - PAR[PAR_RY]*cstht*snphi;
    3032        }
    3133        return (value);
    3234    } 
    3335
    34     // value is M
     36    // value is Y
    3537    if (coord->data.F32[1] == 1) {
    36         float value = PAR[PAR_M0] + PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi;
     38        float value = PAR[PAR_Y0] + PAR[PAR_RY]*sntht*csphi - PAR[PAR_RX]*cstht*snphi;
    3739
    3840        if (deriv) {
    3941            psF32 *dPAR = deriv->data.F32;
    40             dPAR[PAR_L0]  = 0.0;
    41             dPAR[PAR_M0]  = 1.0;
    42             dPAR[PAR_RL]  = -cstht*snphi;
    43             dPAR[PAR_RM]  = +sntht*csphi;
    44             dPAR[PAR_P0]  = -PAR[PAR_RM]*sntht*snphi - PAR[PAR_RL]*cstht*csphi;
    45             dPAR[PAR_T0]  = -PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi;
     42            dPAR[PAR_X0]  = 0.0;
     43            dPAR[PAR_Y0]  = 1.0;
     44            dPAR[PAR_RX]  = -cstht*snphi;
     45            dPAR[PAR_RY]  = +sntht*csphi;
     46            dPAR[PAR_P0]  = -PAR[PAR_RY]*sntht*snphi - PAR[PAR_RX]*cstht*csphi;
     47            dPAR[PAR_T0]  = -PAR[PAR_RY]*sntht*csphi - PAR[PAR_RX]*cstht*snphi;
    4648        }
    4749        return (value);
    … …  
    4951    psAbort ("programming error: invalid coordinate");
    5052}
    51 
    52 # undef PAR_L0
    53 # undef PAR_RL
    54 # undef PAR_PHI
    55 # define PAR_L0  0
    56 # define PAR_RL  1
    57 # define PAR_PHI 2
    58 
    59 psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord) {
    60 
    61     psF32 *PAR = params->data.F32;
    62 
    63     float csphi = cos(PAR[PAR_PHI]);
    64     float snphi = sin(PAR[PAR_PHI]);
    65 
    66     float theta = coord->data.F32[0];
    67     float cstht = cos(theta);
    68     float sntht = sin(theta);
    69 
    70     float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi - PAR[PAR_RL]*sntht*snphi;
    71 
    72     if (deriv) {
    73         psF32 *dPAR = deriv->data.F32;
    74         dPAR[PAR_L0]  = 1.0;
    75         dPAR[PAR_RL]  = cstht*csphi - sntht*snphi;
    76         dPAR[PAR_PHI] = -PAR[PAR_RL]*cstht*snphi - PAR[PAR_RL]*sntht*csphi;
    77     }
    78     return (value);
    79 }
    80 
    81 psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord) {
    82 
    83     psF32 *PAR = params->data.F32;
    84 
    85     float csphi = cos(PAR[PAR_PHI]);
    86     float snphi = sin(PAR[PAR_PHI]);
    87 
    88     float theta = coord->data.F32[0];
    89     float cstht = cos(theta);
    90     float sntht = sin(theta);
    91 
    92     float value = PAR[PAR_L0] + PAR[PAR_RL]*sntht*csphi + PAR[PAR_RL]*cstht*snphi;
    93 
    94     if (deriv) {
    95         psF32 *dPAR = deriv->data.F32;
    96         dPAR[PAR_L0]  = 1.0;
    97         dPAR[PAR_RL]  = +sntht*csphi + cstht*snphi;
    98         dPAR[PAR_PHI] = -PAR[PAR_RL]*sntht*snphi + PAR[PAR_RL]*cstht*csphi;
    99     }
    100     return (value);
    101 }

  • trunk/psastro/src/psastroModelDataLoad.c

    r15880 r15891  
    3333
    3434    // select the input data sources
    35     pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
    36     if (!output) psAbort ("PSASTRO.MODEL not listed in config->files");
     35    pmFPAfile *output = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUT.MODEL");
     36    if (!output) psAbort ("PSASTRO.OUT.MODEL not listed in config->files");
    3737
    3838    pmFPAview *view = pmFPAviewAlloc (0);

  • trunk/psastro/src/psastroModelFitBoresite.c

    r15887 r15891  
    22
    33// we now have a set of observed L,M values.  fit these to the boresite model
    4 bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
     4psVector *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po) {
    55
    6     assert (Lo->n > 2);
    7     assert (Lo->n == Mo->n);
    8     assert (Lo->n == Po->n);
     6    assert (Xo->n > 2);
     7    assert (Xo->n == Yo->n);
     8    assert (Xo->n == Po->n);
    99
    1010    // arrays to hold the data to be fitted
    11     psArray *x = psArrayAlloc(2*Lo->n);
    12     psVector *y = psVectorAlloc(2*Lo->n, PS_TYPE_F32);
     11    psArray *x = psArrayAlloc(2*Xo->n);
     12    psVector *y = psVectorAlloc(2*Xo->n, PS_TYPE_F32);
    1313
    1414    int n = 0;
    15     for (int i = 0; i < Lo->n; i++) {
     15    for (int i = 0; i < Xo->n; i++) {
    1616
    1717        psVector *coord = NULL;
    1818
    19         // L coordinate value
     19        // X coordinate value
    2020        coord = psVectorAlloc (2, PS_TYPE_F32);
    2121        coord->data.F32[1] = 0.0;
    2222        coord->data.F32[0] = Po->data.F32[i];
    2323        x->data[n] = coord;
    24         y->data.F32[n] = Lo->data.F32[i];
     24        y->data.F32[n] = Xo->data.F32[i];
    2525        n++;
    2626       
    27         // M coordinate value
     27        // Y coordinate value
    2828        coord = psVectorAlloc (2, PS_TYPE_F32);
    2929        coord->data.F32[1] = 1.0;
    3030        coord->data.F32[0] = Po->data.F32[i];
    3131        x->data[n] = coord;
    32         y->data.F32[n] = Mo->data.F32[i];
     32        y->data.F32[n] = Yo->data.F32[i];
    3333        n++;
    3434    }   
    … …  
    4747    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    4848
    49     // center (Lo) = mean(Lo), RL = range / 2
    50     psVectorStats (stats, Lo, NULL, NULL, 0);
    51     params->data.F32[PAR_L0] = stats->sampleMean;
    52     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
     49    // center (Xo) = mean(Xo), RX = range / 2
     50    psVectorStats (stats, Xo, NULL, NULL, 0);
     51    params->data.F32[PAR_X0] = stats->sampleMean;
     52    params->data.F32[PAR_RX] = (stats->max - stats->min) / 2.0;
    5353
    54     // center (Mo) = mean(Mo), RM = range / 2
    55     psVectorStats (stats, Mo, NULL, NULL, 0);
    56     params->data.F32[PAR_M0] = stats->sampleMean;
    57     params->data.F32[PAR_RM] = (stats->max - stats->min) / 2.0;
     54    // center (Yo) = mean(Yo), RY = range / 2
     55    psVectorStats (stats, Yo, NULL, NULL, 0);
     56    params->data.F32[PAR_Y0] = stats->sampleMean;
     57    params->data.F32[PAR_RY] = (stats->max - stats->min) / 2.0;
    5858
    5959    params->data.F32[PAR_P0] = 0.0;
    … …  
    6363    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    6464   
    65     fprintf (stderr, "guess values:\n");
    66     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    67     fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
    68     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    69     fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
    70     fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
    71     fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
     65    // fprintf (stderr, "guess values:\n");
     66    // fprintf (stderr, "Xo:  %f\n", params->data.F32[PAR_X0]);
     67    // fprintf (stderr, "Yo:  %f\n", params->data.F32[PAR_Y0]);
     68    // fprintf (stderr, "RX:  %f\n", params->data.F32[PAR_RX]);
     69    // fprintf (stderr, "RY:  %f\n", params->data.F32[PAR_RY]);
     70    // fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
     71    // fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
    7272
    7373    // XXX skip the weights for now
    … …  
    7575
    7676    fprintf (stderr, "fitted values:\n");
    77     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    78     fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
    79     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    80     fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
     77    fprintf (stderr, "Xo:  %f\n", params->data.F32[PAR_X0]);
     78    fprintf (stderr, "Yo:  %f\n", params->data.F32[PAR_Y0]);
     79    fprintf (stderr, "RX:  %f\n", params->data.F32[PAR_RX]);
     80    fprintf (stderr, "RY:  %f\n", params->data.F32[PAR_RY]);
    8181    fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
    8282    fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
    8383
    84     return true;
     84    return params;
    8585}
    86 
    87 # if (0)
    88 # undef PAR_L0
    89 # undef PAR_RL
    90 # undef PAR_PHI
    91 # define PAR_L0  0
    92 # define PAR_RL  1
    93 # define PAR_PHI 2
    94 
    95 // we now have a set of observed L,M values.  fit these to the boresite model
    96 bool psastroModelFitBoresite_L (psVector *Lo, psVector *Mo, psVector *Po) {
    97 
    98     assert (Lo->n > 2);
    99     assert (Lo->n == Mo->n);
    100     assert (Lo->n == Po->n);
    101 
    102     // arrays to hold the data to be fitted
    103     psArray *x = psArrayAlloc(Lo->n);
    104     psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
    105 
    106     for (int i = 0; i < Lo->n; i++) {
    107 
    108         psVector *coord = NULL;
    109 
    110         // L coordinate value
    111         coord = psVectorAlloc (1, PS_TYPE_F32);
    112         coord->data.F32[0] = Po->data.F32[i];
    113         x->data[i] = coord;
    114         y->data.F32[i] = Lo->data.F32[i];
    115     }   
    116 
    117     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
    118    
    119     // create the minimization constraints
    120     psMinConstraint *constraint = psMinConstraintAlloc();
    121 
    122     // XXX for now, no parameter masks, skip checkLimits
    123     // constraint->checkLimits = psastroModelBoresiteLimits;
    124 
    125     // make an initial guess:
    126     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    127 
    128     // center (Lo) = mean(Lo), RL = range / 2
    129     psVectorStats (stats, Lo, NULL, NULL, 0);
    130     params->data.F32[PAR_L0] = stats->sampleMean;
    131     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
    132     params->data.F32[PAR_PHI] = 0.0;
    133 
    134     psMinimization *myMin = psMinimizationAlloc (15, 0.001);
    135     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    136    
    137     fprintf (stderr, "guess values:\n");
    138     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    139     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    140     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    141 
    142     // XXX skip the weights for now
    143     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteL);
    144 
    145     fprintf (stderr, "fitted values:\n");
    146     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    147     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    148     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    149 
    150     return true;
    151 }
    152 
    153 // we now have a set of observed L,M values.  fit these to the boresite model
    154 bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
    155 
    156     assert (Lo->n > 2);
    157     assert (Lo->n == Mo->n);
    158     assert (Lo->n == Po->n);
    159 
    160     // arrays to hold the data to be fitted
    161     psArray *x = psArrayAlloc(Lo->n);
    162     psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
    163 
    164     for (int i = 0; i < Lo->n; i++) {
    165 
    166         psVector *coord = NULL;
    167 
    168         // L coordinate value
    169         coord = psVectorAlloc (1, PS_TYPE_F32);
    170         coord->data.F32[0] = Po->data.F32[i];
    171         x->data[i] = coord;
    172         y->data.F32[i] = Mo->data.F32[i];
    173     }   
    174 
    175     psVector *params = psVectorAlloc (3, PS_TYPE_F32);
    176    
    177     // create the minimization constraints
    178     psMinConstraint *constraint = psMinConstraintAlloc();
    179 
    180     // XXX for now, no parameter masks, skip checkLimits
    181     // constraint->checkLimits = psastroModelBoresiteLimits;
    182 
    183     // make an initial guess:
    184     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
    185 
    186     // center (Lo) = mean(Lo), RL = range / 2
    187     psVectorStats (stats, Mo, NULL, NULL, 0);
    188     params->data.F32[PAR_L0] = stats->sampleMean;
    189     params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
    190     params->data.F32[PAR_PHI] = 0.0;
    191 
    192     psMinimization *myMin = psMinimizationAlloc (15, 0.001);
    193     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    194    
    195     fprintf (stderr, "guess values:\n");
    196     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    197     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    198     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    199 
    200     // XXX skip the weights for now
    201     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteM);
    202 
    203     fprintf (stderr, "fitted values:\n");
    204     fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
    205     fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
    206     fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
    207 
    208     return true;
    209 }
    210 # endif

  • trunk/psastro/src/psastroModelParseCamera.c

    r15880 r15891  
    1818       
    1919    // set up an output fpa structure & file based on the last input file
    20     pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.MODEL");
     20    pmFPAfile *output = pmFPAfileDefineOutput (config, input->fpa, "PSASTRO.OUT.MODEL");
    2121    if (!output) {
    22         psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA for PSASTRO.MODEL from input");
     22        psError(PSASTRO_ERR_CONFIG, false, "Failed to build FPA for PSASTRO.OUT.MODEL from input");
    2323        return false;
    2424    }

  • trunk/psastro/src/psastroMosaicAstrom.c

    r15638 r15891  
    2323    // before we do object matches, we need to (optionally) fix failed chips.  We compare chips with
    2424    // the supplied mosaic model.  Adjust significant outliers to match model. 
    25     if (!psastroEnforceChips (config, fpa, recipe)) {
     25    if (!psastroFixChips (config, recipe)) {
    2626        psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
    2727        return false;

  • trunk/psastro/src/psastroStandAlone.h

    r15887 r15891  
    1717bool              psastroDataLoad (pmConfig *config);
    1818
    19 pmConfig *psastroModelArguments (int argc, char **argv);
    20 bool psastroModelParseCamera (pmConfig *config);
    21 bool psastroModelDataLoad (pmConfig *config);
    22 bool psastroModelAnalysis (pmConfig *config);
     19pmConfig         *psastroModelArguments (int argc, char **argv);
     20bool              psastroModelParseCamera (pmConfig *config);
     21bool              psastroModelDataLoad (pmConfig *config);
     22bool              psastroModelAnalysis (pmConfig *config);
     23bool              psastroModelAdjust (pmConfig *config);
     24bool              psastroModelDataSave (pmConfig *config);
    2325
    24 bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po);
    25 psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord);
    26 
    27 psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord);
    28 psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord);
     26psVector         *psastroModelFitBoresite (psVector *Xo, psVector *Yo, psVector *Po);
     27psF32             psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord);
    2928
    3029// these are used to define the boresite model parameters
    31 # define PAR_L0  0  // Lo = params[0]
    32 # define PAR_M0  1  // Mo = params[1]
    33 # define PAR_RL  2  // RL = params[2]
    34 # define PAR_RM  3  // RM = params[3]
     30# define PAR_X0  0  // Xo = params[0]
     31# define PAR_Y0  1  // Yo = params[1]
     32# define PAR_RX  2  // RX = params[2]
     33# define PAR_RY  3  // RY = params[3]
    3534# define PAR_P0  4  // P0 = params[4]
    3635# define PAR_T0  5  // phi = params[4]
Note: See TracChangeset for help on using the changeset viewer.