Index: /branches/eam_branches/ipp-20111122/Ohana/src/addstar/test/relphot.reject.dvo =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/addstar/test/relphot.reject.dvo (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/addstar/test/relphot.reject.dvo (revision 33391) @@ -34,4 +34,6 @@ $CHIP_DY = 1000 $DEFECT_FRAC = 0.03 + $OFFSET_FRAC_UC = 0.1 + $OFFSET_FRAC_NC = 0.5 if (not($?VERBOSE)) set VERBOSE = 0 @@ -54,5 +56,5 @@ # sequence to count images (only used in this function) # XXX TEST create tmpseq 0 9 - create tmpseq 0 3 + create tmpseq 0 9 ##### mjd and zpt values for ubercal'ed images @@ -63,10 +65,10 @@ $filt_uc:1 = g; filtN_uc[1] = 1 $filt_uc:2 = g; filtN_uc[2] = 1 - # $filt_uc:3 = r; filtN_uc[3] = 2 - # $filt_uc:4 = r; filtN_uc[4] = 2 - # $filt_uc:5 = r; filtN_uc[5] = 2 - # $filt_uc:6 = i; filtN_uc[6] = 3 - # $filt_uc:7 = i; filtN_uc[7] = 3 - # $filt_uc:8 = i; filtN_uc[8] = 3 + $filt_uc:3 = g; filtN_uc[3] = 1 + $filt_uc:4 = g; filtN_uc[4] = 1 + $filt_uc:5 = g; filtN_uc[5] = 1 + $filt_uc:6 = g; filtN_uc[6] = 1 + $filt_uc:7 = g; filtN_uc[7] = 1 + $filt_uc:8 = g; filtN_uc[8] = 1 # airmass slopes for this sequnece @@ -87,10 +89,10 @@ mjd_uc[1] = 55000.02 mjd_uc[2] = 55000.03 - # mjd_uc[3] = 55015.01 - # mjd_uc[4] = 55015.02 - # mjd_uc[5] = 55015.03 - # mjd_uc[6] = 55025.01 - # mjd_uc[7] = 55025.02 - # mjd_uc[8] = 55025.03 + mjd_uc[3] = 55015.01 + mjd_uc[4] = 55015.02 + mjd_uc[5] = 55015.03 + mjd_uc[6] = 55025.01 + mjd_uc[7] = 55025.02 + mjd_uc[8] = 55025.03 ##### mjd and zpt values for not-ubercal'ed images @@ -101,10 +103,10 @@ $filt_nc:1 = g; filtN_nc[1] = 1 $filt_nc:2 = g; filtN_nc[2] = 1 - # $filt_nc:3 = r; filtN_nc[3] = 2 - # $filt_nc:4 = r; filtN_nc[4] = 2 - # $filt_nc:5 = r; filtN_nc[5] = 2 - # $filt_nc:6 = i; filtN_nc[6] = 3 - # $filt_nc:7 = i; filtN_nc[7] = 3 - # $filt_nc:8 = i; filtN_nc[8] = 3 + $filt_nc:3 = g; filtN_nc[3] = 1 + $filt_nc:4 = g; filtN_nc[4] = 1 + $filt_nc:5 = g; filtN_nc[5] = 1 + $filt_nc:6 = g; filtN_nc[6] = 1 + $filt_nc:7 = g; filtN_nc[7] = 1 + $filt_nc:8 = g; filtN_nc[8] = 1 # airmass slopes for these sequnece @@ -126,10 +128,10 @@ mjd_nc[1] = 55000.12 mjd_nc[2] = 55000.13 - # mjd_nc[3] = 55015.11 - # mjd_nc[4] = 55015.12 - # mjd_nc[5] = 55015.13 - # mjd_nc[6] = 55025.11 - # mjd_nc[7] = 55025.12 - # mjd_nc[8] = 55025.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 end @@ -152,10 +154,10 @@ init mkzptfile - mkstars 500 + mkstars 10000 tapPLAN {(((mjd_uc[] + mjd_nc[])*4 + (mjd_uc[] + mjd_nc[])) + (mjd_uc[] + mjd_nc[]))} # generate the basic images and check they were correctly ingested by dvo - mksequence $fileroot $catdir + mksequence.drift $fileroot $catdir for i 0 mjd_uc[] ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw @@ -175,6 +177,7 @@ # run relphot on the serial db and check that the images now match the expected values - 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 - 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 + # 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 + 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 + 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 for i 0 mjd_uc[] ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot @@ -315,5 +318,31 @@ end -macro mksequence +macro mksequence.drift + if ($0 != 3) + echo "mksequence.drift (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + 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 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + 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 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mksequence.1pt if ($0 != 3) echo "mksequence (fileroot) (catdir)" @@ -495,5 +524,5 @@ # size of region of interest in linear arcseconds $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE - $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE * (1 + tmpseq[]*$OFFSET_FRAC_NC) create tmp 0 $1 Index: /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/libdvo/include/dvo.h (revision 33391) @@ -293,4 +293,5 @@ unsigned int imageID; unsigned int dbFlags; + unsigned int photFlags; int catID; // unsigned int? unsigned short photcode; Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/BrightCatalog.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/BrightCatalog.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/BrightCatalog.c (revision 33391) @@ -62,19 +62,20 @@ // need to create and assign to flat-field correction - GET_COLUMN(dR, "RA_OFF", float); - GET_COLUMN(dD, "DEC_OFF", float); - GET_COLUMN(M, "MAG_SYS", float); - GET_COLUMN(Mcal, "MAG_CAL", float); - GET_COLUMN(dM, "MAG_ERR", float); - GET_COLUMN(airmass, "AIRMASS", float); - GET_COLUMN(Xccd, "X_CCD", float); - GET_COLUMN(Yccd, "Y_CCD", float); - GET_COLUMN(dt, "EXPTIME", float); - GET_COLUMN(t, "TIME", int); - GET_COLUMN(averef, "AVE_REF", int); // XXX signed vs unsigned? - GET_COLUMN(imageID, "IMAGE_ID", int); - GET_COLUMN(dbFlags, "DB_FLAGS", int); - GET_COLUMN(catID, "CAT_ID", int); - GET_COLUMN(photcode, "PHOTCODE", short); + GET_COLUMN(dR, "RA_OFF", float); + GET_COLUMN(dD, "DEC_OFF", float); + GET_COLUMN(M, "MAG_SYS", float); + GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(airmass, "AIRMASS", float); + GET_COLUMN(Xccd, "X_CCD", float); + GET_COLUMN(Yccd, "Y_CCD", float); + GET_COLUMN(dt, "EXPTIME", float); + GET_COLUMN(t, "TIME", int); + GET_COLUMN(averef, "AVE_REF", int); // XXX signed vs unsigned? + GET_COLUMN(imageID, "IMAGE_ID", int); + GET_COLUMN(dbFlags, "DB_FLAGS", int); + GET_COLUMN(photFlags, "PHOT_FLAGS", int); + GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(photcode, "PHOTCODE", short); // XXX free the fits table data here @@ -82,19 +83,20 @@ ALLOCATE (measure, MeasureTiny, Nrow); for (i = 0; i < Nrow; i++) { - measure[i].dR = dR[i]; - measure[i].dD = dD[i]; - measure[i].M = M[i]; - measure[i].Mcal = Mcal[i]; - measure[i].dM = dM[i]; - measure[i].airmass = airmass[i]; - measure[i].Xccd = Xccd[i]; - measure[i].Yccd = Yccd[i]; - measure[i].dt = dt[i]; - measure[i].t = t[i]; - measure[i].averef = averef[i]; - measure[i].imageID = imageID[i]; - measure[i].dbFlags = dbFlags[i]; - measure[i].catID = catID[i]; - measure[i].photcode = photcode[i]; + measure[i].dR = dR[i]; + measure[i].dD = dD[i]; + measure[i].M = M[i]; + measure[i].Mcal = Mcal[i]; + measure[i].dM = dM[i]; + measure[i].airmass = airmass[i]; + measure[i].Xccd = Xccd[i]; + measure[i].Yccd = Yccd[i]; + measure[i].dt = dt[i]; + measure[i].t = t[i]; + measure[i].averef = averef[i]; + measure[i].imageID = imageID[i]; + measure[i].dbFlags = dbFlags[i]; + measure[i].photFlags = photFlags[i]; + measure[i].catID = catID[i]; + measure[i].photcode = photcode[i]; } // fprintf (stderr, "loaded data for %lld measures\n", (long long) Nrow); @@ -113,4 +115,5 @@ free (imageID ); free (dbFlags ); + free (photFlags ); free (catID ); free (photcode); @@ -261,19 +264,20 @@ // create intermediate storage arrays - float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); - float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); - float *M ; ALLOCATE (M , float, catalog->Nmeasure); - float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); - float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); - float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); - float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); - float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); - float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); - int *t ; ALLOCATE (t , int , catalog->Nmeasure); - int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); - int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); - int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); - int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); - short *photcode ; ALLOCATE (photcode, short, catalog->Nmeasure); + float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); + float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); + float *M ; ALLOCATE (M , float, catalog->Nmeasure); + float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); + float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); + float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); + float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); + float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); + int *t ; ALLOCATE (t , int , catalog->Nmeasure); + int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); + int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); + int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); + int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); + int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); + short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); // assign the storage arrays @@ -294,23 +298,25 @@ imageID[i] = measure[i].imageID ; dbFlags[i] = measure[i].dbFlags ; + photFlags[i]= measure[i].photFlags; photcode[i] = measure[i].photcode ; } // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "X_CCD", Xccd, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "Y_CCD", Yccd, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "EXPTIME", dt, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "TIME", t, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "AVE_REF", averef, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Nmeasure); - gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "X_CCD", Xccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "Y_CCD", Yccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", dt, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "TIME", t, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AVE_REF", averef, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); free (dR ); @@ -327,4 +333,5 @@ free (imageID ); free (dbFlags ); + free (photFlags ); free (catID ); free (photcode); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/ImageOps.c (revision 33391) @@ -456,8 +456,14 @@ } - minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt].M_80, minUbercalDist); // M_80 == ubercalDist + PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode); + if (!code) goto skip; + if (code->equiv < 1) goto skip; + int Nsec = GetPhotcodeNsec (code->equiv); + if (Nsec == -1) goto skip; + minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].M_80, minUbercalDist); // M_80 == ubercalDist // 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); + skip: list[N] = Msys - Mrel - Mmos - Mgrid + Mflat; dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR); @@ -497,4 +503,5 @@ image[i].dummy3 = minUbercalDist; + fprintf (stderr, "%s : %d\n", image[i].name, image[i].dummy3); } free (list); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/StarOps.c (revision 33391) @@ -56,4 +56,6 @@ int Nsecfilt = GetPhotcodeNsecfilt (); Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0; + + FILE *f = fopen ("few.v1.dat", "w"); for (i = 0; i < Ncatalog; i++) { @@ -135,4 +137,5 @@ // when performing the grid analysis, STAR_TOOFEW will be set to 1; if (N <= STAR_TOOFEW) { /* too few measurements */ + 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); catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW; Nfew ++; @@ -151,4 +154,5 @@ } } + fclose (f); 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); @@ -560,4 +564,6 @@ TOOFEW = MAX (5, STAR_TOOFEW); + FILE *f = fopen ("poor.dat", "w"); + Ndel = Nave = 0; Ncal = Nmos = Ngrid = Nfew = 0; @@ -598,5 +604,8 @@ N++; } - if (N <= TOOFEW) { Nfew ++; continue; } + if (N <= TOOFEW) { + Nfew ++; + continue; + } /* 3-sigma clip based on stats of inner 50% */ @@ -648,5 +657,10 @@ /* mark bad measures (> 3 sigma deviant) */ for (k = 0; k < N; k++) { - if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) { + // treat the scatter of the star as a systematic term. this is a bit of an + // over-estimage (a perfect Gauss distribution with perfect errors would have + // mySigma = sqrt(2) too large) + float mySigma = hypot (stats.sigma, dlist[k]); + if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) { + 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); catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM; if (final) { @@ -655,4 +669,8 @@ } Ndel ++; + } else { + if (catalog[i].measureT[ilist[k]].photFlags & 0x02) { + 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); + } } } @@ -661,4 +679,5 @@ } } + fclose (f); initstats (STATMODE); if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave); @@ -667,4 +686,6 @@ free (ilist); free (tlist); + + } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/bcatalog.c (revision 33391) @@ -18,18 +18,19 @@ int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) { - measureT[0].dR = measure[0].dR; - measureT[0].dD = measure[0].dD; - measureT[0].M = measure[0].M; - measureT[0].Mcal = measure[0].Mcal; - measureT[0].dM = measure[0].dM; - measureT[0].airmass = measure[0].airmass; - measureT[0].Xccd = measure[0].Xccd; - measureT[0].Yccd = measure[0].Yccd; - measureT[0].t = measure[0].t; - measureT[0].dt = measure[0].dt; - measureT[0].averef = measure[0].averef; - measureT[0].imageID = measure[0].imageID; - measureT[0].dbFlags = measure[0].dbFlags; - measureT[0].photcode = measure[0].photcode; + measureT[0].dR = measure[0].dR; + measureT[0].dD = measure[0].dD; + measureT[0].M = measure[0].M; + measureT[0].Mcal = measure[0].Mcal; + measureT[0].dM = measure[0].dM; + measureT[0].airmass = measure[0].airmass; + measureT[0].Xccd = measure[0].Xccd; + measureT[0].Yccd = measure[0].Yccd; + measureT[0].t = measure[0].t; + measureT[0].dt = measure[0].dt; + measureT[0].averef = measure[0].averef; + measureT[0].imageID = measure[0].imageID; + measureT[0].dbFlags = measure[0].dbFlags; + measureT[0].photFlags = measure[0].photFlags; + measureT[0].photcode = measure[0].photcode; return (TRUE); @@ -39,8 +40,8 @@ off_t i, j, offset; - int ecode, found, Ns; + int found, Ns; off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; float mag; - int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan; + int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor; int Nsecfilt = GetPhotcodeNsecfilt (); @@ -56,5 +57,5 @@ Nmeasure = Naverage = 0; - Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = 0; + Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0; // copy the following fields to the subcatalog: @@ -105,8 +106,9 @@ /* select measurements by photcode */ - ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode); + PhotCode *code = GetPhotcodebyCode (catalog[0].measure[offset].photcode); + if (!code) continue; found = FALSE; for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) { - if (ecode == photcodes[Ns][0].code) found = TRUE; + if (code->equiv == photcodes[Ns][0].code) found = TRUE; } if (!found) { @@ -127,4 +129,7 @@ if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; } if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; } + + if (catalog[0].measure[offset].photFlags & code->photomBadMask) { Nbad++; continue;} + if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;} // check for galaxies @@ -206,6 +211,6 @@ fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", subcatalog[0].Naverage, subcatalog[0].Nmeasure, i, catalog[0].filename); - fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies\n", - Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy); + 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", + Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor); } return (TRUE); @@ -336,18 +341,19 @@ for (i = 0; i < catalog[0].Nmeasure; i++) { // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]); - measureT[i].dR = measure[i].dR; - measureT[i].dD = measure[i].dD; - measureT[i].M = measure[i].M; - measureT[i].Mcal = measure[i].Mcal; - measureT[i].dM = measure[i].dM; - measureT[i].airmass = measure[i].airmass; - measureT[i].Xccd = measure[i].Xccd; - measureT[i].Yccd = measure[i].Yccd; - measureT[i].t = measure[i].t; - measureT[i].dt = measure[i].dt; - measureT[i].averef = measure[i].averef; - measureT[i].imageID = measure[i].imageID; - measureT[i].dbFlags = measure[i].dbFlags; - measureT[i].photcode = measure[i].photcode; + measureT[i].dR = measure[i].dR; + measureT[i].dD = measure[i].dD; + measureT[i].M = measure[i].M; + measureT[i].Mcal = measure[i].Mcal; + measureT[i].dM = measure[i].dM; + measureT[i].airmass = measure[i].airmass; + measureT[i].Xccd = measure[i].Xccd; + measureT[i].Yccd = measure[i].Yccd; + measureT[i].t = measure[i].t; + measureT[i].dt = measure[i].dt; + measureT[i].averef = measure[i].averef; + measureT[i].imageID = measure[i].imageID; + measureT[i].dbFlags = measure[i].dbFlags; + measureT[i].photFlags = measure[i].photFlags; + measureT[i].photcode = measure[i].photcode; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/liststats.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/liststats.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/liststats.c (revision 33391) @@ -106,5 +106,5 @@ dS += M; } - X2 = X2 / Nm; + X2 = X2 / (Nm - 1); dS = sqrt (dS / Nm); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c (revision 33390) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relphot/src/setMrelFinal.c (revision 33391) @@ -7,5 +7,5 @@ void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) { - off_t i, j, m; + off_t i; int ecode; @@ -36,13 +36,14 @@ for (i = 0; i < catalog[0].Naverage; i++) { - catalog[0].secfilt[Nsecfilt*i+Nsec].M = NAN; - catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN; - catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nsecfilt*j+Nsec].Ncode = 0; - catalog[0].secfilt[Nsecfilt*j+Nsec].Nused = 0; - - catalog[0].secfilt[Nsecfilt*i+Nsec].flags &= ~PHOTOM_PASSES; - - m = catalog[0].average[i].measureOffset; + off_t N = Nsecfilt*i+Nsec; + catalog[0].secfilt[N].M = NAN; + catalog[0].secfilt[N].dM = NAN; + catalog[0].secfilt[N].Xm = NAN_S_SHORT; + catalog[0].secfilt[N].Ncode = 0; + catalog[0].secfilt[N].Nused = 0; + catalog[0].secfilt[N].flags &= ~PHOTOM_PASSES; + + off_t m = catalog[0].average[i].measureOffset; + off_t j; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { @@ -89,9 +90,10 @@ catalog[0].average[i].flags &= ~ID_STAR_FEW; catalog[0].average[i].flags &= ~ID_STAR_POOR; + off_t j; for (j = 0; j < Nsecfilt; j++) { catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW; catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR; } - m = catalog[0].average[i].measureOffset; + off_t m = catalog[0].average[i].measureOffset; for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) { catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;