IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33391


Ignore:
Timestamp:
Mar 1, 2012, 8:49:24 PM (14 years ago)
Author:
eugene
Message:

add photFlags to MeasureTiny, pass along to BrightStars; updates to the relphot.reject.dvo test set; tests prints for ubercalDist; add rejection of bad and poor stars in bcatalog (make "poor" optional?); reject outlier measures based on hypot(star_sigma, meas_sigma)

Location:
branches/eam_branches/ipp-20111122/Ohana/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20111122/Ohana/src/addstar/test/relphot.reject.dvo

    r33388 r33391  
    3434  $CHIP_DY = 1000
    3535  $DEFECT_FRAC = 0.03
     36  $OFFSET_FRAC_UC = 0.1
     37  $OFFSET_FRAC_NC = 0.5
    3638  if (not($?VERBOSE)) set VERBOSE = 0
    3739
    … …  
    5456  # sequence to count images (only used in this function)
    5557  # XXX  TEST create tmpseq 0 9
    56   create tmpseq 0 3
     58  create tmpseq 0 9
    5759
    5860  ##### mjd and zpt values for ubercal'ed images
    … …  
    6365  $filt_uc:1 = g;  filtN_uc[1] = 1
    6466  $filt_uc:2 = g;  filtN_uc[2] = 1
    65   # $filt_uc:3 = r;  filtN_uc[3] = 2
    66   # $filt_uc:4 = r;  filtN_uc[4] = 2
    67   # $filt_uc:5 = r;  filtN_uc[5] = 2
    68   # $filt_uc:6 = i;  filtN_uc[6] = 3
    69   # $filt_uc:7 = i;  filtN_uc[7] = 3
    70   # $filt_uc:8 = i;  filtN_uc[8] = 3
     67  $filt_uc:3 = g;  filtN_uc[3] = 1
     68  $filt_uc:4 = g;  filtN_uc[4] = 1
     69  $filt_uc:5 = g;  filtN_uc[5] = 1
     70  $filt_uc:6 = g;  filtN_uc[6] = 1
     71  $filt_uc:7 = g;  filtN_uc[7] = 1
     72  $filt_uc:8 = g;  filtN_uc[8] = 1
    7173
    7274  # airmass slopes for this sequnece
    … …  
    8789  mjd_uc[1] = 55000.02
    8890  mjd_uc[2] = 55000.03
    89   # mjd_uc[3] = 55015.01
    90   # mjd_uc[4] = 55015.02
    91   # mjd_uc[5] = 55015.03
    92   # mjd_uc[6] = 55025.01
    93   # mjd_uc[7] = 55025.02
    94   # mjd_uc[8] = 55025.03
     91  mjd_uc[3] = 55015.01
     92  mjd_uc[4] = 55015.02
     93  mjd_uc[5] = 55015.03
     94  mjd_uc[6] = 55025.01
     95  mjd_uc[7] = 55025.02
     96  mjd_uc[8] = 55025.03
    9597
    9698  ##### mjd and zpt values for not-ubercal'ed images
    … …  
    101103  $filt_nc:1 = g;  filtN_nc[1] = 1
    102104  $filt_nc:2 = g;  filtN_nc[2] = 1
    103   # $filt_nc:3 = r;  filtN_nc[3] = 2
    104   # $filt_nc:4 = r;  filtN_nc[4] = 2
    105   # $filt_nc:5 = r;  filtN_nc[5] = 2
    106   # $filt_nc:6 = i;  filtN_nc[6] = 3
    107   # $filt_nc:7 = i;  filtN_nc[7] = 3
    108   # $filt_nc:8 = i;  filtN_nc[8] = 3
     105  $filt_nc:3 = g;  filtN_nc[3] = 1
     106  $filt_nc:4 = g;  filtN_nc[4] = 1
     107  $filt_nc:5 = g;  filtN_nc[5] = 1
     108  $filt_nc:6 = g;  filtN_nc[6] = 1
     109  $filt_nc:7 = g;  filtN_nc[7] = 1
     110  $filt_nc:8 = g;  filtN_nc[8] = 1
    109111
    110112  # airmass slopes for these sequnece
    … …  
    126128  mjd_nc[1] = 55000.12
    127129  mjd_nc[2] = 55000.13
    128   # mjd_nc[3] = 55015.11
    129   # mjd_nc[4] = 55015.12
    130   # mjd_nc[5] = 55015.13
    131   # mjd_nc[6] = 55025.11
    132   # mjd_nc[7] = 55025.12
    133   # mjd_nc[8] = 55025.13
     130  mjd_nc[3] = 55015.11
     131  mjd_nc[4] = 55015.12
     132  mjd_nc[5] = 55015.13
     133  mjd_nc[6] = 55025.11
     134  mjd_nc[7] = 55025.12
     135  mjd_nc[8] = 55025.13
    134136end
    135137
    … …  
    152154  init
    153155  mkzptfile
    154   mkstars 500
     156  mkstars 10000
    155157
    156158  tapPLAN {(((mjd_uc[] + mjd_nc[])*4 + (mjd_uc[] + mjd_nc[])) + (mjd_uc[] + mjd_nc[]))}
    157159
    158160  # generate the basic images and check they were correctly ingested by dvo
    159   mksequence $fileroot $catdir
     161  mksequence.drift $fileroot $catdir
    160162  for i 0 mjd_uc[]
    161163    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw
    … …  
    175177
    176178  # run relphot on the serial db and check that the images now match the expected values
    177   if ($VERBOSE) echo relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0
    178   exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0
     179  # if ($VERBOSE) echo relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0
     180  if ($VERBOSE) echo relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0
     181  exec relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0
    179182  for i 0 mjd_uc[]
    180183    ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot
    … …  
    315318end
    316319
    317 macro mksequence
     320macro mksequence.drift
     321  if ($0 != 3)
     322    echo "mksequence.drift (fileroot) (catdir)"
     323    break
     324  end
     325
     326  echo "starting mksequence..."
     327  local i T1 T2
     328
     329  ctimes -abs now -var T1
     330  # we have a second set of mjd and zero point values for non-ubercal'ed images
     331  echo -no-return "making mjd_uc[] UC images: "
     332  for i 0 mjd_uc[]
     333    echo -no-return "$i.. "
     334    mkexposure $1.uc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_UC/3600.0} zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2
     335  end
     336  echo ""
     337  echo -no-return "making mjd_nc[] NC images: "
     338  for i 0 mjd_nc[]
     339    echo -no-return "$i.. "
     340    mkexposure $1.nc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_NC/3600.0} zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2
     341  end
     342  ctimes -abs now -var T2
     343  echo "done ({86400*($T2 - $T1)} sec)"
     344end
     345
     346macro mksequence.1pt
    318347  if ($0 != 3)
    319348    echo "mksequence (fileroot) (catdir)"
    … …  
    495524  # size of region of interest in linear arcseconds
    496525  $RA_RANGE  = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE
    497   $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE
     526  $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE * (1 + tmpseq[]*$OFFSET_FRAC_NC)
    498527
    499528  create tmp 0 $1

  • branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h

    r33375 r33391  
    293293  unsigned int   imageID;
    294294  unsigned int   dbFlags;
     295  unsigned int   photFlags;
    295296  int            catID; // unsigned int?
    296297  unsigned short photcode;

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/BrightCatalog.c

    r33364 r33391  
    6262 
    6363    // need to create and assign to flat-field correction
    64     GET_COLUMN(dR,       "RA_OFF",   float);
    65     GET_COLUMN(dD,       "DEC_OFF",  float);
    66     GET_COLUMN(M,        "MAG_SYS",  float);
    67     GET_COLUMN(Mcal,     "MAG_CAL",  float);
    68     GET_COLUMN(dM,       "MAG_ERR",  float);
    69     GET_COLUMN(airmass,  "AIRMASS",  float);
    70     GET_COLUMN(Xccd,     "X_CCD",    float);
    71     GET_COLUMN(Yccd,     "Y_CCD",    float);
    72     GET_COLUMN(dt,       "EXPTIME",  float);
    73     GET_COLUMN(t,        "TIME",     int);
    74     GET_COLUMN(averef,   "AVE_REF",  int); // XXX signed vs unsigned?
    75     GET_COLUMN(imageID,  "IMAGE_ID", int);
    76     GET_COLUMN(dbFlags,  "DB_FLAGS", int);
    77     GET_COLUMN(catID,    "CAT_ID",   int);
    78     GET_COLUMN(photcode, "PHOTCODE", short);
     64    GET_COLUMN(dR,        "RA_OFF",     float);
     65    GET_COLUMN(dD,        "DEC_OFF",    float);
     66    GET_COLUMN(M,         "MAG_SYS",    float);
     67    GET_COLUMN(Mcal,      "MAG_CAL",    float);
     68    GET_COLUMN(dM,        "MAG_ERR",    float);
     69    GET_COLUMN(airmass,   "AIRMASS",    float);
     70    GET_COLUMN(Xccd,      "X_CCD",      float);
     71    GET_COLUMN(Yccd,      "Y_CCD",      float);
     72    GET_COLUMN(dt,        "EXPTIME",    float);
     73    GET_COLUMN(t,         "TIME",       int);
     74    GET_COLUMN(averef,    "AVE_REF",    int); // XXX signed vs unsigned?
     75    GET_COLUMN(imageID,   "IMAGE_ID",   int);
     76    GET_COLUMN(dbFlags,   "DB_FLAGS",   int);
     77    GET_COLUMN(photFlags, "PHOT_FLAGS", int);
     78    GET_COLUMN(catID,     "CAT_ID",     int);
     79    GET_COLUMN(photcode,  "PHOTCODE",   short);
    7980    // XXX free the fits table data here
    8081
    … …  
    8283    ALLOCATE (measure, MeasureTiny, Nrow);
    8384    for (i = 0; i < Nrow; i++) {
    84       measure[i].dR       = dR[i];
    85       measure[i].dD       = dD[i];
    86       measure[i].M        = M[i];
    87       measure[i].Mcal     = Mcal[i];
    88       measure[i].dM       = dM[i];
    89       measure[i].airmass  = airmass[i];
    90       measure[i].Xccd     = Xccd[i];
    91       measure[i].Yccd     = Yccd[i];
    92       measure[i].dt       = dt[i];
    93       measure[i].t        = t[i];
    94       measure[i].averef   = averef[i];
    95       measure[i].imageID  = imageID[i];
    96       measure[i].dbFlags  = dbFlags[i];
    97       measure[i].catID    = catID[i];
    98       measure[i].photcode = photcode[i];
     85      measure[i].dR        = dR[i];
     86      measure[i].dD        = dD[i];
     87      measure[i].M         = M[i];
     88      measure[i].Mcal      = Mcal[i];
     89      measure[i].dM        = dM[i];
     90      measure[i].airmass   = airmass[i];
     91      measure[i].Xccd      = Xccd[i];
     92      measure[i].Yccd      = Yccd[i];
     93      measure[i].dt        = dt[i];
     94      measure[i].t         = t[i];
     95      measure[i].averef    = averef[i];
     96      measure[i].imageID   = imageID[i];
     97      measure[i].dbFlags   = dbFlags[i];
     98      measure[i].photFlags = photFlags[i];
     99      measure[i].catID     = catID[i];
     100      measure[i].photcode  = photcode[i];
    99101    }
    100102    // fprintf (stderr, "loaded data for %lld measures\n", (long long) Nrow);
    … …  
    113115    free (imageID );
    114116    free (dbFlags );
     117    free (photFlags );
    115118    free (catID   );
    116119    free (photcode);
    … …  
    261264
    262265    // create intermediate storage arrays
    263     float *dR        ; ALLOCATE (dR      ,  float, catalog->Nmeasure);
    264     float *dD        ; ALLOCATE (dD      ,  float, catalog->Nmeasure);
    265     float *M         ; ALLOCATE (M       ,  float, catalog->Nmeasure);
    266     float *Mcal      ; ALLOCATE (Mcal    ,  float, catalog->Nmeasure);
    267     float *dM        ; ALLOCATE (dM      ,  float, catalog->Nmeasure);
    268     float *airmass   ; ALLOCATE (airmass ,  float, catalog->Nmeasure);
    269     float *Xccd      ; ALLOCATE (Xccd    ,  float, catalog->Nmeasure);
    270     float *Yccd      ; ALLOCATE (Yccd    ,  float, catalog->Nmeasure);
    271     float *dt        ; ALLOCATE (dt      ,  float, catalog->Nmeasure);
    272     int   *t         ; ALLOCATE (t       ,  int  , catalog->Nmeasure);
    273     int   *averef    ; ALLOCATE (averef  ,  int  , catalog->Nmeasure);
    274     int   *imageID   ; ALLOCATE (imageID ,  int  , catalog->Nmeasure);
    275     int   *dbFlags   ; ALLOCATE (dbFlags ,  int  , catalog->Nmeasure);
    276     int   *catID     ; ALLOCATE (catID   ,  int  , catalog->Nmeasure);
    277     short *photcode  ; ALLOCATE (photcode,  short, catalog->Nmeasure);
     266    float *dR        ; ALLOCATE (dR       ,  float, catalog->Nmeasure);
     267    float *dD        ; ALLOCATE (dD       ,  float, catalog->Nmeasure);
     268    float *M         ; ALLOCATE (M        ,  float, catalog->Nmeasure);
     269    float *Mcal      ; ALLOCATE (Mcal     ,  float, catalog->Nmeasure);
     270    float *dM        ; ALLOCATE (dM       ,  float, catalog->Nmeasure);
     271    float *airmass   ; ALLOCATE (airmass  ,  float, catalog->Nmeasure);
     272    float *Xccd      ; ALLOCATE (Xccd     ,  float, catalog->Nmeasure);
     273    float *Yccd      ; ALLOCATE (Yccd     ,  float, catalog->Nmeasure);
     274    float *dt        ; ALLOCATE (dt       ,  float, catalog->Nmeasure);
     275    int   *t         ; ALLOCATE (t        ,  int  , catalog->Nmeasure);
     276    int   *averef    ; ALLOCATE (averef   ,  int  , catalog->Nmeasure);
     277    int   *imageID   ; ALLOCATE (imageID  ,  int  , catalog->Nmeasure);
     278    int   *dbFlags   ; ALLOCATE (dbFlags  ,  int  , catalog->Nmeasure);
     279    int   *photFlags ; ALLOCATE (photFlags,  int  , catalog->Nmeasure);
     280    int   *catID     ; ALLOCATE (catID    ,  int  , catalog->Nmeasure);
     281    short *photcode  ; ALLOCATE (photcode ,  short, catalog->Nmeasure);
    278282
    279283    // assign the storage arrays
    … …  
    294298      imageID[i]  = measure[i].imageID  ;
    295299      dbFlags[i]  = measure[i].dbFlags  ;
     300      photFlags[i]= measure[i].photFlags;
    296301      photcode[i] = measure[i].photcode ;
    297302    }
    298303
    299304    // add the columns to the output array
    300     gfits_set_bintable_column (&theader, &ftable, "RA_OFF",   dR,       catalog->Nmeasure);
    301     gfits_set_bintable_column (&theader, &ftable, "DEC_OFF",  dD,       catalog->Nmeasure);
    302     gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",  M,        catalog->Nmeasure);
    303     gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",  Mcal,     catalog->Nmeasure);
    304     gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM,       catalog->Nmeasure);
    305     gfits_set_bintable_column (&theader, &ftable, "AIRMASS",  airmass,  catalog->Nmeasure);
    306     gfits_set_bintable_column (&theader, &ftable, "X_CCD",    Xccd,     catalog->Nmeasure);
    307     gfits_set_bintable_column (&theader, &ftable, "Y_CCD",    Yccd,     catalog->Nmeasure);
    308     gfits_set_bintable_column (&theader, &ftable, "EXPTIME",  dt,       catalog->Nmeasure);
    309     gfits_set_bintable_column (&theader, &ftable, "TIME",     t,        catalog->Nmeasure);
    310     gfits_set_bintable_column (&theader, &ftable, "AVE_REF",  averef,   catalog->Nmeasure);
    311     gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID,  catalog->Nmeasure);
    312     gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags,  catalog->Nmeasure);
    313     gfits_set_bintable_column (&theader, &ftable, "CAT_ID",   catID,    catalog->Nmeasure);
    314     gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure);
     305    gfits_set_bintable_column (&theader, &ftable, "RA_OFF",     dR,        catalog->Nmeasure);
     306    gfits_set_bintable_column (&theader, &ftable, "DEC_OFF",    dD,        catalog->Nmeasure);
     307    gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
     308    gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
     309    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",    dM,        catalog->Nmeasure);
     310    gfits_set_bintable_column (&theader, &ftable, "AIRMASS",    airmass,   catalog->Nmeasure);
     311    gfits_set_bintable_column (&theader, &ftable, "X_CCD",      Xccd,      catalog->Nmeasure);
     312    gfits_set_bintable_column (&theader, &ftable, "Y_CCD",      Yccd,      catalog->Nmeasure);
     313    gfits_set_bintable_column (&theader, &ftable, "EXPTIME",    dt,        catalog->Nmeasure);
     314    gfits_set_bintable_column (&theader, &ftable, "TIME",       t,         catalog->Nmeasure);
     315    gfits_set_bintable_column (&theader, &ftable, "AVE_REF",    averef,    catalog->Nmeasure);
     316    gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID",   imageID,   catalog->Nmeasure);
     317    gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS",   dbFlags,   catalog->Nmeasure);
     318    gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure);
     319    gfits_set_bintable_column (&theader, &ftable, "CAT_ID",     catID,     catalog->Nmeasure);
     320    gfits_set_bintable_column (&theader, &ftable, "PHOTCODE",   photcode,  catalog->Nmeasure);
    315321
    316322    free (dR      );
    … …  
    327333    free (imageID );
    328334    free (dbFlags );
     335    free (photFlags );
    329336    free (catID   );
    330337    free (photcode);

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c

    r33374 r33391  
    456456      }
    457457
    458       minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt].M_80, minUbercalDist); // M_80 == ubercalDist
     458      PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     459      if (!code) goto skip;
     460      if (code->equiv < 1) goto skip;
     461      int Nsec = GetPhotcodeNsec (code->equiv);
     462      if (Nsec == -1) goto skip;
     463      minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].M_80, minUbercalDist); // M_80 == ubercalDist
    459464
    460465      // fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat);
    461466
     467    skip:
    462468      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
    463469      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    … …  
    497503
    498504    image[i].dummy3 = minUbercalDist;
     505    fprintf (stderr, "%s : %d\n", image[i].name, image[i].dummy3);
    499506  }
    500507  free (list);

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c

    r33374 r33391  
    5656  int Nsecfilt = GetPhotcodeNsecfilt ();
    5757  Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
     58
     59  FILE *f = fopen ("few.v1.dat", "w");
    5860
    5961  for (i = 0; i < Ncatalog; i++) {
    … …  
    135137        // when performing the grid analysis, STAR_TOOFEW will be set to 1;
    136138        if (N <= STAR_TOOFEW) { /* too few measurements */
     139          fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[catalog[i].averageT[j].measureOffset].imageID, N, STAR_TOOFEW);
    137140          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
    138141          Nfew ++;
    … …  
    151154    }
    152155  }
     156  fclose (f);
    153157  fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys);
    154158
    … …  
    560564  TOOFEW = MAX (5, STAR_TOOFEW);
    561565
     566  FILE *f = fopen ("poor.dat", "w");
     567
    562568  Ndel = Nave = 0;
    563569  Ncal = Nmos = Ngrid = Nfew = 0;
    … …  
    598604          N++;
    599605        }
    600         if (N <= TOOFEW) { Nfew ++; continue; }
     606        if (N <= TOOFEW) {
     607          Nfew ++;
     608          continue;
     609        }
    601610
    602611        /* 3-sigma clip based on stats of inner 50% */
    … …  
    648657        /* mark bad measures (> 3 sigma deviant) */
    649658        for (k = 0; k < N; k++) {
    650           if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
     659          // treat the scatter of the star as a systematic term.  this is a bit of an
     660          // over-estimage (a perfect Gauss distribution with perfect errors would have
     661          // mySigma = sqrt(2) too large)
     662          float mySigma = hypot (stats.sigma, dlist[k]);
     663          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
     664            fprintf (f, "%10.6f %10.6f %d %f %f %f %f 0x%08x\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[ilist[k]].imageID, list[k], stats.median, stats.sigma, mySigma, catalog[i].measureT[ilist[k]].photFlags);
    651665            catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    652666            if (final) {
    … …  
    655669            }
    656670            Ndel ++;
     671          } else {
     672            if (catalog[i].measureT[ilist[k]].photFlags & 0x02) {
     673              fprintf (stderr, "failed to mark %10.6f %10.6f %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[ilist[k]].imageID);
     674            }
    657675          }
    658676        }
    … …  
    661679    }
    662680  }
     681  fclose (f);
    663682  initstats (STATMODE);
    664683  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
    … …  
    667686  free (ilist);
    668687  free (tlist);
     688
     689
    669690}
    670691

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c

    r33372 r33391  
    1818int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) {
    1919
    20   measureT[0].dR       = measure[0].dR;
    21   measureT[0].dD       = measure[0].dD;
    22   measureT[0].M        = measure[0].M;
    23   measureT[0].Mcal     = measure[0].Mcal;
    24   measureT[0].dM       = measure[0].dM;
    25   measureT[0].airmass  = measure[0].airmass;
    26   measureT[0].Xccd     = measure[0].Xccd;
    27   measureT[0].Yccd     = measure[0].Yccd;
    28   measureT[0].t        = measure[0].t;
    29   measureT[0].dt       = measure[0].dt;
    30   measureT[0].averef   = measure[0].averef;
    31   measureT[0].imageID  = measure[0].imageID;
    32   measureT[0].dbFlags  = measure[0].dbFlags;
    33   measureT[0].photcode = measure[0].photcode;
     20  measureT[0].dR        = measure[0].dR;
     21  measureT[0].dD        = measure[0].dD;
     22  measureT[0].M         = measure[0].M;
     23  measureT[0].Mcal      = measure[0].Mcal;
     24  measureT[0].dM        = measure[0].dM;
     25  measureT[0].airmass   = measure[0].airmass;
     26  measureT[0].Xccd      = measure[0].Xccd;
     27  measureT[0].Yccd      = measure[0].Yccd;
     28  measureT[0].t         = measure[0].t;
     29  measureT[0].dt        = measure[0].dt;
     30  measureT[0].averef    = measure[0].averef;
     31  measureT[0].imageID   = measure[0].imageID;
     32  measureT[0].dbFlags   = measure[0].dbFlags;
     33  measureT[0].photFlags = measure[0].photFlags;
     34  measureT[0].photcode  = measure[0].photcode;
    3435
    3536  return (TRUE);
    … …  
    3940 
    4041  off_t i, j, offset;
    41   int ecode, found, Ns;
     42  int found, Ns;
    4243  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    4344  float mag;
    44   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan;
     45  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
    4546
    4647  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    5657  Nmeasure = Naverage = 0;
    5758
    58   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = 0;
     59  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
    5960
    6061  // copy the following fields to the subcatalog:
    … …  
    105106
    106107      /* select measurements by photcode */
    107       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode);
     108      PhotCode *code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
     109      if (!code) continue;
    108110      found = FALSE;
    109111      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
    110         if (ecode == photcodes[Ns][0].code) found = TRUE;
     112        if (code->equiv == photcodes[Ns][0].code) found = TRUE;
    111113      }
    112114      if (!found) {
    … …  
    127129      if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; }
    128130      if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
     131
     132      if (catalog[0].measure[offset].photFlags & code->photomBadMask) { Nbad++; continue;}
     133      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
    129134
    130135      // check for galaxies
    … …  
    206211    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n",
    207212             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
    208     fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies\n",
    209              Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy);
     213    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     214             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
    210215  }
    211216  return (TRUE);
    … …  
    336341  for (i = 0; i < catalog[0].Nmeasure; i++) {
    337342    // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
    338     measureT[i].dR       = measure[i].dR;
    339     measureT[i].dD       = measure[i].dD;
    340     measureT[i].M        = measure[i].M;
    341     measureT[i].Mcal     = measure[i].Mcal;
    342     measureT[i].dM       = measure[i].dM;
    343     measureT[i].airmass  = measure[i].airmass;
    344     measureT[i].Xccd     = measure[i].Xccd;
    345     measureT[i].Yccd     = measure[i].Yccd;
    346     measureT[i].t        = measure[i].t;
    347     measureT[i].dt       = measure[i].dt;
    348     measureT[i].averef   = measure[i].averef;
    349     measureT[i].imageID  = measure[i].imageID;
    350     measureT[i].dbFlags  = measure[i].dbFlags;
    351     measureT[i].photcode = measure[i].photcode;
     343    measureT[i].dR        = measure[i].dR;
     344    measureT[i].dD        = measure[i].dD;
     345    measureT[i].M         = measure[i].M;
     346    measureT[i].Mcal      = measure[i].Mcal;
     347    measureT[i].dM        = measure[i].dM;
     348    measureT[i].airmass   = measure[i].airmass;
     349    measureT[i].Xccd      = measure[i].Xccd;
     350    measureT[i].Yccd      = measure[i].Yccd;
     351    measureT[i].t         = measure[i].t;
     352    measureT[i].dt        = measure[i].dt;
     353    measureT[i].averef    = measure[i].averef;
     354    measureT[i].imageID   = measure[i].imageID;
     355    measureT[i].dbFlags   = measure[i].dbFlags;
     356    measureT[i].photFlags = measure[i].photFlags;
     357    measureT[i].photcode  = measure[i].photcode;
    352358  }
    353359

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/liststats.c

    r16040 r33391  
    106106    dS += M;
    107107  }
    108   X2 = X2 / Nm;
     108  X2 = X2 / (Nm - 1);
    109109  dS = sqrt (dS / Nm);
    110110

  • branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c

    r33387 r33391  
    77void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    88
    9   off_t i, j, m;
     9  off_t i;
    1010  int ecode;
    1111
    … …  
    3636
    3737      for (i = 0; i < catalog[0].Naverage; i++) {
    38         catalog[0].secfilt[Nsecfilt*i+Nsec].M  = NAN;
    39         catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN;
    40         catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT;
    41         catalog[0].secfilt[Nsecfilt*j+Nsec].Ncode = 0;
    42         catalog[0].secfilt[Nsecfilt*j+Nsec].Nused = 0;
    43 
    44         catalog[0].secfilt[Nsecfilt*i+Nsec].flags &= ~PHOTOM_PASSES;
    45 
    46         m = catalog[0].average[i].measureOffset;
     38        off_t N = Nsecfilt*i+Nsec;
     39        catalog[0].secfilt[N].M  = NAN;
     40        catalog[0].secfilt[N].dM = NAN;
     41        catalog[0].secfilt[N].Xm = NAN_S_SHORT;
     42        catalog[0].secfilt[N].Ncode = 0;
     43        catalog[0].secfilt[N].Nused = 0;
     44        catalog[0].secfilt[N].flags &= ~PHOTOM_PASSES;
     45
     46        off_t m = catalog[0].average[i].measureOffset;
     47        off_t j;
    4748        for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    4849       
    … …  
    8990    catalog[0].average[i].flags &= ~ID_STAR_FEW;
    9091    catalog[0].average[i].flags &= ~ID_STAR_POOR;
     92    off_t j;
    9193    for (j = 0; j < Nsecfilt; j++) {
    9294        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW;
    9395        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR;
    9496    }
    95     m = catalog[0].average[i].measureOffset;
     97    off_t m = catalog[0].average[i].measureOffset;
    9698    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    9799      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
Note: See TracChangeset for help on using the changeset viewer.