IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26384


Ignore:
Timestamp:
Dec 11, 2009, 1:57:01 PM (16 years ago)
Author:
bills
Message:

Several changes to support managing PSPS object and detction ids.

  1. Support for 64 bit columns in fits tables. "K"
  2. if TSCAL == 1 and TZERO == 0 omit those keywords from the tables
  3. compute psps object and detection ids
Location:
trunk/Ohana/src
Files:
12 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/Makefile

    r26280 r26384  
    8181$(SRC)/ReadImageHeader.$(ARCH).o \
    8282$(SRC)/UpdateImageIDs.$(ARCH).o \
    83 $(SRC)/update_coords.$(ARCH).o
     83$(SRC)/update_coords.$(ARCH).o \
     84$(SRC)/psps_ids.$(ARCH).o
    8485
    8586ADDSTARD = \
    … …  
    118119$(SRC)/opening_angle.$(ARCH).o \
    119120$(SRC)/replace_match.$(ARCH).o \
    120 $(SRC)/update_coords.$(ARCH).o
     121$(SRC)/update_coords.$(ARCH).o \
     122$(SRC)/psps_ids.$(ARCH).o
    121123
    122124ADDSTART = \
    … …  
    157159$(SRC)/NewImage_Thread.$(ARCH).o \
    158160$(SRC)/NewReflist_Thread.$(ARCH).o \
    159 $(SRC)/NewRefcat_Thread.$(ARCH).o
     161$(SRC)/NewRefcat_Thread.$(ARCH).o \
     162$(SRC)/psps_ids.$(ARCH).o
    160163
    161164ADDSTARC = \
    … …  
    206209$(SRC)/ConfigInit.$(ARCH).o \
    207210$(SRC)/Shutdown.$(ARCH).o \
    208 $(SRC)/SetSignals.$(ARCH).o
     211$(SRC)/SetSignals.$(ARCH).o \
     212$(SRC)/psps_ids.$(ARCH).o
    209213
    210214SEDSTAR = \

  • trunk/Ohana/src/addstar/include/addstar.h

    r26286 r26384  
    289289PhotCode *pmm_get_photcode (char *emulsion, char *filter);
    290290
     291#define PSPS_ID TRUE
     292uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
     293uint64_t CreatePSPSObjectID(double ra, double dec);
     294
    291295// this is a gnu extension?? caution!
    292296void *memrchr(const void *s, int c, size_t n);

  • trunk/Ohana/src/addstar/src/FilterStars.c

    r26338 r26384  
    8888    }
    8989
    90 # if (0)
    9190    if (PSPS_ID) {
    9291      double mjd;
    9392      mjd = ohana_sec_to_mjd (image[0].tzero);
    94       // XXX this is clearly wrong : what does PSPS want?
    95       stars[N].measure.extID = PSPS_create_detectid (mjd, imageID);
     93      stars[N].measure.extID = CreatePSPSDetectionID(mjd, image[0].ccdnum, stars[N].measure.detID);
     94    } else {
     95      stars[N].measure.extID = 0;
    9696    }
    97 # endif
    9897
    99     stars[N].measure.imageID = imageID; // this value is update in UpdateImageIDs
     98    stars[N].measure.imageID = imageID; // this value is updated in UpdateImageIDs
    10099
    101100    N ++;

  • trunk/Ohana/src/addstar/src/find_matches.c

    r26283 r26384  
    271271    catalog[0].average[Nave].catID         = catID;
    272272    catalog[0].average[Nave].flags         = 0;
     273    if (PSPS_ID) {
     274        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
     275                                                            catalog[0].average[Nave].D);
     276    } else {
     277        catalog[0].average[Nave].extID         = 0;
     278    }
    273279
    274280    objID ++;

  • trunk/Ohana/src/addstar/src/find_matches_closest.c

    r26283 r26384  
    266266    catalog[0].average[Nave].catID         = catID;
    267267    catalog[0].average[Nave].flags         = 0;
     268    if (PSPS_ID) {
     269        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
     270                                                            catalog[0].average[Nave].D);
     271    } else {
     272        catalog[0].average[Nave].extID         = 0;
     273    }
     274
    268275
    269276    objID ++;

  • trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c

    r26283 r26384  
    291291    catalog[0].average[Nave].catID         = catID;
    292292    catalog[0].average[Nave].flags         = 0;
     293    if (PSPS_ID) {
     294        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
     295                                                            catalog[0].average[Nave].D);
     296    } else {
     297        catalog[0].average[Nave].extID         = 0;
     298    }
    293299
    294300    objID ++;

  • trunk/Ohana/src/addstar/src/find_matches_refstars.c

    r26283 r26384  
    261261    catalog[0].average[Nave].catID         = catID;
    262262    catalog[0].average[Nave].flags         = 0;
     263    if (PSPS_ID) {
     264        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
     265                                                            catalog[0].average[Nave].D);
     266    } else {
     267        catalog[0].average[Nave].extID         = 0;
     268    }
     269
    263270
    264271    objID ++;

  • trunk/Ohana/src/addstar/src/psps_ids.c

    • Property svn:mergeinfo set to (toggle deleted branches)
      /branches/eam_branches/20091113/Ohana/src/addstar/src/PSPS_IDs.c26119-26255
    r26370 r26384  
    11# include "addstar.h"
    22
    3 uint64_7 PSPS_create_detectID (double MJDobs, int imageID) {
     3// Compute PSPS defined IDs
    44
    5   // MJD has a PS lifetime range of : 3288 (2009/01/01) to 9132 (2025/01/01) (< 2^14)
     5uint64_t
     6CreatePSPSDetectionID(double tobs, int ccdid, int detID)
     7{
     8//    double floor();
     9    // t0 for detection id is 2007-01-01 00:00:00 utc
     10    double  t0 = 54101.0;
     11    double diff;
     12    int itmp;
     13    uint64_t detectid;
     14       
     15    diff = floor( 100000. * (tobs - t0) );
     16    itmp = diff;
     17    // ccdid must be < 100
     18    detectid = 1000000000*((uint64_t) itmp) + 10000000 * ((uint64_t) ccdid) +
     19             ((uint64_t) detID);
    620
    7   static double t0 = 51544.5; // MJD @ 2000/01/01,12:00:00
     21    return detectid;
     22}
     23   
     24uint64_t
     25CreatePSPSObjectID(double ra, double dec)
     26{
     27    uint64_t part1, part2, part3;
     28    double zid, zresid;
     29    int izone;
     30    double zh = 0.0083333;
     31   
     32    zid = (dec + 90.) / zh;             // 0 - 180*60*2 = 21600 < 15 bits
     33    izone = (int) floor(zid);
     34    zresid = zid -  ((float) izone);    // 0.0 - 1.0
    835
    9   uint64_t iTime;
    10   uint64_t detectID;
     36    part1 = (uint64_t)( izone  * 10000000000000) ;
     37    part2 = ((uint64_t)(ra * 1000000.)) * 10000 ; // 0 - 360*1e6 = 3.6e8 (< 29 bits)
     38    part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec (< 14 bits)
    1139
    12   iTime = floor(tobs - t0);
    13 
    14   // iTime takes the upper 16 bits, imageID can take up to lower 48 bits
    15   detectID = (iTime << 48) | imageID;
    16 
    17   return detectID;
     40    return part1 + part2 + part3;
    1841}
    19 
    20 uint64_t PSPS_create_objID (double ra, double dec) {
    21 
    22   static double zh = 0.0083333; // zone height = 30 arcsec
    23 
    24   uint64_t objID;
    25   uint64_t part1, part2, part3;
    26 
    27   double zid, zresid;
    28   int izone;
    29 
    30   zid = (dec + 90.0) / zh; // 0 - 180*60*2 = 21600 (< 15 bits)
    31   izone = (int) floor(zid);
    32   zresid = zid -  ((float) izone); // 0 - 1.0
    33   part1 = (uint64_t)( izone  * 10000000000000) ;
    34   part2 = ((uint64_t)(ra * 1 000 000.)) * 10000 ; // 0 - 360.0*1e6 = 3.6e8 (< 29 bits)
    35   part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec = 0.003 arcsec) (< 14 bits)
    36 
    37   objID = part1 + part2 + part3;
    38 
    39   return objID;
    40 }
    41 
    42 uint64_t PSPS_create_objID_eam (double ra, double dec) {
    43 
    44   static double zh = 0.0083333; // zone height = 30 arcsec
    45 
    46   uint64_t objID;
    47   uint64_t part1, part2, part3;
    48 
    49   double zid, zresid;
    50   int izone, ira;
    51 
    52   zid = (dec + 90.0) / zh; // 0 - 180*60*2 = 21600 (< 15 bits)
    53   izone = (int) floor(zid);
    54   zresid = zid -  ((float) izone); // 0 - 1.0
    55 
    56   ira = ra * 1 000 000; // 0 - 360.0*1e6 = 3.6e8 (< 29 bits)
    57 
    58   part1 = (izone << 48);
    59   part2 = (ira << 16);
    60   part3 = (int) (zresid * 10000.0) ; // 0 - 10000 (1 bit == 30/10000 arcsec = 0.003 arcsec) (< 14 bits)
    61 
    62   objID = part1 | part2 | part3;
    63 
    64   return objID;
    65 }
    66 
    67 // part1 : 15 bits
    68 // part2 : 29 bits
    69 // part3 : 14 bits
    70 // 15 + 14 + 29 bits = 58 bits...
    71 

  • trunk/Ohana/src/libautocode/generate

    r21508 r26384  
    135135        # FITS tables do not allow for 64bit integers.  we need to
    136136        # write these by splitting the value into high and low 32 bit values
    137         if ($type eq "int64_t")       { $pt1 = "J"; $Np = 2*$Np; }
    138         if ($type eq "uint64_t")      { $pt1 = "J"; $Np = 2*$Np; }
     137        if ($type eq "int64_t")       { $pt1 = "K";}
     138        if ($type eq "uint64_t")      { $pt1 = "K";}
    139139
    140140        if ($type eq "float")         { $pt1 = "E"; }
    … …  
    239239        if ($type eq "uint32_t")      { $pt1 = sprintf "I%s", $length; }
    240240
    241         # FITS tables do not allow for 64bit integers.  we need to
    242         # write these by splitting the value into high and low 32 bit values
    243         if ($type eq "int64_t")       { $pt1 = sprintf "I%s", $length; }
    244         if ($type eq "uint64_t")      { $pt1 = sprintf "I%s", $length; }
     241        if ($type eq "int64_t")       { $pt1 = sprintf "K%s", $length; }
     242        if ($type eq "uint64_t")      { $pt1 = sprintf "K%s", $length; }
    245243
    246244        if ($type eq "float")         { $pt1 = sprintf "F%s", $length; }
    … …  
    324322        if ($type eq "uint32_t")      { $T = "WORD"; $n = 4; }
    325323
    326         # FITS tables do not allow for 64bit integers.  we need to
    327         # write these by splitting the value into high and low 32 bit values
    328         if ($type eq "int64_t")       { $T = "WORD"; $n = 4; $Np = 2*$Np; }
    329         if ($type eq "uint64_t")      { $T = "WORD"; $n = 4; $Np = 2*$Np; }
     324        if ($type eq "int64_t")       { $T = "DBLE"; $n = 8; }
     325        if ($type eq "uint64_t")      { $T = "DBLE"; $n = 8; }
    330326
    331327        if ($type eq "float")         { $T = "WORD"; $n = 4; }
    … …  
    370366        if ($type eq "int32_t")        { $Nbytes += 4*$Np; $valid = 1; }
    371367        if ($type eq "uint32_t")       { $Nbytes += 4*$Np; $valid = 1; }
    372 
    373         # FITS tables do not allow for 64bit integers.  we need to
    374         # write these by splitting the value into high and low 32 bit values
    375368        if ($type eq "int64_t")        { $Nbytes += 8*$Np; $valid = 1; }
    376369        if ($type eq "uint64_t")       { $Nbytes += 8*$Np; $valid = 1; }

  • trunk/Ohana/src/libfits/table/F_define_column.c

    r7054 r26384  
    2323  sprintf (field, "TFORM%d", Nfields);
    2424  gfits_modify (header, field, "%s", 1, format);
    25   sprintf (field, "TSCAL%d", Nfields);
    26   gfits_modify (header, field, "%lf", 1, bscale);
    27   sprintf (field, "TZERO%d", Nfields);
    28   gfits_modify (header, field, "%lf", 1, bzero);
     25
     26  // add scaling parameters unless they amount to a noop
     27  if ((bscale != 1.0) || (bzero != 0.0)) {
     28      sprintf (field, "TSCAL%d", Nfields);
     29      gfits_modify (header, field, "%lf", 1, bscale);
     30      sprintf (field, "TZERO%d", Nfields);
     31      gfits_modify (header, field, "%lf", 1, bzero);
     32  }
    2933
    3034  /* update TFIELDS & NAXIS1 */

  • trunk/Ohana/src/libfits/table/F_get_column.c

    r16810 r26384  
    11# include <ohana.h>
    22# include <gfitsio.h>
     3#include <inttypes.h>
    34# define SWAP_BYTE { \
    45  char tmp; \
    … …  
    9394    }
    9495  }
     96  if (!strcmp (type, "int64_t")) {
     97    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
     98# ifdef BYTE_SWAP
     99      SWAP_DBLE;
     100# endif
     101      *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero;
     102    }
     103  }
    95104  if (!strcmp (type, "float")) {
    96105    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
    … …  
    218227# endif
    219228      *(int *)Pout = *(int *)Pin*Bscale + Bzero;
     229    }
     230  }
     231  if (!strcmp (type, "int64_t")) {
     232    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
     233# ifdef BYTE_SWAP
     234      SWAP_DBLE;
     235# endif
     236      *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero;
    220237    }
    221238  }
    … …  
    349366    }
    350367  }
     368  if (!strcmp (type, "int64_t")) {
     369    ALLOCATE (array, char, Ny*8);
     370    Pout = array;
     371    for (i = 0; i < Ny; i++, Pin+=Nx, Pout+=8) {
     372      memcpy (line, Pin, Nval*Nbytes);
     373      sscanf (line, "%" PRId64, (int64_t *)Pout);
     374    }
     375  }
    351376  if (!strcmp (type, "float")) {
    352377    ALLOCATE (array, char, Ny*4);

  • trunk/Ohana/src/libfits/table/F_table_format.c

    r15487 r26384  
    1717
    1818  *Nbytes = 0;
    19   if (*Fchar == 'X') { *Nbytes = 1;  strcpy (type, "char");   *Nval = 1 + (int) Nv / 8; }
    20   if (*Fchar == 'L') { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
    21   if (*Fchar == 'A') { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
    22   if (*Fchar == 'B') { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
    23   if (*Fchar == 'I') { *Nbytes = 2;  strcpy (type, "short");  *Nval = Nv;               }
    24   if (*Fchar == 'J') { *Nbytes = 4;  strcpy (type, "int");    *Nval = Nv;               }
    25   if (*Fchar == 'E') { *Nbytes = 4;  strcpy (type, "float");  *Nval = Nv;               }
    26   if (*Fchar == 'D') { *Nbytes = 8;  strcpy (type, "double"); *Nval = Nv;               }
    27   if (*Fchar == 'P') { *Nbytes = 8;  strcpy (type, "var");    *Nval = 2*Nv;             }
    28   if (*Fchar == 'C') { *Nbytes = 8;  strcpy (type, "float");  *Nval = 2*Nv;             }
    29   if (*Fchar == 'M') { *Nbytes = 16; strcpy (type, "double"); *Nval = 2*Nv;             }
    30   if (!*Nbytes) { return (FALSE); }
     19  switch (*Fchar) {
     20  case  'X':
     21    { *Nbytes = 1;  strcpy (type, "char");   *Nval = 1 + (int) Nv / 8; }
     22    break;
     23  case  'L':
     24    { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
     25    break;
     26  case  'A':
     27    { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
     28    break;
     29  case  'B':
     30    { *Nbytes = 1;  strcpy (type, "char");   *Nval = Nv;               }
     31    break;
     32  case  'I':
     33    { *Nbytes = 2;  strcpy (type, "short");  *Nval = Nv;               }
     34    break;
     35  case  'J':
     36    { *Nbytes = 4;  strcpy (type, "int");    *Nval = Nv;               }
     37    break;
     38  case  'K':
     39    { *Nbytes = 8;  strcpy (type, "int64_t"); *Nval = Nv;              }
     40    break;
     41  case  'E':
     42    { *Nbytes = 4;  strcpy (type, "float");  *Nval = Nv;               }
     43    break;
     44  case  'D':
     45    { *Nbytes = 8;  strcpy (type, "double"); *Nval = Nv;               }
     46    break;
     47  case  'P':
     48    { *Nbytes = 8;  strcpy (type, "var");    *Nval = 2*Nv;             }
     49    break;
     50  case  'C':
     51    { *Nbytes = 8;  strcpy (type, "float");  *Nval = 2*Nv;             }
     52    break;
     53  case  'M':
     54    { *Nbytes = 16; strcpy (type, "double"); *Nval = 2*Nv;             }
     55    break;
     56  default:
     57    return (FALSE);
     58  }
    3159
    3260  return (TRUE);
    … …  
    3967   I - 16 bit int
    4068   J - 32 bit int
     69   K - 64 bit int
    4170   A - char
    4271   E - float
    … …  
    199228  short *tmpShort;
    200229  int *tmpInt;
     230  int64_t *tmpInt64;
    201231
    202232  off = 0;
    … …  
    229259        continue;
    230260    }
     261    if ((tscale == 1.0) && (tzero == 0.0)) {
     262        off += Nchar;
     263        continue;
     264    }
    231265
    232266    if (!strcmp (type, "char"))   {
    … …  
    252286          *tmpInt = *tmpInt * tscale + tzero;
    253287        }
     288      }
     289    }
     290    if (!strcmp (type, "int64_t"))   {
     291      for (j = 0; j < Ny; j++) {
     292        for (n = 0; n < Nval; n++) {
     293          tmpInt64 = (int64_t *)&ftable[0].buffer[j*Nx + n*Nbytes + off];
     294          *tmpInt64 = *tmpInt64 * tscale + tzero;
     295        }
    254296      }
    255297    }
    … …  
    269311  short *tmpShort;
    270312  int *tmpInt;
     313  int64_t *tmpInt64;
    271314
    272315  off = 0;
    … …  
    299342        continue;
    300343    }
     344    if ((tscale == 1.0) && (tzero == 0.0)) {
     345        off += Nchar;
     346        continue;
     347    }
    301348
    302349    if (!strcmp (type, "char"))   {
    … …  
    324371      }
    325372    }
     373    if (!strcmp (type, "int64_t"))   {
     374      for (j = 0; j < Ny; j++) {
     375        for (n = 0; n < Nval; n++) {
     376          tmpInt64 = (int64_t *)&ftable[0].buffer[j*Nx + n*Nbytes + off];
     377          *tmpInt64 = *tmpInt64 * tscale + tzero;
     378        }
     379      }
     380    }
    326381    off += Nchar;
    327382  }

  • trunk/Ohana/src/tools/src/ftable.c

    r16387 r26384  
    361361            memcpy (line, &data[i*Nv*Nb + Nb*j], Nb);
    362362            fprintf (stdout, "%d ", *(int *)line);
     363          }
     364          if (!strcmp (type, "int64_t")) {
     365            memcpy (line, &data[i*Nv*Nb + Nb*j], Nb);
     366            fprintf (stdout, "%d ", *(int64_t*)line);
    363367          }
    364368          if (!strcmp (type, "float")) {
Note: See TracChangeset for help on using the changeset viewer.