IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36819


Ignore:
Timestamp:
Jun 5, 2014, 4:06:26 PM (12 years ago)
Author:
eugene
Message:

add Koppenhoefer correction to psastro

Location:
branches/eam_branches/ipp-20140423/psastro/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140423/psastro/src/psastro.h

    r36737 r36819  
    6666bool              psastroAnalysis (pmConfig *config, psMetadata *stats);
    6767
    68 bool              psastroConvertFPA (pmFPA *fpa, psMetadata *recipe);
    69 bool              psastroConvertReadout (pmReadout *readout, psMetadata *recipe);
     68bool              psastroConvertFPA (pmConfig *config, pmFPA *fpa, psMetadata *recipe);
     69bool              psastroConvertReadout (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe);
     70bool              psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars);
     71
    7072psArray          *pmSourceToAstromObj (psArray *sources);
    7173bool              psastroAstromGuess (int *nStars, pmConfig *config);

  • branches/eam_branches/ipp-20140423/psastro/src/psastroArguments.c

    r35715 r36819  
    6060        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.USE.MODEL", PS_META_REPLACE, "", true);
    6161    }
    62     // define the reference astrometry file
     62    // define the reference astrometry file (add a container ASTROM.MODEL to config->arguments if -astrommodel (file) is found)
    6363    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist");
     64
     65    // define the koppenhoefer correction file (add a container KH.CORRECT to config->arguments if -kh-correct (file) is found)
     66    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "KH.CORRECT", "-kh-correct", NULL);
    6467
    6568    // define the reference astrometry file

  • branches/eam_branches/ipp-20140423/psastro/src/psastroConvert.c

    r31661 r36819  
    1616static psArray *chooseStars(psArray *inStars, char *listName, psArray *sources, psVector *index, int nMax, float iMagMin, float iMagMax, pmSourceMode skip);
    1717
    18 bool psastroConvertFPA (pmFPA *fpa, psMetadata *recipe) {
     18bool psastroConvertFPA (pmConfig *config, pmFPA *fpa, psMetadata *recipe) {
    1919
    2020    pmChip *chip;
    … …  
    3535                if (! readout->data_exists) { continue; }
    3636
    37                 psastroConvertReadout (readout, recipe);
     37                if (!psastroConvertReadout (config, view, readout, recipe)) {
     38                  psError(PSASTRO_ERR_CONFIG, true, "problem converting readout\n");
     39                  return false;
     40                }
    3841            }
    3942        }
    … …  
    4447
    4548// pass/apply the WCS information?
    46 bool psastroConvertReadout (pmReadout *readout, psMetadata *recipe) {
     49bool psastroConvertReadout (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe) {
    4750
    4851    bool status;
    … …  
    5760    // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?)
    5861    psArray *inStars = pmSourceToAstromObj (sources);
     62
     63    // apply Koppenhoefer correction if needed
     64    if (!psastroCorrectKH (config, view, readout, recipe, inStars)){
     65      psError(PSASTRO_ERR_CONFIG, true, "failed to correct Koppenhoefer Effect\n");
     66      return false;
     67    }
    5968
    6069    // sort in ascending magnitude order
    … …  
    183192        psF32 *dPAR = model->dparams->data.F32;
    184193
     194        // Note from EAM: I measured the KH correction from the nightly science database
     195        // The KH effect was active until 2011/05/11 when the camera voltages were modified
     196        // Between 2011/09/06 and 2013/04/30, there was a bug in pmPSF_ModelToAxes with the
     197        // result that the reported values were too large by a factor of sqrt(2).  Fortunately,
     198        // this did not affect the data in the nightly science DVO used to measure the
     199        // effect (all affected data was processed BEFORE the bug was introduced), and it
     200        // does not affect any of the PV2 or PV3 data processed AFTER the bug was fixed.
     201        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
     202
    185203        pmAstromObj *obj = pmAstromObjAlloc ();
    186204
    … …  
    192210        obj->Mag = source->psfMag;
    193211        obj->dMag = source->psfMagErr;
     212        obj->SBinst = source->psfMag + 5.0*log10(axes.major);
    194213
    195214        // XXX do we have the information giving the readout and cell offset?
    … …  
    221240}
    222241
     242static float pltScale[] = {
     243+0.0,
     244-0.255347177386,
     245-0.255585625172,
     246-0.255733260512,
     247-0.255747352242,
     248-0.255611001253,
     249-0.255311933756,
     250 +0.0,
     251 +0.0,
     252 +0.0,
     253-0.255296201766,
     254-0.255688636720,
     255-0.256161080003,
     256-0.256464074850,
     257-0.256324860930,
     258-0.256147363186,
     259-0.255734743476,
     260-0.255422729909,
     261 +0.0,
     262 +0.0,
     263-0.255508323163,
     264-0.256142899513,
     265-0.255493895173,
     266-0.256997376740,
     267-0.256935591102,
     268-0.256644588679,
     269-0.256123321533,
     270-0.255515242517,
     271 +0.0,
     272 +0.0,
     273-0.255633033872,
     274-0.256356916189,
     275-0.256929253668,
     276-0.257250442505,
     277-0.257260403216,
     278-0.256962024927,
     279-0.256387320280,
     280-0.255681164920,
     281 +0.0,
     282 +0.0,
     283-0.255671992302,
     284-0.256371312678,
     285-0.256936249495,
     286-0.257237892330,
     287-0.257248004675,
     288-0.256911942035,
     289-0.256365648448,
     290-0.255632205397,
     291 +0.0,
     292 +0.0,
     293-0.255455072403,
     294-0.256080205411,
     295-0.256570428520,
     296-0.256908982217,
     297-0.256899240971,
     298-0.256607972383,
     299-0.256107793808,
     300-0.255482232273,
     301 +0.0,
     302 +0.0,
     303-0.255237901688,
     304-0.255669494152,
     305-0.256101093233,
     306-0.256372770309,
     307-0.256373448133,
     308-0.256118197262,
     309-0.255689391643,
     310-0.255279636681,
     311 +0.0,
     312 +0.0,
     313 +0.0,
     314-0.255286350489,
     315-0.255517200172,
     316-0.255689459741,
     317-0.255668706447,
     318-0.255539349556,
     319-0.255363308907,
     320};
     321
     322bool psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars) {
     323
     324  bool status;
     325
     326  psAssert (readout, "missing readout");
     327  psAssert (readout->parent, "missing cell");
     328  psAssert (readout->parent->parent, "missing chip");
     329
     330  // should we stay or should we go?
     331  bool apply = psMetadataLookupBool (&status, recipe, "KH.CORRECT.APPLY.EXP");
     332  if (!apply) return true;
     333
     334  // get the chip ID from the chip name:
     335  pmChip *chip = readout->parent->parent;
     336  char *chipName = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
     337
     338  psAssert (strlen(chipName) == 4, "error in chip name");
     339  psAssert (chipName[0] == 'X', "error in chip name");
     340  psAssert (chipName[1] == 'Y', "error in chip name");
     341 
     342  int chipID = atoi (&chipName[2]);
     343
     344  // XXX hardwired list of chips to correct
     345  // some notes:
     346  // XY24 surpisingly does not need correction
     347  // XY27 has poor astrometry with a weird correction for all mags
     348  // XY36 surpisingly does not need correction
     349  // XY67 has poor astrometry with a weird correction for all mags
     350
     351  // correct this chip? (raise an error for unexpected chipID values)
     352  switch (chipID) {
     353    case  1: // do not correct
     354    case  2: // do not correct
     355    case  3: // do not correct
     356    case 10: // do not correct
     357    case 11: // do not correct
     358    case 12: // do not correct
     359    case 13: // do not correct
     360    case 20: // do not correct
     361    case 21: // do not correct
     362    case 22: // do not correct
     363    case 23: // do not correct
     364    case 24: // do not correct
     365    case 27: // do not correct
     366    case 30: // do not correct
     367    case 31: // do not correct
     368    case 32: // do not correct
     369    case 33: // do not correct
     370    case 36: // do not correct
     371    case 44: // do not correct
     372    case 45: // do not correct
     373    case 46: // do not correct
     374    case 47: // do not correct
     375    case 54: // do not correct
     376    case 55: // do not correct
     377    case 56: // do not correct
     378    case 57: // do not correct
     379    case 64: // do not correct
     380    case 65: // do not correct
     381    case 66: // do not correct
     382    case 67: // do not correct
     383    case 74: // do not correct
     384    case 75: // do not correct
     385    case 76: // do not correct
     386      return true;
     387    case  4: // correct
     388    case  5: // correct
     389    case  6: // correct
     390    case 14: // correct
     391    case 15: // correct
     392    case 16: // correct
     393    case 17: // correct
     394    case 25: // correct
     395    case 26: // correct
     396    case 34: // correct
     397    case 35: // correct
     398    case 37: // correct
     399    case 40: // correct
     400    case 41: // correct
     401    case 42: // correct
     402    case 43: // correct
     403    case 50: // correct
     404    case 51: // correct
     405    case 52: // correct
     406    case 53: // correct
     407    case 60: // correct
     408    case 61: // correct
     409    case 62: // correct
     410    case 63: // correct
     411    case 71: // correct
     412    case 72: // correct
     413    case 73: // correct
     414      break;
     415    default:
     416      psAbort ("chipID is invalid");
     417  }
     418
     419  // grab the KH correction file
     420  pmFPAfile *KHfile = psMetadataLookupPtr (NULL, config->files, "PSASTRO.KH.CORRECT");
     421  if (!KHfile) {
     422    psError(PM_ERR_CONFIG, false, "KH correction file not found");
     423    return false;
     424  }
     425
     426  // grab the corresponding chip
     427  pmChip *KHchip = pmFPAviewThisChip (view, KHfile->fpa);
     428  psAssert (KHchip, "found KH file, but not chip?");
     429
     430  // grab the correction spline
     431  KHcorrectData *spline = psMetadataLookupPtr (&status, KHchip->analysis, "KH.CORRECT");
     432  if (!spline) {
     433    psError(PM_ERR_CONFIG, false, "KH correction not found");
     434    return false;
     435  }
     436
     437  // the hard-wired array of plate-scales per chip above come from
     438  // a single smf.  in there, a negative plate scale means parity is flipped
     439  // on the sky.
     440  float myPltScale = fabs(pltScale[chipID]);
     441
     442  // apply the correction to the detections
     443  for (int i = 0; i < inStars->n; i++) {
     444    pmAstromObj *obj = inStars->data[i];
     445
     446    float Xraw = obj->pix->x;
     447    float SBinst = obj->SBinst;
     448
     449    float dX = KHcorrectApply (spline, SBinst);
     450
     451    float Xfix = Xraw + dX / myPltScale;
     452
     453    obj->pix->x = Xfix;
     454  }
     455
     456  return true;
     457
     458}
     459

  • branches/eam_branches/ipp-20140423/psastro/src/psastroDataLoad.c

    r21409 r36819  
    4242    }
    4343
     44    // apply Koppenhoeffer correction to this exposure?
     45    bool applyKH = psMetadataLookupBool (NULL, recipe, "KH.CORRECT.APPLY.EXP");
     46    if (applyKH) {
     47        // make sure the KH correction file is loaded.  de-activate all files except PSASTRO.KH.CORRECT
     48        pmFPAfileActivate (config->files, false, NULL);
     49        pmFPAfileActivate (config->files, true, "PSASTRO.KH.CORRECT");
     50
     51        pmFPAview *viewKH = pmFPAviewAlloc (0);
     52
     53        // files associated with the science image
     54        if (!pmFPAfileIOChecks (config, viewKH, PM_FPA_BEFORE)) {
     55            psError (PS_ERR_IO, false, "Can't load the Koppenhoeffer Correction file");
     56            return false;
     57        }
     58        psFree (viewKH);
     59    }
     60
    4461    // select the input data sources
    4562    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    … …  
    7390                if (!readout->data_exists) { continue; }
    7491
    75                 if (!psastroConvertReadout (readout, recipe)) ESCAPE;
     92                if (!psastroConvertReadout (config, view, readout, recipe)) ESCAPE;
    7693
    7794                if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;

  • branches/eam_branches/ipp-20140423/psastro/src/psastroDefineFiles.c

    r36441 r36819  
    4343        if (!psastroDefineFile (config, input->fpa, "PSASTRO.MODEL", "ASTROM.MODEL", PM_FPA_FILE_ASTROM_MODEL, PM_DETREND_TYPE_ASTROM)) {
    4444            psError (PS_ERR_IO, false, "Can't find an astrometry model file");
    45             return NULL;
    46         }
    47     }
     45            return false;
     46        }
     47    }
     48
     49    bool KHapply = psMetadataLookupBool (&status, recipe, "KH.CORRECT.APPLY");
     50    if (KHapply) {
     51        // Get the time from FPA.TIME
     52        psTime *time = psMetadataLookupPtr(NULL, input->fpa->concepts, "FPA.TIME");
     53        if (time->sec == 0 && time->nsec == 0) {
     54            psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, but date/time of exposure not found");
     55            return false;
     56        }
     57        // what is the appropriate date range for the correction
     58        char *KHendDate = psMetadataLookupStr (&status, recipe, "KH.CORRECT.ENDDATE");
     59        if (!KHendDate) {
     60            psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, but KH.CORRECT.ENDDATE not in recipe");
     61            return false;
     62        }
     63        psTime *endtime = psTimeFromString (KHendDate, PS_TIME_TAI);
     64        if (!endtime) {
     65            psError (PM_ERR_CONFIG, false, "KH.CORRECT.APPLY requested, KH.CORRECT.ENDDATE has an invalid date/time");
     66            return false;
     67        }
     68        if (psTimeDelta (endtime, time) < 0) {
     69            KHapply = false;
     70            psLogMsg ("psastro", PS_LOG_INFO, "Koppenhoefer correction requested, but exposure after end date, correction skipped");
     71        } else {
     72            psLogMsg ("psastro", PS_LOG_INFO, "Koppenhoefer correction requested, exposure in valid date range, correction may be applied");
     73        }   
     74    }
     75    if (KHapply) {
     76        if (!psastroDefineFile (config, input->fpa, "PSASTRO.KH.CORRECT", "KH.CORRECT", PM_FPA_FILE_KH_CORRECT, PM_DETREND_TYPE_KH_CORRECT)) {
     77            psError (PS_ERR_IO, false, "Can't find a Koppenhoefer Correction file");
     78            return NULL;
     79        }
     80    }
     81    // the bool above in the recipe says "apply generally for this camera".  here we are asserting that
     82    // we can and should apply for this exposure XXX change the name?
     83    psMetadataAddBool (recipe, PS_LIST_TAIL, "KH.CORRECT.APPLY.EXP", PS_META_REPLACE, "", KHapply);
    4884
    4985    // this step is optional

  • branches/eam_branches/ipp-20140423/psastro/src/psastroLoadCrosstalk.c

    r36671 r36819  
    175175
    176176            float x_t_chip,y_t_chip;
    177             int faint_ct = 0;
     177            // int faint_ct = 0; NOTE: not currently used
    178178            pmCellCoordsForChip(&x_cell,&y_cell,CTsourceCell,ref->chip->x,ref->chip->y);
    179179
    … …  
    229229                Xt = 3;
    230230                Ut = 2;
    231                 faint_ct = 1;
     231                // faint_ct = 1; NOTE: not currently used
    232232              }
    233233            }
    … …  
    251251                Xt = 4;
    252252                Ut = 2;
    253                 faint_ct = 1;
     253                // faint_ct = 1; NOTE: not currently used
    254254              }
    255255            }
Note: See TracChangeset for help on using the changeset viewer.