IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2014, 9:37:38 PM (12 years ago)
Author:
eugene
Message:

dcr + kh correction

Location:
branches/eam_branches/ipp-20140904/Ohana/src/uniphot
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/Makefile

    r37576 r37744  
    112112$(SRC)/setastrom.$(ARCH).o          \
    113113$(SRC)/initialize_setastrom.$(ARCH).o \
     114$(SRC)/dcr_correction.$(ARCH).o \
     115$(SRC)/kh_correction.$(ARCH).o \
    114116$(SRC)/astrom_correction.$(ARCH).o \
    115117$(SRC)/update_dvo_setastrom.$(ARCH).o \
    … …  
    125127$(SRC)/update_dvo_setastrom.$(ARCH).o \
    126128$(SRC)/update_catalog_setastrom.$(ARCH).o \
     129$(SRC)/dcr_correction.$(ARCH).o \
     130$(SRC)/kh_correction.$(ARCH).o \
    127131$(SRC)/astrom_correction.$(ARCH).o \
    128132$(SRC)/initialize_setastrom.$(ARCH).o

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/include/setastrom.h

    r35103 r37744  
    2525char        *SINGLE_CPT;
    2626
     27char        *KH_FILE;
     28char        *DCR_FILE;
     29
     30int          KH_RESET;
     31int          DCR_RESET;
     32
    2733SkyRegion    UserPatch;
    2834
    … …  
    4652// int           SetSignals                      PROTO((void));
    4753
    48 int           update_dvo_setastrom            PROTO((char *corrfile));
    49 int           update_dvo_setastrom_parallel   PROTO((SkyTable *sky, char *corrfile));
     54int           update_dvo_setastrom            PROTO(());
     55int           update_dvo_setastrom_parallel   PROTO((SkyTable *sky));
    5056
    5157int           update_catalog_setastrom        PROTO((Catalog *catalog));
    5258
    5359// spline correction functions
    54 int load_astrom_correction (char *filename);
    55 int get_astrom_correction (int sub, int chip, double *dX, double *dY, float Minst);
    5660double spline_apply_dbl (double *x, double *y, double *y2, int N, double X);
     61
     62int load_kh_correction (char *filename);
     63int get_kh_correction (int sub, int chip, double *dX, double *dY, float Minst);
     64
     65int load_dcr_correction (char *filename);
     66int get_dcr_correction (int filter, double *dX, double *dY, float Color);

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/src/astrom_correction.c

    r35103 r37744  
    22
    33// astrometry correction functions:
    4 // load the correction set from a file
    5 
    6 // the correction is saved as a set of splines, one per file extension
    7 // extensions have names of the form XYnn.Dx.Td
    8 
    9 // where nn = 00 - 77, the chip number
    10 // x is X or Y, the correction direction (dir)
    11 // d is 0 or 1, pre or post camera correction (sub)
    12 
    13 // these are stored in an N-D array of (Spline *) splineset[sub][dir][chip]
    14 
    15 # define N_SUB 2
    16 # define N_DIR 2
    17 # define N_CHIP 78
    18 // note that we are allocating pointers for XY00 - XY77, but only the octal elements are set (and not the 00,07,70,77 ones)
    19 
    20 Spline ****splineset;
    21 
    22 int load_astrom_correction (char *filename) {
    23 
    24   int Ncol;
    25   off_t Nrow;
    26 
    27   char type[16];
    28   char extname[80];
    29 
    30   Header header;
    31   Header theader;
    32   FTable ftable;
    33 
    34   /* open file for input */
    35   FILE *f = fopen (filename, "r");
    36   if (f == (FILE *) NULL) {
    37     gprint (GP_ERR, "can't open file for read : %s\n", filename);
    38     return FALSE;
    39   }
    40 
    41   // read the PHU header
    42   if (!gfits_load_header (f, &header)) return (FALSE);
    43  
    44   int Nbytes = gfits_data_size (&header);
    45   fseeko (f, Nbytes, SEEK_CUR);
    46   // XXX gfits_free_header (&header);
    47 
    48   ftable.header = &theader;
    49 
    50   int sub, dir, chip;
    51 
    52   // we have a N-D array of splines:
    53   // splineset[sub][dir][chip]
    54 
    55   ALLOCATE (splineset, Spline ***, N_SUB);
    56   for (sub = 0; sub < N_SUB; sub++) {
    57     ALLOCATE (splineset[sub], Spline **, N_DIR);
    58     for (dir = 0; dir < N_DIR; dir++) {
    59       ALLOCATE (splineset[sub][dir], Spline *, N_CHIP);
    60       for (chip = 0; chip < N_CHIP; chip ++) {
    61         splineset[sub][dir][chip] = NULL;
    62       }
    63     }
    64   }
    65 
    66   while (TRUE) {
    67     // load data for this header : if not found, assume we hit the end of the file
    68     if (!gfits_load_header (f, &theader)) break;
    69    
    70     if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (FALSE);
    71    
    72     if (!gfits_scan (&theader, "EXTNAME", "%s", 1, extname)) return (FALSE);
    73 
    74     // figure out the chip ID, diretion, and sequence from the EXTNAME
    75 
    76     // EXTNAME is of the form XY%d%d.D%s.T%s where
    77     char chipname[3];
    78     chipname[0] = extname[2];
    79     chipname[1] = extname[3];
    80     chipname[2] = 0;
    81     int chip = atoi(chipname);
    82     myAssert (chip >   0, "error in extname %s not(chip >   0)\n", extname);
    83     myAssert (chip <  77, "error in extname %s not(chip <  77)\n", extname);
    84     myAssert (chip !=  7, "error in extname %s not(chip !=  7)\n", extname);
    85     myAssert (chip != 70, "error in extname %s not(chip != 70)\n", extname);
    86 
    87     int dir = -1;
    88     if (extname[6] == 'X') {
    89       dir = 0;
    90     }
    91     if (extname[6] == 'Y') {
    92       dir = 1;
    93     }
    94     myAssert (dir != -1, "error in extname %s (dir)\n", extname);
    95 
    96     int sub = -1;
    97     if (extname[9] == '0') {
    98       sub = 0;
    99     }
    100     if (extname[9] == '1') {
    101       sub = 1;
    102     }
    103     myAssert (sub != -1, "error in extname %s (sub)\n", extname);
    104 
    105     // need to create and assign to flat-field correction
    106     double *xk = gfits_get_bintable_column_data (&theader, &ftable, "X_KNOT", type, &Nrow, &Ncol);
    107     myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    108    
    109     // need to create and assign to flat-field correction
    110     double *yk = gfits_get_bintable_column_data (&theader, &ftable, "Y_KNOT", type, &Nrow, &Ncol);
    111     myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    112    
    113     // need to create and assign to flat-field correction
    114     double *y2 = gfits_get_bintable_column_data (&theader, &ftable, "DY2_DX", type, &Nrow, &Ncol);
    115     myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    116    
    117     ALLOCATE (splineset[sub][dir][chip], Spline, 1);
    118     splineset[sub][dir][chip]->xk = xk;
    119     splineset[sub][dir][chip]->yk = yk;
    120     splineset[sub][dir][chip]->y2 = y2;
    121     splineset[sub][dir][chip]->Nknots = Nrow;
    122   }
    123   return (TRUE);
    124 }
    125 
    126 int get_astrom_correction (int sub, int chip, double *dX, double *dY, float Minst) {
    127 
    128   Spline *spline = NULL;
    129 
    130   *dX = 0.0;
    131   *dY = 0.0;
    132 
    133   // XXX need to include some limits on Minst?
    134 
    135   spline = splineset[sub][0][chip];
    136   if (spline) {
    137     *dX = spline_apply_dbl (spline->xk, spline->yk, spline->y2, spline->Nknots, Minst);
    138   }
    139 
    140   spline = splineset[sub][1][chip];
    141   if (spline) {
    142     *dY = spline_apply_dbl (spline->xk, spline->yk, spline->y2, spline->Nknots, Minst);
    143   }
    144   return (TRUE);
    145 }
     4// KH correction
     5// DCR correction
    1466
    1477/* evaluate spline for x, y, y2 at X */

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/src/initialize_setastrom.c

    r35103 r37744  
    22
    33void usage_setastrom () {
    4     fprintf (stderr, "USAGE: setastrom (corrections) [options]\n");
     4    fprintf (stderr, "USAGE: setastrom [options]\n");
    55    fprintf (stderr, "  options:\n");
     6    fprintf (stderr, "    -DCR (file) : supply DCR correction splines\n");
     7    fprintf (stderr, "    -KH (file) : supply Koppenhoefer correction splines\n");
    68    fprintf (stderr, "    -v : verbose mode\n");
    79    fprintf (stderr, "    -region Rmin Rmax Dmin Dmax\n");
    … …  
    4951  int N;
    5052
     53  KH_FILE = NULL;
     54  if ((N = get_argument (argc, argv, "-KH"))) {
     55    remove_argument (N, &argc, argv);
     56    char *tmpfile = strcreate (argv[N]);
     57    KH_FILE = abspath (tmpfile, DVO_MAX_PATH);
     58    remove_argument (N, &argc, argv);
     59  }
     60  KH_RESET = FALSE;
     61  if ((N = get_argument (argc, argv, "-KH-reset"))) {
     62    remove_argument (N, &argc, argv);
     63    KH_RESET = TRUE;
     64  }
     65
     66  DCR_FILE = NULL;
     67  if ((N = get_argument (argc, argv, "-DCR"))) {
     68    remove_argument (N, &argc, argv);
     69    char *tmpfile = strcreate (argv[N]);
     70    DCR_FILE = abspath (tmpfile, DVO_MAX_PATH);
     71    remove_argument (N, &argc, argv);
     72  }
     73  DCR_RESET = FALSE;
     74  if ((N = get_argument (argc, argv, "-DCR-reset"))) {
     75    remove_argument (N, &argc, argv);
     76    DCR_RESET = TRUE;
     77  }
     78
    5179  SINGLE_CPT = NULL;
    5280  if ((N = get_argument (argc, argv, "-cpt"))) {
    … …  
    118146  }
    119147
    120   if (argc != 2) usage_setastrom();
     148  if (!KH_FILE && !DCR_FILE) {
     149    fprintf (stderr, "at least one of -KH and -DCR must be supplied\n");
     150    exit (1);
     151  }
     152
     153  if (argc != 1) usage_setastrom();
    121154
    122155  return (TRUE);
    … …  
    124157
    125158void usage_setastrom_client () {
    126   fprintf (stderr, "USAGE: setastrom_client (correction) -hostID (hostID) -catdir (catdir) -hostdir (hostdir) [options]\n");
     159  fprintf (stderr, "USAGE: setastrom_client -hostID (hostID) -catdir (catdir) -hostdir (hostdir) [options]\n");
    127160  fprintf (stderr, "  options:\n");
     161  fprintf (stderr, "    -DCR (file) : supply DCR correction splines\n");
     162  fprintf (stderr, "    -KH (file) : supply Koppenhoefer correction splines\n");
    128163  fprintf (stderr, "    -region Rmin Rmax Dmin Dmax\n");
    129164  fprintf (stderr, "    -update-catformat (format) : change database schema on output\n");
    … …  
    163198
    164199  int N;
     200
     201  KH_FILE = NULL;
     202  if ((N = get_argument (argc, argv, "-KH"))) {
     203    remove_argument (N, &argc, argv);
     204    KH_FILE = strcreate (argv[N]);
     205    remove_argument (N, &argc, argv);
     206  }
     207  KH_RESET = FALSE;
     208  if ((N = get_argument (argc, argv, "-KH-reset"))) {
     209    remove_argument (N, &argc, argv);
     210    KH_RESET = TRUE;
     211  }
     212
     213  DCR_FILE = NULL;
     214  if ((N = get_argument (argc, argv, "-DCR"))) {
     215    remove_argument (N, &argc, argv);
     216    DCR_FILE = strcreate (argv[N]);
     217    remove_argument (N, &argc, argv);
     218  }
     219  DCR_RESET = FALSE;
     220  if ((N = get_argument (argc, argv, "-DCR-reset"))) {
     221    remove_argument (N, &argc, argv);
     222    DCR_RESET = TRUE;
     223  }
    165224
    166225  /* specify portion of the sky : allow default of all sky? */
    … …  
    223282  if (!CATDIR) usage_setastrom_client();
    224283
    225   if (argc != 2) usage_setastrom_client();
     284  if (!KH_FILE && !DCR_FILE) {
     285    fprintf (stderr, "at least one of -KH and -DCR must be supplied\n");
     286    exit (1);
     287  }
     288
     289  if (argc != 1) usage_setastrom_client();
    226290  return (TRUE);
    227291}

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/src/update_catalog_setastrom.c

    r37685 r37744  
    11# include "setastrom.h"
     2
     3// this code does both DCR and KH corrections, if both are specified. 
    24
    35// XXX this function is dependent in two ways on the specific numerical values of the photcodes
    … …  
    79int update_catalog_setastrom (Catalog *catalog) {
    810
    9   off_t i, found;
     11  off_t i;
    1012
    1113  time_t timeRef = ohana_date_to_sec("2011/05/11,00:00:00");
    1214
    13   found = 0;
     15  int Nsecfilt = GetPhotcodeNsecfilt ();
     16  PhotCode *code_g = GetPhotcodebyName ("g");
     17  PhotCode *code_i = GetPhotcodebyName ("i");
     18  PhotCode *code_z = GetPhotcodebyName ("z");
     19  PhotCode *code_y = GetPhotcodebyName ("y");
     20
     21  int Nsec_g = GetPhotcodeNsec(code_g->code);
     22  int Nsec_i = GetPhotcodeNsec(code_i->code);
     23  int Nsec_z = GetPhotcodeNsec(code_z->code);
     24  int Nsec_y = GetPhotcodeNsec(code_y->code);
     25
     26  off_t found = 0;
    1427
    1528  for (i = 0; i < catalog[0].Nmeasure; i++) {
    … …  
    1730    Measure *measure = &catalog[0].measure[i];
    1831
    19     // XXX only do GPC1 data for now
     32    // only do GPC1 data for now (true for both DCR and KH)
    2033    if (measure[0].photcode < 10000) continue;
    2134    if (measure[0].photcode > 10600) continue;
    2235     
    23     // XXX hardwire the list of chips to correct
     36    int averef = measure->averef;
     37    Average *average = &catalog[0].average[averef];
     38    SecFilt *secfilt = &catalog[0].secfilt[averef*Nsecfilt];
     39
     40    int chipID = measure[0].photcode % 100;
     41
     42    // XXX hardwire the list of chips to KH correct
    2443    // some notes:
    2544    // XY24 surprisingly does not need correction
    … …  
    2847    // XY67 has poor astrometry with a weird correction for all mags
    2948
    30     int chipID = measure[0].photcode % 100;
     49    int doKH = FALSE;
    3150    switch (chipID) {
    3251      case 4:
    … …  
    5776      case 72:
    5877      case 73:
     78        doKH = TRUE;
    5979        break;
    6080      default:
    61         measure[0].Xfix = measure[0].Xccd;
    62         measure[0].Yfix = measure[0].Yccd;
    63         continue;
    64     }
    65 
    66     double Xccd = measure[0].Xccd;
    67     double Yccd = measure[0].Yccd;
     81        break;
     82    }
     83
    6884    double pltScale = fabs(measure[0].pltscale);
    69 
    70     // this is not needed, is it?
    71     // double posAngle = FromShortDegrees(measure[0].posangle);
    72 
    73     // I need to know which chip we have.  use photcode assuming range
    74     int chip = measure[0].photcode % 100;
    75 
    76     // correction may be Minst or Minst + 5.0*log(fwhm
    77     float fwhm_maj = FromShortPixels(measure[0].FWx);
    78     float Minst = PhotInst (measure, MAG_CLASS_PSF);
    79     float MinstSB = Minst + 5.0*log10(fwhm_maj);
    80 
    81     // I should validate that the spline I'm using is SB corrected not Minst...
    82 
    83     // XXX this is very ad hoc.  the splines are measured based on an extraction from dvo
    84     // the dvo values of fwhm_maj are sigmas (NOT FWHM!) in pixels. 
     85    double posAngle = FromShortDegrees(measure[0].posangle);
     86
     87    /**** KH section ****/
     88
     89    double dX_KH = 0.0;
     90    double dY_KH = 0.0;
     91    if (KH_FILE && doKH && (measure[0].t < timeRef)) {
     92
     93      // correction may be Minst or Minst + 5.0*log(fwhm
     94      float fwhm_maj = FromShortPixels(measure[0].FWx);
     95      float Minst = PhotInst (measure, MAG_CLASS_PSF);
     96      // float MinstSB = Minst + 5.0*log10(fwhm_maj);
     97
     98      // NOTE: the spline correction was measured on the nightly-science database using measurements
     99      // from the period in which fwhm_maj & fwhm_min were a factor of sqrt(2) too small.  I am making
     100      // this correction here (but ideally I would re-measure the KH trend from PV2 data and re-generate
     101      // the splines
     102      float MinstSB = Minst + 5.0*log10(M_SQRT2*fwhm_maj);
     103
     104      // I should validate that the spline I'm using is SB corrected not Minst..
     105      // NOTE: it is SB, not raw Minst
     106
     107      // XXX this is very ad hoc.  the splines are measured based on an extraction from dvo
     108      // the dvo values of fwhm_maj are sigmas (NOT FWHM!) in pixels. 
    85109   
    86     // saturate fwhm_maj in the range 1.0 < fwhm_maj < 5.0
    87     // (images with really odd psf probably have poor astrometry anyway...)
    88     fwhm_maj = MIN(fwhm_maj,5.0);
    89     fwhm_maj = MAX(fwhm_maj,1.0);
    90 
    91     double dX;
    92     double dY;
    93 
    94     // correction is in arcseconds
    95     if (measure[0].t < timeRef) {
    96       get_astrom_correction (0, chip, &dX, &dY, MinstSB);
    97     } else {
    98       // XXX do not correct new data
    99       // get_astrom_correction (1, chip, &dX, &dY, MinstSB);
    100       dX = 0.0;
    101       dY = 0.0;
    102     }
     110      // saturate fwhm_maj in the range 1.0 < fwhm_maj < 5.0
     111      // (images with really odd psf probably have poor astrometry anyway...)
     112      fwhm_maj = MIN(fwhm_maj,5.0);
     113      fwhm_maj = MAX(fwhm_maj,1.0);
     114
     115      // correction is in arcseconds
     116      get_kh_correction (0, chipID, &dX_KH, &dY_KH, MinstSB);
     117    }
     118
     119    /**** DCR section ****/
     120
     121    double dX_DCR = 0.0;
     122    double dY_DCR = 0.0;
     123    if (DCR_FILE) {
     124
     125      // check if the color*airmass term is not-nan (Skip otherwise)
     126
     127      // XXX hard-wire gpc1 photcodes for now:
     128      // I should use photcode ops: PhotCode *code = PhotCodeByCode(measure->photcode)
     129      // int equiv = code->equiv
     130      // Nsec = GetNsec(equiv)
     131
     132      float dColor = 0.0;
     133      int filtCode = (int)(measure->photcode / 100); // eg, 101 = r
     134      switch (filtCode) {
     135        case 100:
     136        case 101:
     137        case 102: {
     138          float gmag = secfilt[Nsec_g].M;
     139          float imag = secfilt[Nsec_i].M;
     140          dColor = average->refColorBlue - (gmag - imag);
     141          break;
     142        }
     143        case 103:
     144        case 104: {
     145          float zmag = secfilt[Nsec_z].M;
     146          float ymag = secfilt[Nsec_y].M;
     147          dColor = average->refColorRed - (zmag - ymag);
     148          break;
     149        }
     150      }
     151      int filtSeq = filtCode % 100;
     152      if (!isfinite(dColor)) goto skip_DCR;
     153
     154      // I need to get the parallactic angle, but I only have alt & az, s:o
     155   
     156      // airmass = 1.0 / cos(90 - alt)
     157      // cos(90 - alt) = cos(zd) = sin(alt) = 1.0 / airmass
     158      // airmass = 1.0 / sin(alt)
     159
     160      // airmass = 1.0 / cos(zd)
     161      // cos(zd) = 1.0 / airmass == sn_alt
     162      // sin(zd) = sqrt(1.0 - cos(zd)^2) == cs_alt
     163      // tan(zd) = sin(zd) / cos(zd) = cs_alt / sn_alt
     164
     165      double sn_alt = 1.0 / measure->airmass;
     166      double cs_alt = sqrt(1.0 - SQ(sn_alt));
     167      double tan_zd = cs_alt / sn_alt;
     168
     169      double cs_az  = cos(RAD_DEG*measure->az);
     170      double sn_az  = sin(RAD_DEG*measure->az);
     171
     172# define LATITUDE (20.7070999146*RAD_DEG)
     173      double sn_lat = sin(LATITUDE);
     174      double cs_lat = cos(LATITUDE);
     175
     176      // double sind = sn_alt * sn_lat + cs_alt * cs_az * cs_lat;
     177      // *dec  = DEG_RAD * asin (sind);
     178
     179      // WARNING: these are NOT sin(ha),cos(ha) but are scaled
     180      // double ha = DEG_RAD * atan2 (sinh, cosh);
     181      double sinh = -cs_alt * sn_az;
     182      double cosh =  sn_alt * cs_lat - cs_alt * cs_az * sn_lat;
     183      double r_ha = hypot(sinh, cosh);
     184      double sn_ha = sinh / r_ha;
     185      double cs_ha = cosh / r_ha;
     186
     187      sinh = -cs_az * sn_alt * sn_ha * sn_lat + sn_az * sn_alt * cs_ha - sn_ha * cs_alt * cs_lat;
     188      cosh = -sn_az          * sn_ha * sn_lat - cs_az          * cs_ha;
     189      double r_rot = hypot(sinh, cosh);
     190
     191      double sn_rot = -sinh / r_rot;
     192      double cs_rot = +cosh / r_rot;
     193
     194      double dPx = 0.0;
     195      double dPy = 0.0;
     196
     197      float colorOffset = dColor * tan_zd;
     198      get_dcr_correction (filtSeq, &dPx, &dPy, colorOffset);
     199
     200      // if we had x-terms, we would have:
     201      // double dPx = dRx *cs_rot + dD_s*sn_rot;
     202      // double dPy = dD_s*cs_rot - dRx *sn_rot;
     203      // double dR = dPx*cs_rot - dPy*sn_rot;
     204      // double dD = dPx*sn_rot + dPy*cs_rot;
     205
     206      // since dPx = 0.0, instead we have: dPx, dPy as measured in dvo.dcr.sh
     207      double dR = -dPy*sn_rot;
     208      double dD = +dPy*cs_rot;
     209
     210      // dR = (measure->R - average->R) * 3600 * dcos(dec) [ie, in the projected plane]
     211      // dD = (measure->D - average->D) * 3600
     212      // average->D = measure->D - dD/3600
     213
     214      // now we need to rotate dR,dD into the chip frame and correction Xccd,Yccd
     215
     216      // are the plate-scale and posangle correctly set in the db?
     217
     218      // XXX the koppenhoefer code implies pos if angle from y -> D, (ie negative of the
     219      // below).  if so, apply sin(pos) -> -sin(pos)
     220
     221      // if parity = sky:
     222      // double dR =  dX*cos(pos) + dY*sin(pos);
     223      // double dD = -dX*sin(pos) + dY*cos(pos);
     224
     225      // XXX signs set here to reflect pos = -pos
     226      float posAngRad = posAngle * RAD_DEG;
     227      dX_DCR =  dR*cos(posAngRad) + dD*sin(posAngRad);
     228      dY_DCR = -dR*sin(posAngRad) + dD*cos(posAngRad);
     229    }
     230
     231  skip_DCR:
    103232
    104233    // do not modify the original Xccd,Yccd values
    105     measure[0].XoffKH = Xccd + dX / pltScale;
    106     measure[0].YoffKH = Yccd + dY / pltScale;
    107     measure[0].Xfix   = measure[0].XoffKH;
    108     measure[0].Yfix   = measure[0].YoffKH;
    109 
     234    if (KH_FILE) {
     235      measure[0].XoffKH = dX_KH / pltScale;
     236      measure[0].YoffKH = dY_KH / pltScale;
     237    }
     238    if (KH_RESET) {
     239      measure[0].XoffKH = 0.0;
     240      measure[0].YoffKH = 0.0;
     241    }
     242    if (DCR_FILE) {
     243      measure[0].XoffDCR = -dX_DCR / pltScale;
     244      measure[0].YoffDCR = -dY_DCR / pltScale;
     245    }
     246    if (DCR_RESET) {
     247      measure[0].XoffDCR = -dX_DCR / pltScale;
     248      measure[0].YoffDCR = -dY_DCR / pltScale;
     249    }
     250
     251    measure[0].Xfix = measure[0].Xccd;
     252    measure[0].Yfix = measure[0].Yccd;
     253    if (isfinite(measure[0].XoffKH) && isfinite(measure[0].YoffKH)) {
     254      measure[0].Xfix += measure[0].XoffKH;
     255      measure[0].Yfix += measure[0].YoffKH;
     256    }
     257    if (isfinite(measure[0].XoffDCR) && isfinite(measure[0].YoffDCR)) {
     258      measure[0].Xfix += measure[0].XoffDCR;
     259      measure[0].Yfix += measure[0].YoffDCR;
     260    }
    110261    found ++;
    111262  }

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/src/update_catalog_setastrom_dcr.c

    r37732 r37744  
    11# include "setastrom.h"
    22
    3 // XXX this function is dependent in two ways on the specific numerical values of the photcodes
    4 // 1) the examined measurements are limited to GPC1 values based on the photcode range 10000 - 10600
    5 // 2) the chip ID is determined from the numerical value of the photcode (chip = photcode % 100)
     3// this DCR code is largely hard-wired for gpc1 data
    64
    75int update_catalog_setastrom_dcr (Catalog *catalog) {
    … …  
    2422    Measure *measure = &catalog[0].measure[i];
    2523
    26     int averef = measure->averef;
    27     Average *average = &catalog[0].average[averef];
    28     Secfilt *secfilt = &catalog[0].secfilt[averef*Nsecfilt];
    29 
    30     int chipID = measure[0].photcode % 100;
    31 
    3224    double Xccd = measure[0].Xccd;
    3325    double Yccd = measure[0].Yccd;
    … …  
    3527    double posAngle = FromShortDegrees(measure[0].posangle);
    3628
    37     // I need to get the parallactic angle, but I only have alt & az, no
     29    // I need to get the parallactic angle, but I only have alt & az, s:o
    3830   
    3931    // airmass = 1.0 / cos(90 - alt)
    40     // cos(90 - alt) = 1.0 / airmass
    41     // cos(90)cos(alt) + sin(90)sin(alt)
    42     // 0*cos(alt) + sin(alt)
     32    // cos(90 - alt) = cos(zd) = sin(alt) = 1.0 / airmass
    4333    // airmass = 1.0 / sin(alt)
    4434
    … …  
    5646
    5747# define LATITUDE (20.7070999146*RAD_DEG)
     48
    5849    double sn_ha = -cs_alt * sn_az;
    5950    double cs_ha =  sn_alt * cos(LATITUDE) - cs_alt * cs_az * sin(LATITUDE);
    … …  
    6657
    6758    // XXX hard-wire gpc1 photcodes for now:
     59    // I should use photcode ops: PhotCode *code = PhotCodeByCode(measure->photcode)
     60    // int equiv = code->equiv
     61    // Nsec = GetNsec(equiv)
     62
    6863    int filtCode = (int)(measure->photcode / 100); // eg, 101 = r
    6964    switch (filtCode) {
    … …  
    8277        break;
    8378    }
     79    int filtSeq = filtCode % 100;
     80    float colorOffset = dColor * tan_dz;
    8481
    85     double dPy = dColor * tan_zd * DCRslope;
     82    get_dcr_correction (filtSeq, double *dPx, double *dPy, colorOffset);
    8683
    8784    // if we had x-terms, we would have:
    … …  
    9188    // double dD = dPx*sn_rot + dPy*cs_rot;
    9289
    93     // since dPx = 0.0, instead we have:
     90    // since dPx = 0.0, instead we have: dPx, dPy as measured in dvo.dcr.sh
    9491    double dR = -dPy*sn_rot;
    9592    double dD = +dPy*cs_rot;
    … …  
    107104
    108105    // if parity = sky:
    109     double dD =  dY*cos(pos) - dX*sin(pos);
    110     double dR =  dY*sin(pos) + dX*cos(pos);
     106    // double dR =  dX*cos(pos) + dY*sin(pos);
     107    // double dD = -dX*sin(pos) + dY*cos(pos);
     108
    111109    double dX =  dR*cos(pos) - dD*sin(pos);
    112110    double dY =  dR*sin(pos) + dD*cos(pos);
    113111
    114112    // if parity != sky:
    115     double dD =  dY*cos(pos) + dX*sin(pos);
    116     double dR =  dY*sin(pos) - dX*cos(pos);
    117     double dX = -dR*cos(pos) + dD*sin(pos);
    118     double dY =  dR*sin(pos) + dD*cos(pos);
     113    // double dD =  dY*cos(pos) + dX*sin(pos);
     114    // double dR =  dY*sin(pos) - dX*cos(pos);
     115    // double dX = -dR*cos(pos) + dD*sin(pos);
     116    // double dY =  dR*sin(pos) + dD*cos(pos);
    119117
    120118    // do not modify the original Xccd,Yccd values
     119
    121120    // why is KH added but DCR subtracted??
     121
     122    // ANSWER: both KH and DCR trends are measured by extracted dR,dD from the database
     123    // (mextract dR dD).  But the KH analysis was done before measure carried R,D when it
     124    // instead carried dR,dD.  In the current definition in dbExtractMeasures, dR,dD =
     125    // (measure.R,D - average.R,D) * 3600, but in the old format, measure.dR,dD were
     126    // defined as (average.R,D - measure.R,D)*3600).
     127
    122128    measure[0].XoffDCR = -dX / pltScale;
    123129    measure[0].YoffDCR = -dY / pltScale;

  • branches/eam_branches/ipp-20140904/Ohana/src/uniphot/src/update_dvo_setastrom.c

    r35103 r37744  
    11# include "setastrom.h"
    22
    3 int update_dvo_setastrom (char *corrfile) {
     3int update_dvo_setastrom () {
    44
    55  SkyTable *sky = NULL;
    … …  
    1313 
    1414  if (PARALLEL && !HOST_ID) {
    15     int status = update_dvo_setastrom_parallel (sky, corrfile);
     15    int status = update_dvo_setastrom_parallel (sky);
    1616    return status;
    1717  }
    1818
    19   load_astrom_correction (corrfile);
     19  if (KH_FILE)  load_kh_correction (KH_FILE);
     20  if (DCR_FILE) load_dcr_correction (DCR_FILE);
    2021
    2122  // determine the populated SkyRegions overlapping the requested area (default depth)
    … …  
    7879# define DEBUG 1
    7980
    80 int update_dvo_setastrom_parallel (SkyTable *sky, char *corrfile) {
     81int update_dvo_setastrom_parallel (SkyTable *sky) {
    8182
    8283  // now launch the setastrom_client jobs to the parallel hosts
    83 
    84   char *abscorrfile = abspath (corrfile, DVO_MAX_PATH);
    8584
    8685  // load the list of hosts
    … …  
    9695
    9796    char command[1024];
    98     snprintf (command, 1024, "setastrom_client %s -hostID %d -catdir %s -hostdir %s -region %f %f %f %f",
    99               abscorrfile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname,
     97    snprintf (command, 1024, "setastrom_client -hostID %d -catdir %s -hostdir %s -region %f %f %f %f",
     98              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname,
    10099              UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
    101 
    102     fprintf (stderr, "command: %s\n", command);
    103100
    104101    char tmpline[1024];
    105102    if (VERBOSE)          { snprintf (tmpline, 1024, "%s -v",                   command);                   strcpy (command, tmpline); }
    106103    if (UPDATE)           { snprintf (tmpline, 1024, "%s -update",              command);                   strcpy (command, tmpline); }
     104    if (KH_FILE)          { snprintf (tmpline, 1024, "%s -KH %s",               command, KH_FILE);          strcpy (command, tmpline); }
     105    if (DCR_FILE)         { snprintf (tmpline, 1024, "%s -DCR %s",              command, DCR_FILE);         strcpy (command, tmpline); }
    107106    if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT); strcpy (command, tmpline); }
     107
     108    if (KH_RESET)         { snprintf (tmpline, 1024, "%s -KH-reset",            command);                   strcpy (command, tmpline); }
     109    if (DCR_RESET)        { snprintf (tmpline, 1024, "%s -DCR-reset",           command);                   strcpy (command, tmpline); }
     110
     111    fprintf (stderr, "command: %s\n", command);
    108112
    109113    if (PARALLEL_MANUAL) continue;
    … …  
    126130    }
    127131  }
    128   free (abscorrfile);
    129132
    130133  // wait for the remote jobs to be completed
Note: See TracChangeset for help on using the changeset viewer.