IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33653 for trunk


Ignore:
Timestamp:
Apr 1, 2012, 3:03:26 PM (14 years ago)
Author:
eugene
Message:

add tests for PS1_V3; add tests for flat-field correction & ubercal; add tests for parallel dvo; fix up WISE ingest; fix up superCOSMOS ingest; improvements to mkcmf for more detailed testing

Location:
trunk/Ohana/src/addstar
Files:
18 edited
14 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/Makefile

    r31496 r33653  
    2424sedstar     : $(BIN)/sedstar.$(ARCH)
    2525load2mass   : $(BIN)/load2mass.$(ARCH)
     26loadwise    : $(BIN)/loadwise.$(ARCH)
     27loadsupercos: $(BIN)/loadsupercos.$(ARCH)
    2628gztest      : $(BIN)/gztest.$(ARCH)
    2729mkcmf       : $(BIN)/mkcmf.$(ARCH)
    2830
    29 all: addstar sedstar load2mass skycells mkcmf
     31all: addstar sedstar load2mass skycells mkcmf loadwise loadsupercos
     32
     33INSTALL = addstar sedstar load2mass skycells mkcmf loadwise loadsupercos
    3034
    3135# I need to fix the client/server version of addstar now that I have dropped Stars
    … …  
    214218$(SRC)/psps_ids.$(ARCH).o
    215219
     220LOAD-WISE = \
     221$(SRC)/loadwise.$(ARCH).o \
     222$(SRC)/loadwise_prelim_rawdata.$(ARCH).o \
     223$(SRC)/loadwise_ops.$(ARCH).o \
     224$(SRC)/loadwise_star_full.$(ARCH).o \
     225$(SRC)/args_loadwise.$(ARCH).o \
     226$(SRC)/ConfigInit.$(ARCH).o \
     227$(SRC)/Shutdown.$(ARCH).o \
     228$(SRC)/StarOps.$(ARCH).o \
     229$(SRC)/find_matches_refstars.$(ARCH).o \
     230$(SRC)/build_links.$(ARCH).o \
     231$(SRC)/replace_match.$(ARCH).o \
     232$(SRC)/SetSignals.$(ARCH).o \
     233$(SRC)/update_coords.$(ARCH).o \
     234$(SRC)/SkyRegionUtils.$(ARCH).o \
     235$(SRC)/psps_ids.$(ARCH).o
     236
     237LOAD-SUPERCOS = \
     238$(SRC)/loadsupercos.$(ARCH).o \
     239$(SRC)/loadsupercos_ops.$(ARCH).o \
     240$(SRC)/loadsupercos_survey.$(ARCH).o \
     241$(SRC)/loadsupercos_plates.$(ARCH).o \
     242$(SRC)/loadsupercos_rawdata.$(ARCH).o \
     243$(SRC)/args_loadsupercos.$(ARCH).o \
     244$(SRC)/ConfigInit.$(ARCH).o \
     245$(SRC)/Shutdown.$(ARCH).o \
     246$(SRC)/StarOps.$(ARCH).o \
     247$(SRC)/find_matches_refstars.$(ARCH).o \
     248$(SRC)/build_links.$(ARCH).o \
     249$(SRC)/replace_match.$(ARCH).o \
     250$(SRC)/SetSignals.$(ARCH).o \
     251$(SRC)/update_coords.$(ARCH).o \
     252$(SRC)/SkyRegionUtils.$(ARCH).o \
     253$(SRC)/parse_time.$(ARCH).o \
     254$(SRC)/psps_ids.$(ARCH).o
     255
    216256SEDSTAR = \
    217257$(SRC)/sedstar.$(ARCH).o \
    … …  
    241281$(SKYCELLS)   : $(INC)/addstar.h
    242282$(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h
     283$(LOAD-WISE)  : $(INC)/addstar.h $(INC)/WISE.h
    243284$(MKCMF)      : $(INC)/mkcmf.h
    244285
    … …  
    250291$(BIN)/sedstar.$(ARCH)     : $(SEDSTAR)
    251292$(BIN)/load2mass.$(ARCH)   : $(LOAD-2MASS)
     293$(BIN)/loadwise.$(ARCH)    : $(LOAD-WISE)
     294$(BIN)/loadsupercos.$(ARCH): $(LOAD-SUPERCOS)
    252295$(BIN)/skycells.$(ARCH)    : $(SKYCELLS)
    253296$(BIN)/mkcmf.$(ARCH)       : $(MKCMF)
    254 
    255 INSTALL = addstar sedstar load2mass skycells mkcmf
    256297
    257298# I need to fix the client/server version of addstar now that I have dropped Stars

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

    r32346 r33653  
    2121short WISE_W1, WISE_W2, WISE_W3, WISE_W4;
    2222
    23 SkyTable *get2mass_acc (SkyRegion *patch, char *path, char *accel);
    24 Stars    *get2mass_2DR_data (SkyRegion *region, char *filename, SkyRegion *patch, int photcode, int *nstars);
    25 Stars    *get2mass_AS_data (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
    26 Stars    *get2mass_AS_rawdata (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars);
     23AddstarClientOptions args_loadwise (int argc, char **argv, AddstarClientOptions options);
    2724
    28 SkyTable *scan2mass_acc (char *path, char *accel);
    29 int       scan2mass_as_data (char *filename);
     25int loadwise_prelim_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options);
    3026
    31 char     *skipNbounds (char *line, char bound, int Nbound, int Nbyte);
    32 e_time    get2mass_time (char *ptr, int Nbound, int Nbyte);
    33 e_time    get2mass_date (char *ptr, int Nbound, int Nmax);
     27int getWISE_setup ();
     28int getWISE_coords (char *line, double *R, double *D, int Nmax);
     29e_time getWISE_date (char *ptr, int Nbound, int Nmax);
     30e_time getWISE_time (char *ptr, int Nbound, int Nmax);
    3431
    35 int       load2mass_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options);
    36 SkyTable *load2mass_acc (char *path, char *accel);
    37 int       get2mass_3star (Stars *star, char *line, int Nmax);
    38 int       load2mass_catalog (Catalog *catalog, Stars *stars, int Nstars);
     32int getWISE_sortStars (WISE_Stars *tstars, int Ntstars);
     33int loadwise_star_full (Stars **star, char *line, int Nmax);
    3934
    40 int       get2mass_setup (int photcode);
    41 int       get2mass_coords (char *line, double *R, double *D, int Nmax);
    42 int       get2mass_star (Stars *star, char *line, int Nmax);
    43 int       get2mass_3star (Stars *star, char *line, int Nmax);
     35char *nextWISEfield (char *line);
     36int setWISE_cc_flag (Stars *star, char qual);
     37int setWISE_ph_qual (Stars *star, char qual);
     38int setWISE_rd_flag (Stars *star, char qual);
     39int setWISE_bl_flag (Stars *star, char qual);
     40int setWISE_gal_flag (Stars *star, char qual);
     41int setWISE_mp_flag (Stars *star, char qual);
     42int setWISE_dup_flag (Stars *star, char qual);
     43int setWISE_use_flag (Stars *star, char qual);
    4444
    45 int get2mass_3star_full (Stars **star, char *line, int Nmax);
    46 char *next2MASSfield (char *line);
    47 int set2MASS_ph_qual (Stars *star, char qual);
    48 int set2MASS_rd_flag (Stars *star, char qual);
    49 int set2MASS_cc_flag (Stars *star, char qual);
    50 int set2MASS_bl_flag (Stars *star, char qual);
    51 int set2MASS_gal_flag (Stars *star, char qual);
    52 int set2MASS_mp_flag (Stars *star, char qual);
    53 int set2MASS_dup_flag (Stars *star, char qual);
    54 int set2MASS_use_flag (Stars *star, char qual);
    55 int get2mass_sortStars (WISE_Stars *tstars, int Ntstars);
     45char *skipNbounds (char *line, char bound, int Nbound, int Nbyte);

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

    r31669 r33653  
    2828    fprintf (stderr, "ERROR:  photcode %s not found in photcode table\n", NAME); \
    2929    exit (0); }
     30
     31# define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
     32# define dSIN(A)   ((double) sin ((double)RAD_DEG*A))
    3033
    3134typedef struct {

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

    r28967 r33653  
    213213
    214214  /* get detection filtering mask */
    215   if(!ScanConfig (config, "DETECTIONFILTER", "%d", 0, &options.detectionFilter)) {
     215  if (!ScanConfig (config, "DETECTIONFILTER", "%d", 0, &options.detectionFilter)) {
    216216      options.detectionFilter = 0;
    217       fprintf (stderr, "Could not find 'DETECTIONFILTER' in config, using '%x'\n", options.detectionFilter);
    218217  }
    219218

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

    r30613 r33653  
    44
    55# define NFILTER 5
    6 # define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
    7 # define dSIN(A)   ((double) sin ((double)RAD_DEG*A))
    86
    97# define GET_COLUMN_5(NAME,TYPE) \
    … …  
    279277}
    280278
    281 // ha/dec -> alt/az
    282 int altaz (double *alt, double *az, double ha, double dec, double latitude) {
    283 
    284   double sind, sinh, cosh;
    285 
    286   sind = dSIN (dec) * dSIN (latitude) + dCOS (dec) * dCOS (ha) * dCOS (latitude);
    287   *alt  = DEG_RAD * asin (sind);
    288 
    289   sinh = - dCOS (dec) * dSIN (ha);
    290   cosh =   dSIN (dec) * dCOS (latitude) - dCOS (dec) * dCOS (ha) * dSIN (latitude);
    291 
    292   *az = DEG_RAD * atan2 (sinh, cosh);
    293 
    294   // I may need to use the parallactic angle to get the right plate rotation...
    295   // sinh = -dCOS(az) * dSIN(alt) * dSIN(ha) * dSIN(latitude) + dSIN(az) * dSIN(alt) * dCOS(ha) - dSIN(ha) * dCOS(alt) * dCOS(latitude);
    296   // cosh = -dSIN(az) * dSIN(ha) * dSIN(latitude) - dCOS(az) * dCOS(ha);
    297   // rot = -DEG_RAD * atan2 (sinh, cosh);
    298   return TRUE;
    299 }
    300 
    301279int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) {
    302280

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

    r28840 r33653  
    291291
    292292    for (j = 0; j < Nsecfilt; j++) {
    293       catalog[0].secfilt[Nave*Nsecfilt+j].M     = NAN;
    294       catalog[0].secfilt[Nave*Nsecfilt+j].dM    = NAN;
    295       catalog[0].secfilt[Nave*Nsecfilt+j].Xm    = NAN_S_SHORT;
    296       catalog[0].secfilt[Nave*Nsecfilt+j].M_20  = NAN_S_SHORT;
    297       catalog[0].secfilt[Nave*Nsecfilt+j].M_80  = NAN_S_SHORT;
    298       catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
    299       catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
     293      catalog[0].secfilt[Nave*Nsecfilt+j].M           = NAN;
     294      catalog[0].secfilt[Nave*Nsecfilt+j].Map         = NAN;
     295      catalog[0].secfilt[Nave*Nsecfilt+j].dM          = NAN;
     296      catalog[0].secfilt[Nave*Nsecfilt+j].Xm          = NAN_S_SHORT;
     297      catalog[0].secfilt[Nave*Nsecfilt+j].M_20        = NAN_S_SHORT;
     298      catalog[0].secfilt[Nave*Nsecfilt+j].M_80        = NAN_S_SHORT;
     299      catalog[0].secfilt[Nave*Nsecfilt+j].Ncode       = 0;
     300      catalog[0].secfilt[Nave*Nsecfilt+j].Nused       = 0;
     301      catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000;
     302      catalog[0].secfilt[Nave*Nsecfilt+j].flags       = 0;
     303      catalog[0].secfilt[Nave*Nsecfilt+j].dummy[0]    = 0;
     304      catalog[0].secfilt[Nave*Nsecfilt+j].dummy[1]    = 0;
    300305    }
    301306

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

    r32632 r33653  
    292292
    293293    for (j = 0; j < Nsecfilt; j++) {
    294       catalog[0].secfilt[Nave*Nsecfilt+j].M  = NAN;
    295       catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN;
    296       catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT;
    297       catalog[0].secfilt[Nave*Nsecfilt+j].M_20  = NAN_S_SHORT;
    298       catalog[0].secfilt[Nave*Nsecfilt+j].M_80  = NAN_S_SHORT;
    299       catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
    300       catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
     294      catalog[0].secfilt[Nave*Nsecfilt+j].M           = NAN;
     295      catalog[0].secfilt[Nave*Nsecfilt+j].Map         = NAN;
     296      catalog[0].secfilt[Nave*Nsecfilt+j].dM          = NAN;
     297      catalog[0].secfilt[Nave*Nsecfilt+j].Xm          = NAN_S_SHORT;
     298      catalog[0].secfilt[Nave*Nsecfilt+j].M_20        = NAN_S_SHORT;
     299      catalog[0].secfilt[Nave*Nsecfilt+j].M_80        = NAN_S_SHORT;
     300      catalog[0].secfilt[Nave*Nsecfilt+j].Ncode       = 0;
     301      catalog[0].secfilt[Nave*Nsecfilt+j].Nused       = 0;
     302      catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000;
     303      catalog[0].secfilt[Nave*Nsecfilt+j].flags       = 0;
     304      catalog[0].secfilt[Nave*Nsecfilt+j].dummy[0]    = 0;
     305      catalog[0].secfilt[Nave*Nsecfilt+j].dummy[1]    = 0;
    301306    }
    302307

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

    r29001 r33653  
    1313  int Nsecfilt;
    1414
    15   if ((NREFSTAR_GROUP != 1) && (NREFSTAR_GROUP != 3)) {
     15  switch (NREFSTAR_GROUP) {
     16    case 1:
     17    case 3:
     18    case 4:
     19      break;
     20    default:
    1621      fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n");
    1722      exit (1);

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

    r32346 r33653  
    22# include "WISE.h"
    33
     4/* This is the DVO program to upload WISE detections into a DVO database.
     5   It is modeled on the load2mass program, but unlike that case, it is expected to be run only
     6   rarely (once?).  The WISE data are delivered as *.bz2 files.  This program assumes the files
     7   have been decompressed.  It does not allow a subset of the sky to be uploaded; entire WISE
     8   files are loaded if supplied on the command line.
     9
     10   USAGE: loadwise -D CATDIR (catdir) (wisefile) [...more files]
     11
     12 */
     13
    414int main (int argc, char **argv) {
    515
    6   char *path;
    716  int i;
    8   SkyTable *sky, *skyWISE;
     17  SkyTable *sky;
    918  SkyList *skylist = NULL;
    10   SkyList *overlap = NULL;
    1119  AddstarClientOptions options;
    1220
    1321  // need to construct these options with args_loadWISE...
    1422  options = ConfigInit (&argc, argv);
    15   options = args_loadWISE (argc, argv, options);
     23  options = args_loadwise (argc, argv, options);
    1624
    1725  // load the full sky description table:
    … …  
    3038  }
    3139
    32   path = WISE_DIR_PRELIM;
    33 
    34   // the accel.dat file has the raw filenames
    35   // test if the file exists, or else try the .gz version
    36   // XXX can we build an accelerator?
    37   skyWISE = loadWISE_acc (path, "accel.dat");
    38  
    39   for (i = 0; i < skyWISE[0].Nregions; i++) {
    40     // check if any of the skylist entries overlap this WISE catalog:
    41     overlap = SkyListByBounds_List (skylist, -1, skyWISE[0].regions[i].Rmin, skyWISE[0].regions[i].Rmax, skyWISE[0].regions[i].Dmin, skyWISE[0].regions[i].Dmax);
    42     if (overlap[0].Nregions == 0) {
    43       SkyListFree (overlap);
    44       continue;
    45     }
    46    
    47     fprintf (stderr, "loading %s\n", skyWISE[0].filename[i]);
    48     loadWISE_as_rawdata (overlap, skyWISE[0].filename[i], options);
    49     SkyListFree (overlap);
     40  for (i = 1; i < argc; i++) {
     41      fprintf (stderr, "loading %s\n", argv[i]);
     42      loadwise_prelim_rawdata (skylist, argv[i], options);
    5043  }
    5144  exit (0);
    5245} 
     46

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

    r32346 r33653  
    11# include "addstar.h"
    22# include "WISE.h"
    3 
    4 static int FilterSkip;
    5 static int TimeSkip;
    6 static int Qentry;
    7 static int Photcode;
    83
    94int getWISE_setup () {

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

    r32346 r33653  
    1010# define DEBUG 0
    1111
    12 int loadWISE_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) {
     12int loadwise_prelim_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) {
    1313 
    1414  int i, j, verbose;
    … …  
    159159        stars[Nstars+3][0].average.R = tstars[j].R;
    160160        stars[Nstars+3][0].average.D = tstars[j].D;
    161         getWISE_star_full (&stars[Nstars], &buffer[offset], Nbyte - offset);
     161        loadwise_star_full (&stars[Nstars], &buffer[offset], Nbyte - offset);
    162162
    163163        tstars[j].flag = TRUE;

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

    r32346 r33653  
    22# include "WISE.h"
    33
    4 // XXX check to see if desired output format is PS1_V1 or later?  (use 16bit version if not?)
     4int setWISE_ph_qual (Stars *star, char qual);
     5int setWISE_var_flag (Stars *star, char qual);
     6int setWISE_ext_flag (Stars **star, char *ptr);
     7int setWISE_cc_flag (Stars *star, char qual);
     8int setWISE_sat_flag (Stars *star, char *ptr);
     9int setWISE_blend_flag (Stars **star, char *ptr);
    510
    611// fill in the data for a WISE quad star.  takes a pointer to the start of the line the
    712// RA and DEC have already been set
    8 int getWISE_star_full (Stars **star, char *line, int Nmax) {
    9 
     13int loadwise_star_full (Stars **star, char *line, int Nmax) {
     14
     15  int i;
    1016  char *ptr;
    11   double dMfull;
    12   double jd;
    1317
    1418  if (line == NULL) Shutdown ("format error in WISE");
    1519
     20  ptr = line;
     21
     22  // I can assign dRA and dDEC to dX and dY if I can consistently set posangle and pltscale
     23  star[0][0].measure.posangle = 0.0;
     24  star[0][0].measure.pltscale = 1.0;
     25
    1626  ptr = skipNbounds (ptr, '|', 3, Nmax); // skip: desig, ra, dec,
    17   star[0][0].measure.FWx = ToShortPixels(strtod (ptr, NULL));
    18   ptr = nextWISEfield (ptr); // sig_dec
    19   star[0][0].measure.FWy = ToShortPixels(strtod (ptr, NULL));
    20   ptr = nextWISEfield (ptr); // sig_radec
    21   star[0][0].measure.theta = ToShortPixels(strtod (ptr, NULL));
    22 
    23   star[1][0].measure.FWx   = star[0][0].measure.FWx;
    24   star[1][0].measure.FWy   = star[0][0].measure.FWy;
    25   star[1][0].measure.theta = star[0][0].measure.theta;
    26   star[2][0].measure.FWx   = star[0][0].measure.FWx;
    27   star[2][0].measure.FWy   = star[0][0].measure.FWy;
    28   star[2][0].measure.theta = star[0][0].measure.theta;
    29   star[3][0].measure.FWx   = star[0][0].measure.FWx;
    30   star[3][0].measure.FWy   = star[0][0].measure.FWy;
    31   star[3][0].measure.theta = star[0][0].measure.theta;
    32 
    33   ptr = skipNbounds (ptr, '|', 10, Nmax); // skip: glon, glat, elon, elat, wx, wy, cntr, source_id, coadd_id, src
     27  star[0][0].measure.dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra
     28  ptr = nextWISEfield (ptr);
     29  star[0][0].measure.dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec
     30  ptr = nextWISEfield (ptr);
     31
     32  // we only know a single set of values for all 4 bands
     33  for (i = 1; i < 4; i++) {
     34      star[i][0].measure.dXccd    = star[0][0].measure.dXccd;
     35      star[i][0].measure.dYccd    = star[0][0].measure.dYccd;
     36      star[i][0].measure.posangle = star[0][0].measure.posangle;
     37      star[i][0].measure.pltscale = star[0][0].measure.pltscale;
     38  }
     39  ptr = skipNbounds (ptr, '|', 5, Nmax); // skip: sig_radec, glon, glat, elon, elat
     40
     41  star[0][0].measure.Xccd = strtod (ptr, NULL); // wx
     42  ptr = nextWISEfield (ptr); // skip wx
     43  star[0][0].measure.Yccd = strtod (ptr, NULL); // wy
     44  ptr = nextWISEfield (ptr); // skip wy
     45
     46  for (i = 1; i < 4; i++) {
     47      star[i][0].measure.Xccd    = star[0][0].measure.Xccd;
     48      star[i][0].measure.Yccd    = star[0][0].measure.Yccd;
     49  }
     50  ptr = skipNbounds (ptr, '|', 4, Nmax); // skip: cntr, source_id, coadd_id, src
    3451
    3552  // W1
    36   star[0][0].measure.M  = strtod (ptr, NULL);
    37   ptr = nextWISEfield (ptr); // j_cmsig
    38   star[0][0].measure.dM = strtod (ptr, NULL);
    39   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w1snr, w1rchi2
    40 
    41   // W2
    42   star[1][0].measure.M  = strtod (ptr, NULL);
    43   ptr = nextWISEfield (ptr); // j_cmsig
    44   star[1][0].measure.dM = strtod (ptr, NULL);
    45   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w2snr, w2rchi2
    46 
    47   // W3
    48   star[2][0].measure.M  = strtod (ptr, NULL);
    49   ptr = nextWISEfield (ptr); // j_cmsig
    50   star[2][0].measure.dM = strtod (ptr, NULL);
    51   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w3snr, w3rchi2
    52 
    53   // W4
    54   star[3][0].measure.M  = strtod (ptr, NULL);
    55   ptr = nextWISEfield (ptr); // j_cmsig
    56   star[3][0].measure.dM = strtod (ptr, NULL);
    57   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w4snr, w4rchi2
    58 
    59   star[0][0].measure.photFlags = 0;
    60   star[1][0].measure.photFlags = 0;
    61   star[2][0].measure.photFlags = 0;
    62   star[3][0].measure.photFlags = 0;
    63   ptr = skipNbounds (ptr, '|', 8, Nmax); // skip fields
    64 
    65   setWISE_cc_flag (star[0], ptr[0]); // cc_flg
    66   setWISE_cc_flag (star[1], ptr[1]); // cc_flg
    67   setWISE_cc_flag (star[2], ptr[2]); // cc_flg
    68   setWISE_cc_flag (star[3], ptr[2]); // cc_flg
    69   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: ext_flg, var_flg
    70 
    71   setWISE_ph_qual (star[0], ptr[0]); // ph_qual
    72   setWISE_ph_qual (star[1], ptr[1]); // ph_qual
    73   setWISE_ph_qual (star[2], ptr[2]); // ph_qual
    74   setWISE_ph_qual (star[3], ptr[3]); // ph_qual
    75   ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: ext_flg, var_flg
    76 
    77   // setWISE_rd_flag (star[0], ptr[0]);
    78   // setWISE_bl_flag (star[0], ptr[0]);
    79   // setWISE_gal_flag (star[0], ptr[0]);
    80   // setWISE_mp_flag (star[0], ptr[0]);
    81 
    82   // ptr = nextWISEfield (ptr); // x_scan
    83   // star[0][0].measure.Xccd = strtod (ptr, NULL);
    84   // star[2][0].measure.Xccd = star[1][0].measure.Xccd = star[0][0].measure.Xccd;
    85   // ptr = nextWISEfield (ptr); // jdate (julian date)
    86 
    87   jd = 2455263.0; // NOTE : WISE prelim release does not contain per-detection time info.
     53  for (i = 0; i < 4; i++) {
     54      star[0][0].measure.M  = strtod (ptr, NULL); // w?mpro
     55      ptr = nextWISEfield (ptr);
     56      star[0][0].measure.dM = strtod (ptr, NULL); // w?sigmpro
     57      ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w?sigmpro, w?snr
     58      star[0][0].measure.psfChisq = strtod (ptr, NULL); // w?rchi2
     59      ptr = nextWISEfield (ptr); // skip : w1rchi2
     60
     61      // init the photFlags field
     62      star[i][0].measure.photFlags = 0;
     63  }
     64
     65  ptr = nextWISEfield (ptr); // skip: rchi2
     66
     67  // set blend flags for all 4 measures
     68  setWISE_blend_flag (star, ptr); // nb & na both used here
     69  ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: nb, na
     70
     71  for (i = 0; i < 4; i++) {
     72      setWISE_sat_flag (star[i], ptr); // w1sat
     73      ptr = nextWISEfield (ptr);
     74  }
     75  ptr = nextWISEfield (ptr); // skip satnum
     76
     77  for (i = 0; i < 4; i++) {
     78      setWISE_cc_flag (star[i], ptr[i]); // cc_flg
     79  }
     80  ptr = nextWISEfield (ptr); // skip cc_flags
     81
     82  // set ext flags for all 4 measures
     83  setWISE_ext_flag (star, ptr); // ext_flg
     84  ptr = nextWISEfield (ptr); // skip ext_flags
     85
     86  for (i = 0; i < 4; i++) {
     87      setWISE_var_flag (star[i], ptr[i]); // var_flg
     88  }
     89  ptr = nextWISEfield (ptr); // skip var_flags
     90
     91  for (i = 0; i < 4; i++) {
     92      setWISE_ph_qual (star[i], ptr[i]); // ph_qual
     93  }
     94
     95  double jd = 2455263.0; // NOTE : WISE prelim release does not contain per-detection time info.
    8896  // the release is based on data taken in the period 14 January 2010 to 29 April 2010
    8997  star[0][0].measure.t = ohana_jd_to_sec (jd);
    … …  
    92100  star[3][0].measure.t = star[0][0].measure.t;
    93101
    94   // ptr = nextWISEfield (ptr); // j_psfchi
    95   // star[0][0].measure.psfChisq = strtod (ptr, NULL);
    96   // ptr = nextWISEfield (ptr); // h_psfchi
    97   // star[1][0].measure.psfChisq = strtod (ptr, NULL);
    98   // ptr = nextWISEfield (ptr); // k_psfchi
    99   // star[2][0].measure.psfChisq = strtod (ptr, NULL);
    100   //
    101   // ptr = nextWISEfield (ptr); // j_m_stdap
    102   // star[0][0].measure.Map = strtod (ptr, NULL);
    103   // ptr = nextWISEfield (ptr); // j_msig_stdap (skip?)
    104   //
    105   // ptr = nextWISEfield (ptr); // h_m_stdap
    106   // star[1][0].measure.Map = strtod (ptr, NULL);
    107   // ptr = nextWISEfield (ptr); // h_msig_stdap (skip?)
    108   //
    109   // ptr = nextWISEfield (ptr); // k_m_stdap
    110   // star[2][0].measure.Map = strtod (ptr, NULL);
    111   // ptr = nextWISEfield (ptr); // k_msig_stdap (skip?)
    112   //
    113   // ptr = nextWISEfield (ptr); // dist_edge_ns (skip)
    114   // ptr = nextWISEfield (ptr); // dist_edge_ew (skip)
    115   // ptr = nextWISEfield (ptr); // dist_edge_flg (skip)
    116   //
    117   // ptr = nextWISEfield (ptr); // dup_src (one flag for all filters)
    118   // setWISE_dup_flag (star[0], ptr[0]);
    119   // setWISE_dup_flag (star[1], ptr[0]);
    120   // setWISE_dup_flag (star[2], ptr[0]);
    121   //
    122   // ptr = nextWISEfield (ptr); // use_src (one flag for all filters)
    123   // setWISE_use_flag (star[0], ptr[0]);
    124   // setWISE_use_flag (star[1], ptr[0]);
    125   // setWISE_use_flag (star[2], ptr[0]);
    126 
    127   star[0][0].measure.photcode  = TM_J;
     102  star[0][0].measure.photcode  = WISE_W1;
    128103  star[0][0].measure.detID   = 0;
    129104  star[0][0].measure.imageID = 0;
    130105
    131   star[1][0].measure.photcode  = TM_H;
     106  star[1][0].measure.photcode  = WISE_W2;
    132107  star[1][0].measure.detID   = 0;
    133108  star[1][0].measure.imageID = 0;
    134109
    135   star[2][0].measure.photcode  = TM_K;
     110  star[2][0].measure.photcode  = WISE_W3;
    136111  star[2][0].measure.detID   = 0;
    137112  star[2][0].measure.imageID = 0;
     113
     114  star[3][0].measure.photcode  = WISE_W4;
     115  star[3][0].measure.detID   = 0;
     116  star[3][0].measure.imageID = 0;
    138117
    139118  return TRUE;
    … …  
    153132}
    154133
     134# define FLAG_PH_A            0x00000001 // quality flag 'A'
     135# define FLAG_PH_B            0x00000002 // quality flag 'B'
     136# define FLAG_PH_C            0x00000004 // quality flag 'C'
     137# define FLAG_PH_U            0x00000008 // quality flag 'U'
     138# define FLAG_PH_X            0x00000010 // quality flag 'X'
     139
     140# define FLAG_SATURATED_PIX   0x00000100 // sat > 0.0
     141
     142# define FLAG_CC_PERSIST      0x00010000 // 'p' or 'P'
     143# define FLAG_CC_HALO         0x00020000 // 'h' or 'H'
     144# define FLAG_CC_GHOST        0x00020000 // 'o' or 'O' (note overloading)
     145# define FLAG_CC_SPIKE        0x00040000 // 'd' or 'D'
     146# define FLAG_CC_SPURIOUS     0x00080000 // Upper Case letters
     147
     148# define FLAG_BLEND_ACTIVE    0x00100000 // nb > 1, na == 0
     149# define FLAG_BLEND_PASSIVE   0x00200000 // nb > 1, na == 1
     150
     151# define FLAG_EXTENDED        0x01000000 // ext == 1
     152
     153# define FLAG_VARIABLE_LEVEL1 0x10000000 // var_flg == 1 to 4
     154# define FLAG_VARIABLE_LEVEL2 0x20000000 // var_flg == 5 to 7
     155# define FLAG_VARIABLE_LEVEL3 0x40000000 // var_flg == 8 or 9
     156
     157int setWISE_blend_flag (Stars **star, char *ptr) {
     158
     159    int nb = atoi (ptr);
     160    if (nb == 1) return TRUE;
     161
     162    ptr = nextWISEfield (ptr); // skip to na
     163    int na = atoi (ptr);
     164
     165    if (na == 0) {
     166        star[0][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
     167        star[1][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
     168        star[2][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
     169        star[3][0].measure.photFlags |= FLAG_BLEND_ACTIVE;
     170    } else {
     171        star[0][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
     172        star[1][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
     173        star[2][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
     174        star[3][0].measure.photFlags |= FLAG_BLEND_PASSIVE;
     175    }
     176    return TRUE;
     177}
     178
     179int setWISE_sat_flag (Stars *star, char *ptr) {
     180
     181    float sat = strtod (ptr, NULL);
     182    if (sat > 0.0) {
     183        star[0].measure.photFlags |= FLAG_SATURATED_PIX;
     184    }
     185    return TRUE;
     186}
     187
    155188int setWISE_cc_flag (Stars *star, char qual) {
    156189
    157190  switch (qual) {
    158     case 'p': star[0].measure.photFlags |= 0x00010000; break;
    159     case 'P': star[0].measure.photFlags |= 0x00010000; break;
    160     case 'h': star[0].measure.photFlags |= 0x00020000; break;
    161     case 'H': star[0].measure.photFlags |= 0x00020000; break;
    162     case 'd': star[0].measure.photFlags |= 0x00040000; break;
    163     case 'D': star[0].measure.photFlags |= 0x00040000; break;
    164     case 'o': star[0].measure.photFlags |= 0x00020000; break;
    165     case 'O': star[0].measure.photFlags |= 0x00020000; break;
     191    case 'p': star[0].measure.photFlags |= FLAG_CC_PERSIST; break;
     192    case 'h': star[0].measure.photFlags |= FLAG_CC_HALO;    break;
     193    case 'd': star[0].measure.photFlags |= FLAG_CC_SPIKE;   break;
     194    case 'o': star[0].measure.photFlags |= FLAG_CC_GHOST;   break;
     195
     196    case 'P': star[0].measure.photFlags |= FLAG_CC_PERSIST | FLAG_CC_SPURIOUS; break;
     197    case 'H': star[0].measure.photFlags |= FLAG_CC_HALO    | FLAG_CC_SPURIOUS; break;
     198    case 'D': star[0].measure.photFlags |= FLAG_CC_SPIKE   | FLAG_CC_SPURIOUS; break;
     199    case 'O': star[0].measure.photFlags |= FLAG_CC_GHOST   | FLAG_CC_SPURIOUS; break;
     200
    166201    case '0': break;
    167202    default:
    … …  
    172207}
    173208
     209int setWISE_ext_flag (Stars **star, char *ptr) {
     210
     211    int ext = atoi (ptr);
     212    if (ext == 1) {
     213        star[0][0].measure.photFlags |= FLAG_EXTENDED;
     214        star[1][0].measure.photFlags |= FLAG_EXTENDED;
     215        star[2][0].measure.photFlags |= FLAG_EXTENDED;
     216        star[3][0].measure.photFlags |= FLAG_EXTENDED;
     217    }
     218    return TRUE;
     219}
     220
     221int setWISE_var_flag (Stars *star, char qual) {
     222
     223    int value = qual - '0';
     224    assert (value >= 0);
     225    assert (value <= 9);
     226    if (value == 0) return TRUE;
     227    if (value < 5) {
     228        star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL1;
     229        return TRUE;
     230    }
     231    if (value < 8) {
     232        star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL2;
     233        return TRUE;
     234    }
     235    star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL3;
     236    return TRUE;
     237}
     238
    174239int setWISE_ph_qual (Stars *star, char qual) {
    175240
    176241  switch (qual) {
    177     case 'A': star[0].measure.photFlags |= 0x00000001; break; // was: 0x0004
    178     case 'B': star[0].measure.photFlags |= 0x00000002; break; // was: 0x0005
    179     case 'C': star[0].measure.photFlags |= 0x00000004; break; // was: 0x0006
    180     case 'U': star[0].measure.photFlags |= 0x00000040; break; // was: 0x0001
    181     case 'X': star[0].measure.photFlags |= 0x00000080; break; // was: 0x0000
     242    case 'A': star[0].measure.photFlags |= FLAG_PH_A; break;
     243    case 'B': star[0].measure.photFlags |= FLAG_PH_B; break;
     244    case 'C': star[0].measure.photFlags |= FLAG_PH_C; break;
     245    case 'U': star[0].measure.photFlags |= FLAG_PH_U; break;
     246    case 'X': star[0].measure.photFlags |= FLAG_PH_X; break;
    182247    default:
    183248      fprintf (stderr, "error!\n");
    … …  
    187252}
    188253
    189 int setWISE_rd_flag (Stars *star, char qual) {
    190 
    191   switch (qual) {
    192     case '0': star[0].measure.photFlags |= 0x00000100; break; // was: 0x0000
    193     case '1': star[0].measure.photFlags |= 0x00000200; break; // was: 0x0010
    194     case '2': star[0].measure.photFlags |= 0x00000400; break; // was: 0x0020
    195     case '3': star[0].measure.photFlags |= 0x00000800; break; // was: 0x0030
    196     case '4': star[0].measure.photFlags |= 0x00001000; break; // was: 0x0040
    197     case '6': star[0].measure.photFlags |= 0x00002000; break; // was: 0x0050
    198     case '9': star[0].measure.photFlags |= 0x00004000; break; // was: 0x0060
    199     default:
    200       fprintf (stderr, "error!\n");
    201       exit (2);
    202   }     
    203   return (TRUE);
    204 }
    205 
    206 int setWISE_bl_flag (Stars *star, char qual) {
    207 
    208   switch (qual) {
    209     case '0': star[0].measure.photFlags &= ~0x00300000; break; // was: ~0x0008
    210     case '1': star[0].measure.photFlags |=  0x00100000; break; // was: ~0x0008
    211     default:  star[0].measure.photFlags |=  0x00200000; break; // was:  0x0008
    212   }     
    213   return (TRUE);
    214 }
    215 
    216 int setWISE_gal_flag (Stars *star, char qual) {
    217 
    218   switch (qual) {
    219     case '0': star[0].measure.photFlags &= ~0x00c00000; break; // was: ~0x0080
    220     case '1': star[0].measure.photFlags |=  0x00400000; break; // was: ~0x0080
    221     default:  star[0].measure.photFlags |=  0x00800000;        // was:  0x0080
    222       star[0].measure.extNsigma = 100.0;
    223       break;
    224   }     
    225   return (TRUE);
    226 }
    227 
    228 int setWISE_mp_flag (Stars *star, char qual) {
    229 
    230   switch (qual) {
    231     case '0': star[0].measure.photFlags &= ~0x03000000; break; // was: ~0x0800
    232     case '1': star[0].measure.photFlags |=  0x01000000; break; // was: ~0x0800
    233     default:  star[0].measure.photFlags |=  0x02000000; break; // was:  0x0800
    234   }     
    235   return (TRUE);
    236 }
    237 
    238 int setWISE_dup_flag (Stars *star, char qual) {
    239 
    240   switch (qual) {
    241     case '0': star[0].measure.photFlags &= ~0x0c000000; break; // was: ~0x1000
    242     case '1': star[0].measure.photFlags |=  0x04000000; break; // was: ~0x1000
    243     default:  star[0].measure.photFlags |=  0x08000000; break; // was:  0x1000
    244   }     
    245   return (TRUE);
    246 }
    247 
    248 int setWISE_use_flag (Stars *star, char qual) {
    249 
    250   switch (qual) {
    251     case '0': star[0].measure.photFlags &= ~0x10000000; break; // was: ~0x2000
    252     case '1': star[0].measure.photFlags |=  0x10000000; break; // was:  0x2000
    253     default:  abort();
    254   }     
    255   return (TRUE);
    256 }
    257 
    258 // unused photFlags:
    259 // 0x0000.8000
    260 // 0x0004.0000
    261 // 0x0008.0000
    262 // 0x2000.0000
    263 // 0x4000.0000
    264 // 0x8000.0000
     254/* unused photFlag bits:
     255
     2560x0000.0020
     2570x0000.0040
     2580x0000.0080
     259
     2600x0000.0200
     2610x0000.0400
     2620x0000.0800
     263
     2640x0000.1000
     2650x0000.2000
     2660x0000.4000
     2670x0000.8000
     268
     2690x0040.0000
     2700x0080.0000
     271
     2720x0200.0000
     2730x0400.0000
     2740x0800.0000
     275
     2760x8000.0000
     277
     278*/

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

    r27582 r33653  
    1414void gauss_init (int Nbin);
    1515double rnd_gauss (double mean, double sigma);
    16 void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
     16void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars);
    1717void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    1818void writeStars_PS1_DEV_1 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    1919void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars);
    2020
     21int ADDNOISE = TRUE;
     22float BAD_PSFQF_FRAC = 0.0;
     23
    2124int main (int argc, char **argv) {
    2225
    … …  
    2528
    2629  int N, Nstars, NSTARS, found;
     30  unsigned int *Flag;
    2731  double *X, *Y, *M, Xmax, Ymax;
    2832
    … …  
    3438  Coords coords;
    3539
     40  static char *photcode = "SIMTEST.r.Chip";
     41  if ((N = get_argument (argc, argv, "-photcode"))) {
     42    remove_argument (N, &argc, argv);
     43    photcode = strcreate (argv[N]);
     44    remove_argument (N, &argc, argv);
     45  }
     46
     47  static double mjd;
     48  mjd = NAN;
     49  if ((N = get_argument (argc, argv, "-mjd"))) {
     50    remove_argument (N, &argc, argv);
     51    mjd = strtod (argv[N], NULL);
     52    remove_argument (N, &argc, argv);
     53  }
     54
    3655  static char *date = "2001-01-01";
    3756  if ((N = get_argument (argc, argv, "-date"))) {
    … …  
    5877  }
    5978
     79  double CRPIX1 = 0.0;
     80  double CRPIX2 = 0.0;
     81  if ((N = get_argument (argc, argv, "-crpix"))) {
     82    remove_argument (N, &argc, argv);
     83    CRPIX1 = atof (argv[N]);
     84    remove_argument (N, &argc, argv);
     85    CRPIX2 = atof (argv[N]);
     86    remove_argument (N, &argc, argv);
     87  }
     88
     89  float exptime = 1.0;
     90  if ((N = get_argument (argc, argv, "-exptime"))) {
     91    remove_argument (N, &argc, argv);
     92    exptime = atof (argv[N]);
     93    remove_argument (N, &argc, argv);
     94  } 
     95
     96  // XXX note that the airmass and ra,dec,mjd can be inconsistent (for a given observatory location)
     97  float airmass = 1.0;
     98  if ((N = get_argument (argc, argv, "-airmass"))) {
     99    remove_argument (N, &argc, argv);
     100    airmass = atof (argv[N]);
     101    remove_argument (N, &argc, argv);
     102  } 
     103
     104  int NX = 0;
     105  int NY = 0;
     106  if ((N = get_argument (argc, argv, "-size"))) {
     107    remove_argument (N, &argc, argv);
     108    NX = atof (argv[N]);
     109    remove_argument (N, &argc, argv);
     110    NY = atof (argv[N]);
     111    remove_argument (N, &argc, argv);
     112  }
     113
    60114  // add support for all cmf types
    61115  static char *type = "PS1_V2";
    … …  
    66120  }
    67121
     122  // add support for all cmf types
     123  int FROM_COORDS = FALSE;
     124  if ((N = get_argument (argc, argv, "-coords"))) {
     125    remove_argument (N, &argc, argv);
     126    FROM_COORDS = TRUE;
     127  }
     128
     129  // expect a field with the photom flags for each detection
     130  int READ_FLAGS = FALSE;
     131  if ((N = get_argument (argc, argv, "-flags"))) {
     132    remove_argument (N, &argc, argv);
     133    READ_FLAGS = TRUE;
     134  }
     135
     136  // add support for all cmf types
     137  ADDNOISE = TRUE;
     138  if ((N = get_argument (argc, argv, "-no-noise"))) {
     139    remove_argument (N, &argc, argv);
     140    ADDNOISE = FALSE;
     141  }
     142
     143  // random bad PSF_QF values
     144  if ((N = get_argument (argc, argv, "-bad-psfqf-frac"))) {
     145    remove_argument (N, &argc, argv);
     146    BAD_PSFQF_FRAC = atof(argv[N]);
     147    remove_argument (N, &argc, argv);
     148  }
     149
    68150  if (argc != 3) {
    69151    fprintf (stderr, "USAGE mkcmf (input) (output) [-date date] [-time time] [-radec ra dec] [-cmftype type]\n");
    … …  
    71153  }
    72154
     155  /* bore site center guess */
     156  strcpy (coords.ctype, "DEC--TAN");
     157  coords.crval1 = RA;
     158  coords.crval2 = DEC;
     159  coords.crpix1 = CRPIX1;
     160  coords.crpix2 = CRPIX2;
     161 
     162  coords.cdelt1 = 0.25/3600.0;
     163  coords.cdelt2 = 0.25/3600.0;
     164
     165  coords.pc1_1  = 1;
     166  coords.pc2_2  = 1;
     167  coords.pc1_2  = 0;
     168  coords.pc2_1  = 0;
     169  coords.Npolyterms = 1;
     170
     171  // load stars and generate complete output fields
    73172  f = fopen (argv[1], "r");
    74173  if (f == NULL) {
    … …  
    83182  ALLOCATE (Y, double, NSTARS);
    84183  ALLOCATE (M, double, NSTARS);
     184  ALLOCATE (Flag, unsigned int, NSTARS);
    85185
    86186  Xmax = Ymax = 0;
    87   while (fscanf (f, "%lf %lf %lf", &X[Nstars], &Y[Nstars], &M[Nstars]) != EOF) {
     187  while (1) {
     188    int status;
     189    double ra, dec, mag, xobs, yobs, xraw, yraw, mraw;
     190    unsigned int flags;
     191    if (FROM_COORDS) {
     192      status = fscanf (f, "%lf %lf %lf %lf %lf %lf", &ra, &dec, &mag, &xraw, &yraw, &mraw);
     193      RD_to_XY (&xobs, &yobs, ra, dec, &coords);
     194    } else {
     195      status = fscanf (f, "%lf %lf %lf", &xobs, &yobs, &mraw);
     196    }
     197    if (status == EOF) break;
     198    if (READ_FLAGS) {
     199      status = fscanf (f, "%x", &flags);
     200      if (status == EOF) {
     201          fprintf (stderr, "error: missing flag for last star?\n");
     202          exit (1);
     203      }
     204      Flag[Nstars] = flags;
     205    } else {
     206      Flag[Nstars] = 0;
     207    }
     208
     209    X[Nstars] = xobs;
     210    Y[Nstars] = yobs;
     211    M[Nstars] = mraw;
     212
    88213    Xmax  = MAX(Xmax, X[Nstars]);
    89214    Ymax  = MAX(Ymax, Y[Nstars]);
    … …  
    94219      REALLOCATE (Y, double, NSTARS);
    95220      REALLOCATE (M, double, NSTARS);
     221      REALLOCATE (Flag, unsigned int, NSTARS);
    96222    }
    97223    Nstars ++;
     224  }
     225  if (NX && NY) {
     226      Xmax = NX;
     227      Ymax = NY;
    98228  }
    99229
    … …  
    106236
    107237  // XXX add minimum needed header fields
    108   gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax + 50));
    109   gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax + 50));
     238  gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax));
     239  gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax));
    110240
    111241  gfits_modify (&header, "NSTARS",   "%d", 1, Nstars);
    112   gfits_modify (&header, "PHOTCODE", "%s", 1, "SIMTEST.r.Chip");
    113   gfits_modify (&header, "DATE-OBS", "%s", 1, date);
    114   gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
     242  gfits_modify (&header, "PHOTCODE", "%s", 1, photcode);
     243  if (isfinite(mjd)) {
     244    gfits_modify (&header, "MJD-OBS", "%lf", 1, mjd);
     245  } else {
     246    gfits_modify (&header, "DATE-OBS", "%s", 1, date);
     247    gfits_modify (&header, "UTC-OBS",  "%s", 1, time);
     248  }
    115249  gfits_modify (&header, "ZERO_PT", "%lf", 1, 25.0);
    116   gfits_modify (&header, "EXPTIME", "%lf", 1, 1.0);
    117   gfits_modify (&header, "AIRMASS", "%lf", 1, 1.0);
     250  gfits_modify (&header, "EXPTIME", "%lf", 1, exptime);
     251  gfits_modify (&header, "AIRMASS", "%lf", 1, airmass);
    118252  gfits_modify (&header, "NASTRO",   "%d", 1, 10);
    119253
    120254  gfits_modify (&header, "IMAGEID",  "%d", 1, (int)(1000*drand48()));
    121255  gfits_modify (&header, "SOURCEID", "%d", 1, (int)(100*drand48()));
    122 
    123   /* bore site center guess */
    124   strcpy (coords.ctype, "DEC--TAN");
    125   coords.crval1 = RA;
    126   coords.crval2 = DEC;
    127   coords.crpix1 = 0.0;
    128   coords.crpix2 = 0.0;
    129  
    130   coords.cdelt1 = 0.25/3600.0;
    131   coords.cdelt2 = 0.25/3600.0;
    132 
    133   coords.pc1_1  = 1;
    134   coords.pc2_2  = 1;
    135   coords.pc1_2  = 0;
    136   coords.pc2_1  = 0;
    137   coords.Npolyterms = 1;
    138256
    139257  PutCoords (&coords, &header);
    … …  
    157275  }
    158276  if (!strcmp(type, "PS1_V2")) {
    159     writeStars_PS1_V2 (&ftable, X, Y, M, Nstars);
     277    writeStars_PS1_V2 (&ftable, X, Y, M, Flag, Nstars);
    160278    found = TRUE;
    161279  }
    … …  
    211329
    212330  A = time(NULL);
    213   for (B = 0; A == time(NULL); B++);
     331  // XXX this is expensive if called a lot (1 sec min)
     332  // for (B = 0; A == time(NULL); B++);
     333  B = A + 10000;
    214334  srand48(B);
    215335 
    … …  
    273393    fSN = 1.0 / sqrt(flux);
    274394
    275     stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0);
    276     stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0);
    277     stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0);
     395    stars[i].X = X[i];
     396    stars[i].Y = Y[i];
     397    stars[i].M = M[i];
     398
     399    if (ADDNOISE) {
     400      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     401      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     402      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
     403    }
    278404
    279405    stars[i].dX = FX * fSN;
    … …  
    312438    fSN = 1.0 / sqrt(flux);
    313439
    314     stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0);
    315     stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0);
    316     stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0);
     440    stars[i].X = X[i];
     441    stars[i].Y = Y[i];
     442    stars[i].M = M[i];
     443
     444    if (ADDNOISE) {
     445      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     446      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     447      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
     448    }
    317449
    318450    stars[i].dX = FX * fSN;
    … …  
    354486    fSN = 1.0 / sqrt(flux);
    355487
    356     stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0);
    357     stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0);
    358     stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0);
     488    stars[i].X = X[i];
     489    stars[i].Y = Y[i];
     490    stars[i].M = M[i];
     491
     492    if (ADDNOISE) {
     493      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     494      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     495      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
     496    }
    359497
    360498    stars[i].dX = FX * fSN;
    … …  
    380518}
    381519
    382 void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, int Nstars) {
     520void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) {
    383521
    384522  int i;
    … …  
    396534    fSN = 1.0 / sqrt(flux);
    397535
    398     stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0);
    399     stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0);
    400     stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0);
    401 
     536    stars[i].X = X[i];
     537    stars[i].Y = Y[i];
     538    stars[i].M = M[i];
     539
     540    if (ADDNOISE) {
     541      stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0);
     542      stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0);
     543      stars[i].M += fSN*rnd_gauss(0.0, 1.0);
     544    }
     545
     546    // randomly give poor PSFQF values
     547    if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) {
     548      stars[i].psfQual   = 0.25;
     549    } else {
     550      stars[i].psfQual   = PSFQUAL;
     551    }
     552   
    402553    stars[i].dX = FX * fSN;
    403554    stars[i].dY = FY * fSN;
    … …  
    413564    stars[i].fy        = FY;
    414565    stars[i].df        = DF;
    415     stars[i].psfQual   = PSFQUAL;
    416566    stars[i].nFrames   = 1;
    417     stars[i].flags     = FLAGS;
     567    stars[i].flags     = Flag[i];
    418568  }
    419569

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

    r17419 r33653  
    145145
    146146}
     147
     148// ha/dec -> alt/az
     149int altaz (double *alt, double *az, double ha, double dec, double latitude) {
     150
     151  double sind, sinh, cosh;
     152
     153  sind = dSIN (dec) * dSIN (latitude) + dCOS (dec) * dCOS (ha) * dCOS (latitude);
     154  *alt  = DEG_RAD * asin (sind);
     155
     156  sinh = - dCOS (dec) * dSIN (ha);
     157  cosh =   dSIN (dec) * dCOS (latitude) - dCOS (dec) * dCOS (ha) * dSIN (latitude);
     158
     159  *az = DEG_RAD * atan2 (sinh, cosh);
     160
     161  // I may need to use the parallactic angle to get the right plate rotation...
     162  // sinh = -dCOS(az) * dSIN(alt) * dSIN(ha) * dSIN(latitude) + dSIN(az) * dSIN(alt) * dCOS(ha) - dSIN(ha) * dCOS(alt) * dCOS(latitude);
     163  // cosh = -dSIN(az) * dSIN(ha) * dSIN(latitude) - dCOS(az) * dCOS(ha);
     164  // rot = -DEG_RAD * atan2 (sinh, cosh);
     165  return TRUE;
     166}
     167

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

    r31160 r33653  
    2828  return;
    2929}
    30 
    31 # define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
    32 # define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } }
    3330
    3431// sort the measure Sequence based on the average Sequence entries

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

    r31239 r33653  
    546546}
    547547
    548 // define the parameters of a single sky projection center
     548// define the parameters of a projection centers for this ring
     549// dec : ~ center of ring in Dec
     550// dDEC : approximate height
     551// nring : number of cells generated for this ring
     552// format : guide to generate the filenames (c-type string format)
    549553SkyRectangle *sky_rectangle_ring (float dec, float dDEC, int *nring, char *format) {
    550554

  • trunk/Ohana/src/addstar/test/simple.dvo

    r27582 r33653  
    88  test.fields PS1_V1    PS1_V1
    99  test.fields PS1_V2    PS1_V1
     10  test.fields PS1_V3    PS1_V1
    1011
    1112  test.fields PS1_DEV_0 PS1_V2
    … …  
    1314  test.fields PS1_V1    PS1_V2
    1415  test.fields PS1_V2    PS1_V2
     16  test.fields PS1_V3    PS1_V2
     17
     18  test.fields PS1_DEV_0 PS1_V3
     19  test.fields PS1_DEV_1 PS1_V3
     20  test.fields PS1_V1    PS1_V3
     21  test.fields PS1_V2    PS1_V3
     22  test.fields PS1_V3    PS1_V3
    1523end 
    1624
    … …  
    3139  mkinput
    3240  exec mkcmf test.in.txt test.cmf -date 2008/1/1 -time 01:00:00 -radec $RA $DEC -type $1
    33   exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2
     41  exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass
    3442
    3543  list testfields -copy testfields_$1
    … …  
    220228  N_FRAMES          : SKIP # no accessor       
    221229end
     230
     231# this list is good for PS1_V3
     232list testfields_PS1_V3
     233  IPP_IDET          : detid     
     234  X_PSF             : xccd     
     235  Y_PSF             : yccd     
     236  X_PSF_SIG         : xccd:err # FAIL
     237  Y_PSF_SIG         : yccd:err # FAIL
     238  RA_PSF            : SKIP # astrometry is not calibrated in the cmf
     239  DEC_PSF           : SKIP # astrometry is not calibrated in the cmf
     240  POSANGLE          : SKIP # astrometry is not calibrated in the cmf
     241  PLTSCALE          : SKIP # astrometry is not calibrated in the cmf
     242  PSF_INST_MAG      : mag:inst 
     243  PSF_INST_MAG_SIG  : mag:err   
     244  AP_MAG_STANDARD   : mag:ap # FAIL
     245  AP_MAG_RADIUS     : SKIP # no accessor
     246  PEAK_FLUX_AS_MAG  : SKIP # no accessor
     247  CAL_PSF_MAG       : SKIP # photometry is not calibrated in the cmf
     248  CAL_PSF_MAG_SIG   : SKIP # photometry is not calibrated in the cmf
     249  SKY               : sky       
     250  SKY_SIG           : sky_err   
     251  PSF_CHISQ         : psf_chisq
     252  CR_NSIGMA         : cr_nsigma
     253  EXT_NSIGMA        : ext_nsigma
     254  PSF_MAJOR         : FWHM_MAJ 
     255  PSF_MINOR         : FWHM_MIN 
     256  PSF_THETA         : THETA # FAIL
     257  PSF_QF            : PSF_QF   
     258  PSF_NDOF          : SKIP # no accessor
     259  PSF_NPIX          : SKIP # no accessor
     260  MOMENTS_XX        : SKIP # no accessor
     261  MOMENTS_XY        : SKIP # no accessor
     262  MOMENTS_YY        : SKIP # no accessor
     263  FLAGS             : phot_flags
     264  N_FRAMES          : SKIP # no accessor       
     265end

  • trunk/Ohana/src/addstar/test/tap.dvo

    r27435 r33653  
    11# -*-sh-*-
     2
     3if (not($?TAP_BREAK)) set TAP_BREAK = 0
    24
    35macro tapOK
    … …  
    2931  $TAP_NTEST = $1
    3032  $TAP_NFAIL = 0
     33  $TAP_NSKIP = 0
    3134  $TAP_NDONE = 0
    3235  if (not($?TAP_BREAK)) set TAP_BREAK = 0
     36end
     37
     38macro tapSKIP
     39  if ($0 != 2)
     40    echo "USAGE: tapSKIP (nskip)"
     41    break
     42  end
     43
     44  $TAP_NSKIP ++
     45  $TAP_NDONE ++
    3346end
    3447
    … …  
    4760  end
    4861
    49   if ($TAP_NFAIL == 0)
     62  if ($TAP_NSKIP)
     63    echo "skipped $TAP_NSKIP of $TAP_NDONE"
     64  end
     65
     66  if ($TAP_NFAIL + $TAP_NSKIP == 0)
    5067    echo "passed $TAP_NDONE tests"
    5168  end
Note: See TracChangeset for help on using the changeset viewer.