IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34868


Ignore:
Timestamp:
Dec 21, 2012, 12:55:14 PM (14 years ago)
Author:
eugene
Message:

setastrom builds

Location:
branches/eam_branches/ipp-20121219/Ohana/src/uniphot
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

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

    r34867 r34868  
    118118$(SRC)/update_catalog_setastrom.$(ARCH).o \
    119119$(SRC)/astrom_correction.$(ARCH).o \
    120 $(SRC)/initialize_setastrom_client.$(ARCH).o
     120$(SRC)/initialize_setastrom.$(ARCH).o
    121121
    122122$(SETASTROM_CLIENT): $(INC)/setastrom.h

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

    r34867 r34868  
    22# include <dvo.h>
    33# include <signal.h>
     4
     5// the interpolating spline information: this is based on the spline used by opihi/cmd.data/spline.c
     6typedef struct {
     7  int Nknots;
     8  double *xk;
     9  double *yk;
     10  double *y2;
     11} Spline;
    412
    513/* global variables set in parameter file */
    … …  
    1422int          PARALLEL_MANUAL;
    1523int          PARALLEL_SERIAL;
     24char        *UPDATE_CATFORMAT;
     25
     26SkyRegion    UserPatch;
    1627
    1728/***** prototypes ****/
    1829int           main                            PROTO((int argc, char **argv));
    1930
     31void          GetConfig                       PROTO((char *config, char *field, char *format, int N, void *ptr));
    2032void          ConfigInit                      PROTO((int *argc, char **argv));
     33
    2134void          usage_setastrom                 PROTO(());
    2235void          initialize_setastrom            PROTO((int argc, char **argv));
    … …  
    3245// int           SetSignals                      PROTO((void));
    3346
    34 int           update_dvo_setastrom            PROTO(());
    35 int           update_dvo_setastrom_parallel   PROTO((SkyTable *sky));
     47int           update_dvo_setastrom            PROTO((char *corrfile));
     48int           update_dvo_setastrom_parallel   PROTO((SkyTable *sky, char *corrfile));
    3649
    37 void          update_catalog_setastrom        PROTO((Catalog *catalog));
     50int           update_catalog_setastrom        PROTO((Catalog *catalog));
    3851
     52// spline correction functions
     53int load_astrom_correction (char *filename);
     54int get_astrom_correction (int sub, int chip, double *dX, double *dY, float Minst);
     55double spline_apply_dbl (double *x, double *y, double *y2, int N, double X);

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

    r34867 r34868  
    1313// these are stored in an N-D array of (Spline *) splineset[sub][dir][chip]
    1414
    15 // the interpolating spline information: this is based on the spline used by opihi/cmd.data/spline.c
    16 typedef struct {
    17   int Nknots;
    18   double *xk;
    19   double *yk;
    20   double *y2;
    21   char *name;
    22 } Spline;
    23 
    2415# define N_SUB 2
    2516# define N_DIR 2
    … …  
    3021
    3122int load_astrom_correction (char *filename) {
     23
     24  int Ncol;
     25  off_t Nrow;
     26
     27  char type[16];
     28  char extname[80];
    3229
    3330  Header header;
    … …  
    5451
    5552  // we have a N-D array of splines:
    56   // splineset[subset][dir][chip]
    57   ALLOCATE (splineset, Spline ***, 2);
    58   for (sub = 0; sub < 2; sub++) {
    59     ALLOCATE (splineset[sub], Spline **, 2);
    60     for (dir = 0; dir < 2; dir++) {
    61       ALLOCATE (splineset[sub][dir], Spline *, 78);
    62       for (chip = 0; chip < 78; chip ++) {
     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 ++) {
    6361        splineset[sub][dir][chip] = NULL;
    6462      }
    … …  
    7068    if (!gfits_load_header (f, &theader)) break;
    7169   
    72     if (!gfits_fread_ftable_data (f, &table, FALSE)) return (FALSE);
     70    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (FALSE);
    7371   
    74     if (!gfits_scan (header, "EXTNAME", "%s", 1, extname)) return (FALSE);
     72    if (!gfits_scan (&theader, "EXTNAME", "%s", 1, extname)) return (FALSE);
    7573
    7674    // figure out the chip ID, diretion, and sequence from the EXTNAME
    7775
    78     int chip = XX;
    79     int dir = XX;
    80     int subset = XX;
     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);
    8186
    82     splineset[sub][dir][chip] = CreateSpline (Nrow);
     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);
    8395
    84     // XXX: need to handle case of spline data not existing...
    85    
     96    int sub = -1;
     97    if (extname[6] == '0') {
     98      sub = 0;
     99    }
     100    if (extname[6] == '1') {
     101      sub = 1;
     102    }
     103    myAssert (sub != -1, "error in extname %s (sub)\n", extname);
     104
    86105    // need to create and assign to flat-field correction
    87106    double *xk = gfits_get_bintable_column_data (&theader, &ftable, "X_KNOT", type, &Nrow, &Ncol);
    88     assert (!strcmp(type, "double"));
     107    myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    89108   
    90109    // need to create and assign to flat-field correction
    91110    double *yk = gfits_get_bintable_column_data (&theader, &ftable, "Y_KNOT", type, &Nrow, &Ncol);
    92     assert (!strcmp(type, "double"));
     111    myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    93112   
    94113    // need to create and assign to flat-field correction
    95114    double *y2 = gfits_get_bintable_column_data (&theader, &ftable, "DY2_DX", type, &Nrow, &Ncol);
    96     assert (!strcmp(type, "double"));
     115    myAssert (!strcmp(type, "double"), "error reading X_KNOT");
    97116   
     117    ALLOCATE (splineset[sub][dir][chip], Spline, 1);
    98118    splineset[sub][dir][chip]->xk = xk;
    99119    splineset[sub][dir][chip]->yk = yk;
    100120    splineset[sub][dir][chip]->y2 = y2;
    101    
     121    splineset[sub][dir][chip]->Nknots = Nrow;
    102122  }
    103123  return (TRUE);
    104124}
     125
     126int 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  spline = splineset[sub][0][chip];
     134  if (spline) {
     135    *dX = spline_apply_dbl (spline->xk, spline->yk, spline->y2, spline->Nknots, Minst);
     136  }
     137
     138  spline = splineset[sub][1][chip];
     139  if (spline) {
     140    *dY = spline_apply_dbl (spline->xk, spline->yk, spline->y2, spline->Nknots, Minst);
     141  }
     142  return (TRUE);
     143}
     144
     145/* evaluate spline for x, y, y2 at X */
     146double spline_apply_dbl (double *x, double *y, double *y2, int N, double X) {
     147
     148  int i, lo, hi;
     149  double dx, a, b, value;
     150 
     151  /* find correct element in array (x must be sorted) */
     152  lo = 0;
     153  hi = N-1;
     154  while (hi - lo > 1) {
     155    i = 0.5*(hi+lo);
     156    if (x[i] > X) {
     157      hi = i;
     158    } else {
     159      lo = i;
     160    }
     161  }
     162
     163  /* error condition: duplicate abssisca */
     164  dx = x[hi] - x[lo];
     165  if (dx == 0.0) {
     166    return (HUGE_VAL);
     167  }
     168
     169  /* evaluate spline */
     170  a = (x[hi] - X) / dx;
     171  b = (X - x[lo]) / dx;
     172
     173  value = a*y[lo] + b*y[hi] + ((a*a*a - a)*y2[lo] + (b*b*b - b)*y2[hi])*(dx*dx) / 6.0;
     174  return (value);
     175
     176}

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

    r34867 r34868  
    22
    33void usage_setastrom () {
    4     fprintf (stderr, "USAGE: setastrom [options]\n");
     4    fprintf (stderr, "USAGE: setastrom (corrections) [options]\n");
    55    fprintf (stderr, "  options:\n");
    66    fprintf (stderr, "    -v : verbose mode\n");
     7    fprintf (stderr, "    -region Rmin Rmax Dmin Dmax\n");
     8    fprintf (stderr, "    -update-catformat (format) : change database schema on output\n");
    79    fprintf (stderr, "    -update : actually write results to detections tables\n");
    810    fprintf (stderr, "    -parallel : run in parallel mode\n");
    … …  
    4648
    4749  int N;
     50
     51  /* specify portion of the sky : allow default of all sky? */
     52  UserPatch.Rmin = 0;
     53  UserPatch.Rmax = 360;
     54  UserPatch.Dmin = -90;
     55  UserPatch.Dmax = +90;
     56  if ((N = get_argument (argc, argv, "-region"))) {
     57    remove_argument (N, &argc, argv);
     58    UserPatch.Rmin = atof (argv[N]);
     59    remove_argument (N, &argc, argv);
     60    UserPatch.Rmax = atof (argv[N]);
     61    remove_argument (N, &argc, argv);
     62    UserPatch.Dmin = atof (argv[N]);
     63    remove_argument (N, &argc, argv);
     64    UserPatch.Dmax = atof (argv[N]);
     65    remove_argument (N, &argc, argv);
     66  }
     67
     68  UPDATE_CATFORMAT = NULL;
     69  if ((N = get_argument (argc, argv, "-update-catformat"))) {
     70    remove_argument (N, &argc, argv);
     71    UPDATE_CATFORMAT = strcreate (argv[N]);
     72    remove_argument (N, &argc, argv);
     73  }
    4874
    4975  VERBOSE = FALSE;
    … …  
    85111  }
    86112
    87   if (argc != 1) usage_setastrom();
     113  if (argc != 2) usage_setastrom();
    88114
    89115  return (TRUE);
    … …  
    91117
    92118void usage_setastrom_client () {
    93   fprintf (stderr, "USAGE: setastrom_client -hostID (hostID) -catdir (catdir) -hostdir (hostdir) -images (images) [options]\n");
     119  fprintf (stderr, "USAGE: setastrom_client (correction) -hostID (hostID) -catdir (catdir) -hostdir (hostdir) [options]\n");
    94120  fprintf (stderr, "  options:\n");
     121  fprintf (stderr, "    -region Rmin Rmax Dmin Dmax\n");
     122  fprintf (stderr, "    -update-catformat (format) : change database schema on output\n");
    95123  fprintf (stderr, "    -v : verbose mode\n");
    96124  fprintf (stderr, "    -update : actually write results to detections tables\n");
    … …  
    129157  int N;
    130158
     159  /* specify portion of the sky : allow default of all sky? */
     160  UserPatch.Rmin = 0;
     161  UserPatch.Rmax = 360;
     162  UserPatch.Dmin = -90;
     163  UserPatch.Dmax = +90;
     164  if ((N = get_argument (argc, argv, "-region"))) {
     165    remove_argument (N, &argc, argv);
     166    UserPatch.Rmin = atof (argv[N]);
     167    remove_argument (N, &argc, argv);
     168    UserPatch.Rmax = atof (argv[N]);
     169    remove_argument (N, &argc, argv);
     170    UserPatch.Dmin = atof (argv[N]);
     171    remove_argument (N, &argc, argv);
     172    UserPatch.Dmax = atof (argv[N]);
     173    remove_argument (N, &argc, argv);
     174  }
     175
     176  UPDATE_CATFORMAT = NULL;
     177  if ((N = get_argument (argc, argv, "-update-catformat"))) {
     178    remove_argument (N, &argc, argv);
     179    UPDATE_CATFORMAT = strcreate (argv[N]);
     180    remove_argument (N, &argc, argv);
     181  }
     182
    131183  VERBOSE = FALSE;
    132184  if ((N = get_argument (argc, argv, "-v"))) {
    … …  
    164216  if (!CATDIR) usage_setastrom_client();
    165217
    166   IMAGES = NULL;
    167   if ((N = get_argument (argc, argv, "-images"))) {
    168     remove_argument (N, &argc, argv);
    169     IMAGES = strcreate (argv[N]);
    170     remove_argument (N, &argc, argv);
    171   }
    172   if (!IMAGES) usage_setastrom_client();
    173 
    174   if (argc != 1) usage_setastrom_client();
     218  if (argc != 2) usage_setastrom_client();
    175219  return (TRUE);
    176220}
    177 
    178 void GetConfig (char *config, char *field, char *format, int N, void *ptr);
    179221
    180222void ConfigInit (int *argc, char **argv) {

  • branches/eam_branches/ipp-20121219/Ohana/src/uniphot/src/setastrom.c

    r34866 r34868  
    44
    55  int status;
    6   FITS_DB db;
    76
    87  /* get configuration info, args, lockfile */
    98  initialize_setastrom (argc, argv);
    109
    11   load_correction ();
    12 
    13   status = update_dvo_setastrom (subset, Nimage);
     10  // argv[1] holds the correction file
     11  status = update_dvo_setastrom (argv[1]);
    1412
    1513  if (!status) exit (1);

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

    r34866 r34868  
    77int update_catalog_setastrom (Catalog *catalog) {
    88
    9   off_t i, j, found;
     9  off_t i, found;
    1010
    11   int timeRef = somefunc("2011/05/11");
     11  time_t timeRef = ohana_date_to_sec("2011/05/11,00:00:00");
    1212
    1313  for (i = 0; i < catalog[0].Nmeasure; i++) {
    … …  
    1919    double Xccd = catalog[0].measure[i].Xccd;
    2020    double Yccd = catalog[0].measure[i].Yccd;
    21     double posAngle = FromShortDegrees(catalog[0].measure[i].posangle);
    2221    double pltScale = fabs(catalog[0].measure[i].pltscale);
     22
     23    // this is not needed, is it?
     24    // double posAngle = FromShortDegrees(catalog[0].measure[i].posangle);
    2325
    2426    // I need to know which chip we have.  use photcode assuming range
    2527    int chip = catalog[0].measure[i].photcode % 100;
    2628
     29    float Minst = PhotInst (&catalog[0].measure[i]);
     30
     31    double dX;
     32    double dY;
     33
    2734    // correction is in arcseconds
    2835    if (catalog[0].measure[i].t < timeRef) {
    29       get_correction (0, chip, &dX, &dY);
     36      get_astrom_correction (0, chip, &dX, &dY, Minst);
    3037    } else {
    31       get_correction (1, chip, &dX, &dY);
     38      get_astrom_correction (1, chip, &dX, &dY, Minst);
    3239    }
    3340

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

    r34866 r34868  
    11# include "setastrom.h"
    22
    3 int update_dvo_setastrom () {
     3int update_dvo_setastrom (char *corrfile) {
    44
    5   SkyRegion UserPatch;
    65  SkyTable *sky = NULL;
    76  SkyList *skylist = NULL;
    87  Catalog catalog;
    9   off_t i, maxID, *index;
     8  off_t i;
    109
    1110  // load the current sky table (layout of all SkyRegions)
    … …  
    1413 
    1514  if (PARALLEL && !HOST_ID) {
    16     update_dvo_setastrom_parallel (sky);
    17     return TRUE;
     15    int status = update_dvo_setastrom_parallel (sky, corrfile);
     16    return status;
    1817  }
    1918
     19  load_astrom_correction (corrfile);
     20
    2021  // determine the populated SkyRegions overlapping the requested area (default depth)
    21   UserPatch.Rmin = 0;
    22   UserPatch.Rmax = 360;
    23   UserPatch.Dmin = -90;
    24   UserPatch.Dmax = +90;
    2522  skylist = SkyListByPatch (sky, -1, &UserPatch);
    2623
    … …  
    6057    }
    6158   
    62   XXX: update format here
     59    // modify the output format as desired (ignore current format on disk)
     60    if (UPDATE_CATFORMAT) {
     61      catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT);
     62    }
    6363
    6464    if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename);
    … …  
    7373# define DEBUG 1
    7474
    75 int update_dvo_setastrom_parallel (SkyTable *sky) {
     75int update_dvo_setastrom_parallel (SkyTable *sky, char *corrfile) {
    7676
    7777  // now launch the setastrom_client jobs to the parallel hosts
     78
     79  char *abscorrfile = abspath (corrfile, DVO_MAX_PATH);
    7880
    7981  // load the list of hosts
    … …  
    8991
    9092    char command[1024];
    91     snprintf (command, 1024, "setastrom_client -hostID %d -catdir %s -hostdir %s",
    92               table->hosts[i].hostID, CATDIR, table->hosts[i].pathname);
     93    snprintf (command, 1024, "setastrom_client %s -hostID %d -catdir %s -hostdir %s -region %f %f %f %f",
     94              abscorrfile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname,
     95              UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
    9396
    9497    fprintf (stderr, "command: %s\n", command);
    9598
    9699    char tmpline[1024];
    97     if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
    98     if (UPDATE)        { snprintf (tmpline, 1024, "%s -update",         command);                    strcpy (command, tmpline); }
     100    if (VERBOSE)          { snprintf (tmpline, 1024, "%s -v",                   command);                   strcpy (command, tmpline); }
     101    if (UPDATE)           { snprintf (tmpline, 1024, "%s -update",              command);                   strcpy (command, tmpline); }
     102    if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT); strcpy (command, tmpline); }
    99103
    100104    if (PARALLEL_MANUAL) continue;
    … …  
    117121    }
    118122  }
     123  free (abscorrfile);
    119124
    120125  // wait for the remote jobs to be completed
Note: See TracChangeset for help on using the changeset viewer.