IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33827


Ignore:
Timestamp:
Apr 26, 2012, 12:30:42 PM (14 years ago)
Author:
eugene
Message:

nFitPhotom was signed, changed to unsigned short; change the iterations to have at least 8 loops before any cleaning; add some verbosity; pull the old (pre-threaded) version of setMmos; setMmos_mosaic should take a pointer to the current mosaic; skip entries with skipCal set

Location:
branches/eam_branches/ipp-20120405/Ohana/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20120405/Ohana/src/libautocode/def/image.d

    r33647 r33827  
    6969
    7070FIELD     dMagSys,          MAG_SYS_ERR,          float,          systematic photometry error
    71 FIELD     nFitAstrom,       N_FIT_ASTROM,         short,          number of stars used for astrometry cal
    72 FIELD     nFitPhotom,       N_FIT_PHOTOM,         short,          number of stars used for photometry cal
     71FIELD     nFitAstrom,       N_FIT_ASTROM,         unsigned short, number of stars used for astrometry cal
     72FIELD     nFitPhotom,       N_FIT_PHOTOM,         unsigned short, number of stars used for photometry cal
    7373
    7474FIELD     photom_map_id,    PHOTOM_MAP_ID,        unsigned int,   reference to 2D zero point map

  • branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.data/gridify.c

    r31450 r33827  
    55  int i, Nx, Ny, Xb, Yb, Normalize, N;
    66  float Xmin, Xmax, dX, Ymin, Ymax, dY, initValue;
    7   float *buf, *val;
     7  float *buf, *val, *cnt;
    88  int *Nval;
    9   Buffer *bf;
    109  Vector *vx, *vy, *vz;
    1110  opihi_flt *x, *y, *z;
     11
     12  Buffer *bf = NULL;
     13  Buffer *ct = NULL;
    1214
    1315  Normalize = TRUE;
    … …  
    1517    remove_argument (N, &argc, argv);
    1618    Normalize = FALSE;
     19    if ((ct = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     20    remove_argument (N, &argc, argv);
    1721  }
    1822
    … …  
    2428  }
    2529
    26   if (argc != 11) {
    27     gprint (GP_ERR, "USAGE: gridify x y z buffer Xmin Xmax dX Ymin Ymax dY\n");
     30  Xmin = Xmax = dX = NAN;
     31  if ((N = get_argument (argc, argv, "-x"))) {
     32    remove_argument (N, &argc, argv);
     33    Xmin = atof (argv[N]);
     34    remove_argument (N, &argc, argv);
     35    Xmax = atof (argv[N]);
     36    remove_argument (N, &argc, argv);
     37    dX   = atof (argv[N]);
     38    remove_argument (N, &argc, argv);
     39  }   
     40
     41  Ymin = Ymax = dY = NAN;
     42  if ((N = get_argument (argc, argv, "-y"))) {
     43    remove_argument (N, &argc, argv);
     44    Ymin = atof (argv[N]);
     45    remove_argument (N, &argc, argv);
     46    Ymax = atof (argv[N]);
     47    remove_argument (N, &argc, argv);
     48    dY   = atof (argv[N]);
     49    remove_argument (N, &argc, argv);
     50  }   
     51
     52  if (argc != 5) {
     53    gprint (GP_ERR, "USAGE: gridify x y z buffer [-x Xmin Xmax dX] [-y Ymin Ymax dY] [-init-value value] [-raw]\n");
    2854    return (FALSE);
    2955  }
    … …  
    4167  REQUIRE_VECTOR_FLT (vz, FALSE);
    4268
    43   Xmin = atof (argv[5]);
    44   Xmax = atof (argv[6]);
    45   dX   = atof (argv[7]);
     69  if (isnan(dX)) {
     70    Xmin = 0;
     71    Xmax = bf[0].matrix.Naxis[0];
     72    dX = 1;
     73  }
    4674
    47   Ymin = atof (argv[8]);
    48   Ymax = atof (argv[9]);
    49   dY   = atof (argv[10]);
     75  if (isnan(dY)) {
     76    Ymin = 0;
     77    Ymax = bf[0].matrix.Naxis[1];
     78    dY = 1;
     79  }
    5080
    51   Nx = (Xmax - Xmin) / dX + 1;
    52   Ny = (Ymax - Ymin) / dY + 1;
     81  Nx = (Xmax - Xmin) / dX;
     82  Ny = (Ymax - Ymin) / dY;
    5383 
    54   gfits_free_matrix (&bf[0].matrix);
    55   gfits_free_header (&bf[0].header);
    56   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
    57   strcpy (bf[0].file, "(empty)");
     84  if ((Nx != bf[0].matrix.Naxis[0]) || (Ny != bf[0].matrix.Naxis[1])) {
     85    gfits_free_matrix (&bf[0].matrix);
     86    gfits_free_header (&bf[0].header);
     87    CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     88    strcpy (bf[0].file, "(empty)");
     89  }
    5890
    5991  ALLOCATE (val, float, Nx*Ny);
    … …  
    68100    Xb = (*x - Xmin) / dX;
    69101    Yb = (*y - Ymin) / dY;
     102    if (Xb < 0) continue;
     103    if (Yb < 0) continue;
    70104    if (Xb >= Nx) continue;
    71105    if (Yb >= Ny) continue;
    … …  
    74108  }
    75109
     110  if (!Normalize) {
     111    gfits_free_matrix (&ct[0].matrix);
     112    gfits_free_header (&ct[0].header);
     113    CreateBuffer (ct, Nx, Ny, -32, 0.0, 1.0);
     114    strcpy (ct[0].file, "(empty)");
     115
     116    buf = (float *) bf[0].matrix.buffer;
     117    cnt = (float *) ct[0].matrix.buffer;
     118    for (i = 0; i < Nx*Ny; i++) {
     119      if (Nval[i] == 0) continue;
     120      buf[i] = val[i];
     121      cnt[i] = Nval[i];
     122    }
     123    free (val);
     124    free (Nval);
     125    return TRUE;
     126  }
     127
    76128  buf = (float *) bf[0].matrix.buffer;
    77129  for (i = 0; i < Nx*Ny; i++) {
    78130    buf[i] = initValue;
    79     if (Normalize) {
    80       if (Nval[i] == 0) {
    81         continue;
    82       }
    83       buf[i] = val[i] / Nval[i];
    84     } else {
    85       buf[i] = val[i];
    86     }
     131    if (Nval[i] == 0) continue;
     132    buf[i] = val[i] / Nval[i];
    87133  }
    88134

  • branches/eam_branches/ipp-20120405/Ohana/src/opihi/dvo/avextract.c

    r33662 r33827  
    3535    remove_argument (N, &argc, argv);
    3636    VERBOSE = TRUE;
     37  }
     38
     39  int VERBOSE2 = FALSE;
     40  if ((N = get_argument (argc, argv, "-vv"))) {
     41    remove_argument (N, &argc, argv);
     42    VERBOSE = TRUE;
     43    VERBOSE2 = TRUE;
    3744  }
    3845
    … …  
    173180    catalog.filename = (HOST_ID || PARALLEL_LOCAL) ? hostfile : skylist[0].filename[i];
    174181    catalog.catflags = LOAD_AVES | LOAD_SECF;
    175     if (needMeasures) {
    176       catalog.catflags |= LOAD_MEAS;
    177     }
     182    catalog.catflags |= needMeasures ? LOAD_MEAS : SKIP_MEAS;
    178183    catalog.Nsecfilt = 0;
    179184
    … …  
    181186     
    182187    // an error exit status here is a significant error
    183     if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
     188    if (!dvo_catalog_open (&catalog, NULL, VERBOSE2, "r")) {
    184189      gprint (GP_ERR, "ERROR: failure to open catalog file %s\n", catalog.filename);
    185190      exit (2);

  • branches/eam_branches/ipp-20120405/Ohana/src/opihi/dvo/avmatch.c

    r33662 r33827  
    129129  }
    130130
     131  // check the requested fields : are all average/secfilt entries, or do we need measures?
     132  int needMeasures = FALSE;
     133  for (i = 0; !needMeasures && (i < Nfields); i++) {
     134    if (fields[i].magMode == MAG_NONE) continue;
     135    if (fields[i].photcode == NULL) continue; // assert this?
     136    if (fields[i].photcode[0].type == PHOT_REF) needMeasures = TRUE;
     137    if (fields[i].photcode[0].type == PHOT_DEP) needMeasures = TRUE;
     138  }
     139
    131140  /* load regions which contain all supplied RA,DEC coordinates */
    132141  if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape;
    … …  
    168177    snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name);
    169178    catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
    170     catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
     179    catalog.catflags = LOAD_AVES | LOAD_SECF;
     180    catalog.catflags |= needMeasures ? LOAD_MEAS : SKIP_MEAS;
    171181    catalog.Nsecfilt = 0;
    172182

  • branches/eam_branches/ipp-20120405/Ohana/src/opihi/dvo/gstar.c

    r33662 r33827  
    116116
    117117  /* lock, load, unlock catalog */
    118   catalog.catflags = GetMeasures ? LOAD_AVES | LOAD_MEAS | LOAD_SECF : LOAD_AVES | LOAD_SECF;
     118  catalog.catflags = LOAD_AVES | LOAD_SECF;
     119  catalog.catflags |= GetMeasures ? LOAD_MEAS : SKIP_MEAS;
    119120  catalog.Nsecfilt = 0;
    120121

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h

    r33823 r33827  
    2727  float dMcal;
    2828  float dMsys;
    29   short nFitPhotom;
     29  unsigned short nFitPhotom;
    3030  short Xm;
    3131  float secz;

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c

    r33820 r33827  
    228228    mosaic[i].photcode = 0;
    229229    mosaic[i].skipCal = FALSE;
     230   
     231    memset (&mosaic[i].coords, 0, sizeof(Coords));
    230232
    231233    MosaicN_IMAGE[i] = 10;
    … …  
    709711}
    710712
    711 # if (0)
    712 int setMmos_old (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    713 
    714   off_t i, j, m, c, n, N, Nmax;
    715   int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright;
    716   float Msys, Mrel, Mcal, Mgrid, Mflat;
    717   double *list, *dlist, *Mlist, *dMlist;
    718   StatType stats;
     713int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     714
     715  off_t i, N, Nmax;
    719716  Image *image;
    720717
    … …  
    722719  if (FREEZE_MOSAICS) return (FALSE);
    723720
     721  if (NTHREADS) {
     722    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
     723    return status;
     724  }
     725
    724726  image = getimages (&N, NULL);
    725727
    726728  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
    727 
    728   int Nsecfilt = GetPhotcodeNsecfilt ();
    729729
    730730  if (PoorImages) {
    … …  
    738738    Nmax = MAX (Nmax, N_onMosaic[i]);
    739739  }
    740   ALLOCATE (list, double, Nmax);
    741   ALLOCATE (dlist, double, Nmax);
    742   ALLOCATE (Mlist, double, Nmax);
    743   ALLOCATE (dMlist, double, Nmax);
    744 
    745   Nfew = Nbad = Ncal = Nrel = Ngrid = Nsys = 0;
    746 
    747   for (i = 0; i < Nmosaic; i++) {
    748    
    749     /* on PoorImages run, skip good images */
    750     if (PoorImages) {
    751       bad = mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
    752       if (!bad) continue;
    753     }     
    754 
    755     // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
    756     // we optionally do not recalibrate images with UBERCAL zero points
    757     if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
    758       mosaic[i].ubercalDist = 0;
    759       // propagate ubercalDist to the images
    760       for (j = 0; j < MosaicN_Image[i]; j++) {
    761         off_t im = MosaicToImage[i][j];
    762         image[im].ubercalDist = mosaic[i].ubercalDist;
    763         // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    764       }
    765       if (KEEP_UBERCAL) continue;
    766     }
    767 
    768     // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE)
    769     if (mosaic[i].skipCal) continue;
    770 
    771     int minUbercalDist = 1000;
    772 
    773     int testImage = FALSE;
    774     testImage |= (abs(mosaic[i].start - 1323003245) < 10);
    775     testImage |= (abs(mosaic[i].start - 1323003069) < 10);
    776     testImage |= (abs(mosaic[i].start - 1323003125) < 10);
    777     testImage |= (abs(mosaic[i].start - 1323003300) < 10);
    778     testImage |= (abs(mosaic[i].start - 1323003365) < 10);
    779     testImage |= (abs(mosaic[i].start - 1323003191) < 10);
    780     testImage |= (abs(mosaic[i].start - 1323003014) < 10);
    781     testImage |= (abs(mosaic[i].start - 1323003484) < 10);
    782     testImage |= (abs(mosaic[i].start - 1323003419) < 10);
    783     testImage |= (abs(mosaic[i].start - 1323002949) < 10);
    784 
    785     FILE *fout = NULL;
    786     if (testImage) {
    787       char filename[64];
    788       snprintf (filename, 64, "test.%05d.%02d.dat", (int) i, npass_output);
    789       fout = fopen (filename, "w");
    790     }
    791 
    792     // number of stars to measure the bright-end scatter
    793     Nbright = 0;
    794 
    795     N = 0;
    796     for (j = 0; j < N_onMosaic[i]; j++) {
    797      
    798       m = MosaicToMeasure[i][j];
    799       c = MosaicToCatalog[i][j];
    800      
    801       if (fout) {
    802         Mcal  = getMcal  (m, c, flatcorr, catalog);
    803         Mgrid = getMgrid (m, c);
    804         Mrel  = getMrel  (catalog, m, c);
    805         Mflat = getMflat (m, c, flatcorr, catalog);
    806 
    807         n = catalog[c].measureT[m].averef;
    808         Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
    809 
    810         float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
    811 
    812         int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
    813 
    814         fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad);
    815       }
    816 
    817       if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    818           Nbad ++;
    819           continue;
    820       }
    821       Mcal  = getMcal  (m, c, flatcorr, catalog);
    822       if (isnan(Mcal)) {
    823           Ncal++;
    824           continue;
    825       }
    826       Mgrid = getMgrid (m, c);
    827       if (isnan(Mgrid)) {
    828           Ngrid ++;
    829           continue;
    830       }
    831       Mrel  = getMrel  (catalog, m, c);
    832       if (isnan(Mrel)) {
    833           Nrel ++;
    834           continue;
    835       }
    836      
    837       // image.Mcal is not supposed to include the flat-field correction, so we need to
    838       // apply that offset as well here for this image (in other words, each detection is
    839       // being compared to the model, excluding the zero point, Mcal.  The model includes
    840       // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
    841 
    842       Mflat = getMflat (m, c, flatcorr, catalog);
    843 
    844       n = catalog[c].measureT[m].averef;
    845       Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
    846       if (isnan(Msys)) {
    847         Nsys++;
    848         continue;
    849       }
    850 
    851       PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
    852       if (!code) goto skip;
    853       if (code->equiv < 1) goto skip;
    854       int Nsec = GetPhotcodeNsec (code->equiv);
    855       if (Nsec == -1) goto skip;
    856       minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
    857 
    858     skip:
    859       list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
    860       dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    861       if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
    862         Mlist[Nbright] = list[N];
    863         dMlist[Nbright] = dlist[N];
    864         Nbright ++;
    865       }
    866       N++;
    867     }
    868     /* N_onMosaic[i] is all measurements, N is good measurements */
    869 
    870     if (fout) {
    871       fclose (fout);
    872     }
    873 
    874     /* too few good measurements or too many bad measurements (skip in PoorImages run) */
    875     if (!PoorImages) {
    876       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[i]);
    877       if (mark) {
    878         if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[MosaicToImage[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[i]); }
    879         mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
    880         Nfew ++;
    881         if (testImage) {
    882           fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[i]));
    883         }
    884       } else {
    885         mosaic[i].flags &= ~ID_IMAGE_PHOTOM_FEW;
    886       }
    887     }
    888     liststats (list, dlist, N, &stats);
    889     if (VERBOSE2 && PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
    890 
    891     mosaic[i].Mcal  = stats.mean;
    892     mosaic[i].dMcal = stats.error;
    893     mosaic[i].nFitPhotom = N;
    894     mosaic[i].Xm    = 100.0*log10(stats.chisq);
    895 
    896     if (testImage) {
    897       fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) i, mosaic[i].start, stats.mean, stats.error, mosaic[i].nFitPhotom);
    898     }
    899 
    900     plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
    901 
    902     // bright end scatter
    903     liststats (Mlist, dMlist, Nbright, &stats);
    904     mosaic[i].dMsys = stats.sigma;
    905 
    906     if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
    907       mosaic[i].Mcal = 0.0;
    908     }
    909 
    910     if (testImage) {
    911       fprintf (stderr, "%f %f\n", mosaic[i].Mcal, mosaic[i].dMsys);
    912     }
    913 
    914     // minUbercalDist calculated here is the min value for any star owned by this image
    915     // since this particular image is tied to that star, bump its distance by 1
    916     mosaic[i].ubercalDist = minUbercalDist + 1;
    917 
    918     // propagate ubercalDist to the images
    919     for (j = 0; j < MosaicN_Image[i]; j++) {
    920       off_t im = MosaicToImage[i][j];
    921       image[im].ubercalDist = mosaic[i].ubercalDist;
    922       // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    923     }
    924   }
    925   free (list);
    926   free (dlist);
    927   free (Mlist);
    928   free (dMlist);
    929 
    930   npass_output ++;
    931 
    932   fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Ngrid, Nrel, Nsys);
    933 
    934   if (PoorImages) {
    935     IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    936     STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    937     MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
    938   }
    939   return (TRUE);
    940 }
    941 # endif
    942  
    943 int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
    944 
    945   off_t i, N, Nmax;
    946   Image *image;
    947 
    948   if (!MOSAIC_ZEROPT) return (FALSE);
    949   if (FREEZE_MOSAICS) return (FALSE);
    950 
    951   if (NTHREADS) {
    952     int status = setMmos_threaded (catalog, PoorImages, flatcorr);
    953     return status;
    954   }
    955 
    956   image = getimages (&N, NULL);
    957 
    958   fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
    959 
    960   if (PoorImages) {
    961     // XXX use bad stars and measurements for PoorImages? or not?
    962     // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    963     IMAGE_BAD = 0;
    964   }
    965 
    966   Nmax = 0;
    967   for (i = 0; i < Nmosaic; i++) {
    968     Nmax = MAX (Nmax, N_onMosaic[i]);
    969   }
    970740
    971741  SetMmosInfo info;
    … …  
    973743
    974744  for (i = 0; i < Nmosaic; i++) {
    975     setMmos_mosaic (mosaic, i, image, catalog, &info, flatcorr);
     745    setMmos_mosaic (&mosaic[i], i, image, catalog, &info, flatcorr);
    976746  }
    977747  SetMmosInfoFree (&info);
    … …  
    1029799
    1030800  int testImage = FALSE;
     801  // testImage |= (abs(mosaic[0].start - 1324104046) < 10);
     802  // testImage |= (abs(mosaic[0].start - 1324103823) < 10);
    1031803  // testImage |= (abs(mosaic[0].start - 1323003245) < 10);
    1032804  // testImage |= (abs(mosaic[0].start - 1323003069) < 10);
    … …  
    1075847    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    1076848      info->Nbad ++;
    1077       return TRUE;
     849      continue;
    1078850    }
    1079851    Mcal  = getMcal  (m, c, flatcorr, catalog);
    1080852    if (isnan(Mcal)) {
    1081853      info->Ncal++;
    1082       return TRUE;
     854      continue;
    1083855    }
    1084856    Mgrid = getMgrid (m, c);
    1085857    if (isnan(Mgrid)) {
    1086858      info->Ngrid ++;
    1087       return TRUE;
     859      continue;
    1088860    }
    1089861    Mrel  = getMrel  (catalog, m, c);
    1090862    if (isnan(Mrel)) {
    1091863      info->Nrel ++;
    1092       return TRUE;
     864      continue;
    1093865    }
    1094866     
    … …  
    1104876    if (isnan(Msys)) {
    1105877      info->Nsys++;
    1106       return TRUE;
     878      continue;
    1107879    }
    1108880
    … …  
    1167939
    1168940  if (testImage) {
    1169     fprintf (stderr, "%f %f\n", mosaic[0].Mcal, mosaic[0].dMsys);
     941    fprintf (stderr, "%f %f  :  %d %f\n", mosaic[0].Mcal, mosaic[0].dMsys, mosaic[0].Xm, pow(10.0, 0.01*mosaic[0].Xm));
    1170942  }
    1171943
    … …  
    1215987  ThreadInfo *threadinfo;
    1216988  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     989
     990  // each time this function is called, we cycle through the available mosaics
     991  // make sure we start at 0
     992  nextMosaic = 0;;
    1217993
    1218994  // launch N worker threads
    … …  
    12661042  free (threadinfo);
    12671043
     1044  npass_output ++;
     1045
     1046  if (PoorImages) {
     1047    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
     1048    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     1049    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
     1050  }
    12681051  return TRUE;
    12691052}
    … …  
    12881071    Image *image = threadinfo->image;
    12891072
    1290     setMmos_mosaic (mosaic, i, image, catalog, &results, flatcorr);
     1073    setMmos_mosaic (&mosaic[i], i, image, catalog, &results, flatcorr);
    12911074    SetMmosInfoAccum (&threadinfo->info, &results);
    12921075  }
    … …  
    15071290  for (i = 0; i < Nmosaic; i++) {
    15081291    if (mosaic[i].flags & IMAGE_BAD) continue;
     1292    if (mosaic[i].skipCal) continue;
    15091293    list[n] = mosaic[i].Mcal;
    15101294    dlist[n] = 1;
    … …  
    15331317  n = 0;
    15341318  for (i = 0; i < Nmosaic; i++) {
    1535 
    15361319    if (mosaic[i].flags & IMAGE_BAD) continue;
     1320    if (mosaic[i].skipCal) continue;
    15371321    list[n] = mosaic[i].dMcal;
    15381322    dlist[n] = 1;
    … …  
    15631347  for (i = 0; i < Nmosaic; i++) {
    15641348    if (mosaic[i].flags & IMAGE_BAD)  continue;
     1349    if (mosaic[i].skipCal) continue;
    15651350
    15661351    N = 0;
    … …  
    16051390  n = 0;
    16061391  for (i = 0; i < Nmosaic; i++) {
    1607 
    16081392    if (mosaic[i].flags & IMAGE_BAD) continue;
     1393    if (mosaic[i].skipCal) continue;
    16091394    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
    16101395    dlist[n] = 1;
    … …  
    16371422  for (i = N = 0; i < Nmosaic; i++) {
    16381423    if (mosaic[i].flags & IMAGE_BAD) continue;
     1424    if (mosaic[i].skipCal) continue;
    16391425    mlist[N] = mosaic[i].Mcal;
    16401426    slist[N] = mosaic[i].dMcal;
    … …  
    16541440    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
    16551441    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
     1442
     1443    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
     1444    if (mosaic[i].skipCal) continue;
    16561445
    16571446    mark = FALSE;

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

    r33823 r33827  
    144144    SetMrelInfoAccum (&summary, &results);
    145145  }
    146   fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     146  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
    147147
    148148  SetMrelInfoFree (&results);
    … …  
    167167    SetMrelInfoAccum (&summary, &results);
    168168  }
    169   fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     169  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
    170170
    171171  SetMrelInfoFree (&results);
    … …  
    191191  ThreadInfo *threadinfo;
    192192  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
     193
     194  // each time this function is called, we cycle through the available catalogs.
     195  // make sure we start at 0
     196  nextCatalog = 0;;
    193197
    194198  // launch N worker threads
    … …  
    221225  // report stats & summary from the threads
    222226  for (i = 0; i < NTHREADS; i++) {
    223     fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n",
     227    if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n",
    224228             i,
    225229             threadinfo[i].summary.Ncode,
    … …  
    232236    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
    233237  }
    234   fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     238  if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
    235239  free (threadinfo);
    236240

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot.c

    r33797 r33827  
    128128      setMmos  (catalog, FALSE, flatcorr);
    129129      setMgrid (catalog, flatcorr);
     130      MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    130131   
    131132      if (PLOTSTUFF) {
    … …  
    139140      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    140141      // if (i % 6 == 1) rationalize_images ();
    141       if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
    142       if (i % 6 == 3) clean_stars (catalog, Ncatalog);
    143       if (i % 6 == 5) clean_mosaics ();
    144       if (i % 6 == 5) clean_images ();
     142
     143      // NOTE : in the past, I was not iterating enough before cleaning.  make sure we do
     144      // at least 8 loops first -- that should get the systematic errors down to the ~1%
     145      // level, even in cases where we have an even split between photometric data and
     146      // data with 1 mag of extinction.
     147
     148      if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
     149      if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
     150      if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
     151      if ((i > 8) && (i % 8 == 5)) clean_images ();
    145152
    146153      // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
    … …  
    148155      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
    149156      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
    150       global_stats (catalog, Ncatalog, flatcorr);
     157      if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr);
    151158      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    152159    }

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/select_images.c

    r33807 r33827  
    250250
    251251    // the sky region RA is defined to be 0 - 360.0
    252     if (RminImage > RmaxBand[iDecBandMin]) {
    253       fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     252    if (RminImage < RminBand[iDecBandMin]) {
     253      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    254254      continue;
    255255    }
    256     if (RminImage > RmaxBand[iDecBandMax]) {
    257       fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     256    if (RminImage < RminBand[iDecBandMax]) {
     257      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    258258      continue;
    259259    }
    260     if (RmaxImage < RminBand[iDecBandMin]) {
    261       fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     260    if (RmaxImage > RmaxBand[iDecBandMin]) {
     261      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    262262      continue;
    263263    }
    264     if (RmaxImage < RminBand[iDecBandMax]) {
    265       fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
     264    if (RmaxImage > RmaxBand[iDecBandMax]) {
     265      if (VERBOSE2) fprintf (stderr, "skip image %s (%f,%f) on boundary\n", timage[i].name, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    266266      continue;
    267267    }
     268
     269    // fprintf (ftest, "%s %lf %lf - %lf %lf : %lf %lf\n", timage[i].name, RminImage, DminImage, RmaxImage, DmaxImage, 0.5*(RminImage + RmaxImage), 0.5*(DminImage + DmaxImage));
    268270
    269271    // image overlaps region, keep it
Note: See TracChangeset for help on using the changeset viewer.