Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/Makefile =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/Makefile (revision 42694) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/Makefile (revision 42695) @@ -96,4 +96,5 @@ $(SRC)/find_matches_refstars.$(ARCH).o \ $(SRC)/find_matches_closest_refstars.$(ARCH).o \ +$(SRC)/photcodes.$(ARCH).o \ $(SRC)/getgsc.$(ARCH).o \ $(SRC)/getusno.$(ARCH).o \ @@ -151,4 +152,5 @@ $(SRC)/find_matches_refstars.$(ARCH).o \ $(SRC)/find_subset.$(ARCH).o \ +$(SRC)/photcodes.$(ARCH).o \ $(SRC)/getgsc.$(ARCH).o \ $(SRC)/getusno.$(ARCH).o \ @@ -181,4 +183,5 @@ $(SRC)/find_matches_refstars.$(ARCH).o \ $(SRC)/find_subset.$(ARCH).o \ +$(SRC)/photcodes.$(ARCH).o \ $(SRC)/getgsc.$(ARCH).o \ $(SRC)/getusno.$(ARCH).o \ Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/include/addstar.h =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/include/addstar.h (revision 42694) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/include/addstar.h (revision 42695) @@ -141,4 +141,5 @@ int OLD_RESORT; int READ_XRAD_DATA; +int CONVOLVED_APERTURES; int DIFF_WITH_INV; @@ -343,4 +344,7 @@ void GetConfig (char *config, char *field, char *format, int N, void *ptr); +int isGPC1warp (int photcode); +int isGPC1stack (int photcode); + /** there is an inconsistency to be resolved: fixed structures (like Image) Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/ReadXradFITS.c =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/ReadXradFITS.c (revision 42694) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/ReadXradFITS.c (revision 42695) @@ -97,80 +97,137 @@ // first entry and compare to the rest - float fwhmValues[3]; - fwhmValues[0] = PSFfwhm[0]; - fwhmValues[1] = PSFfwhm[1]; - fwhmValues[2] = PSFfwhm[2]; - - int i; int Nap = 0; - for (i = 0; i < catalog->Nmeasure; i++) { - dvo_lensing_init (&catalog->lensing[i]); - - if (catalog->measure[i].detID < RadID[Nap]) { - continue; + + if (CONVOLVED_APERTURES) { + // XXX code to support the 3 convolved apertures + + float fwhmValues[3]; + fwhmValues[0] = PSFfwhm[0]; + fwhmValues[1] = PSFfwhm[1]; + fwhmValues[2] = PSFfwhm[2]; + + for (int i = 0; i < catalog->Nmeasure; i++) { + dvo_lensing_init (&catalog->lensing[i]); + + if (catalog->measure[i].detID < RadID[Nap]) { + continue; + // this is a psf measurement which does not have a radial aperture + } + if (catalog->measure[i].detID > RadID[Nap]) { + myAbort("radial apertures for source not in psf list? sources out of order? seems like a bug\n"); + // this could be a radial aperture which does not have a PSF source, but that is not possible + } + + // confirm the 3 FWHM values: + myAssert (fwhmValues[0] == PSFfwhm[Nap+0], "FWHM mismatch %f vs %f", fwhmValues[0], PSFfwhm[Nap+0]); + myAssert (fwhmValues[1] == PSFfwhm[Nap+1], "FWHM mismatch %f vs %f", fwhmValues[1], PSFfwhm[Nap+1]); + myAssert (fwhmValues[2] == PSFfwhm[Nap+2], "FWHM mismatch %f vs %f", fwhmValues[2], PSFfwhm[Nap+2]); + + // EAM 2022.02.17 : here is the comment for the PV3 load: + // XXX this is all hard-wired and should make use of the headers. + // psphot cmfs have 5 radial apertures: + // array 0, 1, 2, 3, 4 + // SDSS 3, 4, 5, 6, 7 + + // EAM 2022.02.17 : here is the situation for UNIONS DR3: + // we have 3 convolutions (raw, 6", 8") + // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec + // I am going to save (4, 16, 32) which have index of (0, 2, 3) + +# define RAD_0 0 +# define RAD_1 2 +# define RAD_2 3 + catalog->lensing[i]. F_ApR5 = AperFlux [(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].dF_ApR5 = AperFluxErr[(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].sF_ApR5 = AperFluxStd[(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].fF_ApR5 = AperFill [(Nap + 0)*Ncol + RAD_0]; + + catalog->lensing[i]. F_ApR6 = AperFlux [(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].dF_ApR6 = AperFluxErr[(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].sF_ApR6 = AperFluxStd[(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].fF_ApR6 = AperFill [(Nap + 0)*Ncol + RAD_1]; + + catalog->lensing[i]. F_ApR7 = AperFlux [(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].dF_ApR7 = AperFluxErr[(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].sF_ApR7 = AperFluxStd[(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].fF_ApR7 = AperFill [(Nap + 0)*Ncol + RAD_2]; + + catalog->lensing[i]. F_ApR5_C1 = AperFlux [(Nap + 1)*Ncol + RAD_0]; + catalog->lensing[i].dF_ApR5_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_0]; + catalog->lensing[i]. F_ApR6_C1 = AperFlux [(Nap + 1)*Ncol + RAD_1]; + catalog->lensing[i].dF_ApR6_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_1]; + catalog->lensing[i]. F_ApR7_C1 = AperFlux [(Nap + 1)*Ncol + RAD_2]; + catalog->lensing[i].dF_ApR7_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_2]; + + catalog->lensing[i]. F_ApR5_C2 = AperFlux [(Nap + 2)*Ncol + RAD_0]; + catalog->lensing[i].dF_ApR5_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_0]; + catalog->lensing[i]. F_ApR6_C2 = AperFlux [(Nap + 2)*Ncol + RAD_1]; + catalog->lensing[i].dF_ApR6_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_1]; + catalog->lensing[i]. F_ApR7_C2 = AperFlux [(Nap + 2)*Ncol + RAD_2]; + catalog->lensing[i].dF_ApR7_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_2]; + + catalog->lensing[i].detID = catalog->measure[i].detID; + + // XXX set the measure, object, etc ID values here + // catalog->lensing[i].objID : set in find_matches_closest.c + // catalog->lensing[i].catID : set in find_matches_closest.c + // catalog->lensing[i].averef : set in find_matches_closest.c + // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c + + Nap += 3; + } + } else { + // PSF values for unconvolve data are undefined ? + + for (int i = 0; i < catalog->Nmeasure; i++) { + dvo_lensing_init (&catalog->lensing[i]); + // this is a psf measurement which does not have a radial aperture + if (catalog->measure[i].detID < RadID[Nap]) { continue; } + + // this could be a radial aperture which does not have a PSF source, but that is not possible + if (catalog->measure[i].detID > RadID[Nap]) { + myAbort("radial apertures for source not in psf list? sources out of order? seems like a bug\n"); + } + + // EAM 2022.02.17 : here is the comment for the PV3 load: + // XXX this is all hard-wired and should make use of the headers. + // psphot cmfs have 5 radial apertures: + // array 0, 1, 2, 3, 4 + // SDSS 3, 4, 5, 6, 7 + + // EAM 2022.02.17 : here is the situation for UNIONS DR4: + // we have 1 convolution (raw) + // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec + // I am going to save (4, 16, 32) which have index of (0, 2, 3) + +# define RAD_0 0 +# define RAD_1 2 +# define RAD_2 3 + catalog->lensing[i]. F_ApR5 = AperFlux [(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].dF_ApR5 = AperFluxErr[(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].sF_ApR5 = AperFluxStd[(Nap + 0)*Ncol + RAD_0]; + catalog->lensing[i].fF_ApR5 = AperFill [(Nap + 0)*Ncol + RAD_0]; + + catalog->lensing[i]. F_ApR6 = AperFlux [(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].dF_ApR6 = AperFluxErr[(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].sF_ApR6 = AperFluxStd[(Nap + 0)*Ncol + RAD_1]; + catalog->lensing[i].fF_ApR6 = AperFill [(Nap + 0)*Ncol + RAD_1]; + + catalog->lensing[i]. F_ApR7 = AperFlux [(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].dF_ApR7 = AperFluxErr[(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].sF_ApR7 = AperFluxStd[(Nap + 0)*Ncol + RAD_2]; + catalog->lensing[i].fF_ApR7 = AperFill [(Nap + 0)*Ncol + RAD_2]; + + catalog->lensing[i].detID = catalog->measure[i].detID; + + // XXX set the measure, object, etc ID values here + // catalog->lensing[i].objID : set in find_matches_closest.c + // catalog->lensing[i].catID : set in find_matches_closest.c + // catalog->lensing[i].averef : set in find_matches_closest.c + // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c + + Nap += 1; } - if (catalog->measure[i].detID > RadID[Nap]) { - myAbort("radial apertures for source not in psf list? sources out of order? seems like a bug\n"); - // this could be a radial aperture which does not have a PSF source, but that is not possible - } - - // confirm the 3 FWHM values: - myAssert (fwhmValues[0] == PSFfwhm[Nap+0], "FWHM mismatch %f vs %f", fwhmValues[0], PSFfwhm[Nap+0]); - myAssert (fwhmValues[1] == PSFfwhm[Nap+1], "FWHM mismatch %f vs %f", fwhmValues[1], PSFfwhm[Nap+1]); - myAssert (fwhmValues[2] == PSFfwhm[Nap+2], "FWHM mismatch %f vs %f", fwhmValues[2], PSFfwhm[Nap+2]); - - // EAM 2022.02.17 : here is the comment for the PV3 load: - // XXX this is all hard-wired and should make use of the headers. - // psphot cmfs have 5 radial apertures: - // array 0, 1, 2, 3, 4 - // SDSS 3, 4, 5, 6, 7 - - // EAM 2022.02.17 : here is the situation for UNIONS DR3: - // we have 3 convolutions (raw, 6", 8") - // for each we have 6 apertures with max radii of (4, 8, 16, 32, 48, 64) pixels = (1, 2, 4, 8, 12, 16) arcsec - // I am going to save (4, 16, 32) which have index of (0, 2, 3) - - # define RAD_0 0 - # define RAD_1 2 - # define RAD_2 3 - catalog->lensing[i]. F_ApR5 = AperFlux [(Nap + 0)*Ncol + RAD_0]; - catalog->lensing[i].dF_ApR5 = AperFluxErr[(Nap + 0)*Ncol + RAD_0]; - catalog->lensing[i].sF_ApR5 = AperFluxStd[(Nap + 0)*Ncol + RAD_0]; - catalog->lensing[i].fF_ApR5 = AperFill [(Nap + 0)*Ncol + RAD_0]; - - catalog->lensing[i]. F_ApR6 = AperFlux [(Nap + 0)*Ncol + RAD_1]; - catalog->lensing[i].dF_ApR6 = AperFluxErr[(Nap + 0)*Ncol + RAD_1]; - catalog->lensing[i].sF_ApR6 = AperFluxStd[(Nap + 0)*Ncol + RAD_1]; - catalog->lensing[i].fF_ApR6 = AperFill [(Nap + 0)*Ncol + RAD_1]; - - catalog->lensing[i]. F_ApR7 = AperFlux [(Nap + 0)*Ncol + RAD_2]; - catalog->lensing[i].dF_ApR7 = AperFluxErr[(Nap + 0)*Ncol + RAD_2]; - catalog->lensing[i].sF_ApR7 = AperFluxStd[(Nap + 0)*Ncol + RAD_2]; - catalog->lensing[i].fF_ApR7 = AperFill [(Nap + 0)*Ncol + RAD_2]; - - catalog->lensing[i]. F_ApR5_C1 = AperFlux [(Nap + 1)*Ncol + RAD_0]; - catalog->lensing[i].dF_ApR5_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_0]; - catalog->lensing[i]. F_ApR6_C1 = AperFlux [(Nap + 1)*Ncol + RAD_1]; - catalog->lensing[i].dF_ApR6_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_1]; - catalog->lensing[i]. F_ApR7_C1 = AperFlux [(Nap + 1)*Ncol + RAD_2]; - catalog->lensing[i].dF_ApR7_C1 = AperFluxErr[(Nap + 1)*Ncol + RAD_2]; - - catalog->lensing[i]. F_ApR5_C2 = AperFlux [(Nap + 2)*Ncol + RAD_0]; - catalog->lensing[i].dF_ApR5_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_0]; - catalog->lensing[i]. F_ApR6_C2 = AperFlux [(Nap + 2)*Ncol + RAD_1]; - catalog->lensing[i].dF_ApR6_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_1]; - catalog->lensing[i]. F_ApR7_C2 = AperFlux [(Nap + 2)*Ncol + RAD_2]; - catalog->lensing[i].dF_ApR7_C2 = AperFluxErr[(Nap + 2)*Ncol + RAD_2]; - - catalog->lensing[i].detID = catalog->measure[i].detID; - - // XXX set the measure, object, etc ID values here - // catalog->lensing[i].objID : set in find_matches_closest.c - // catalog->lensing[i].catID : set in find_matches_closest.c - // catalog->lensing[i].averef : set in find_matches_closest.c - // catalog->lensing[i].imageID : set in FilterStars.c and UpdateImageIDs.c - - Nap += 3; } myAssert (Nap == Nrow, "did we go too far???"); Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/args.c =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/args.c (revision 42694) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/args.c (revision 42695) @@ -91,4 +91,9 @@ remove_argument (N, &argc, argv); READ_XRAD_DATA = TRUE; + } + CONVOLVED_APERTURES = FALSE; + if ((N = get_argument (argc, argv, "-convolved-ap"))) { + remove_argument (N, &argc, argv); + CONVOLVED_APERTURES = TRUE; } DIFF_WITH_INV = FALSE; Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/find_matches_closest.c =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/find_matches_closest.c (revision 42694) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/find_matches_closest.c (revision 42695) @@ -248,8 +248,27 @@ /* in UPDATE mode, this value is not saved; use relphot to recalculate */ if (Nsec > -1) { + if (isGPC1warp(tgtcat[0].measure[Nmeas].photcode)) { + if (FALSE && isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfWrp)) { + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + } + goto set_average; + } + if (isGPC1stack(tgtcat[0].measure[Nmeas].photcode)) { + if (FALSE && isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfStk)) { + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + } + goto set_average; + } if (isnan(tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp)) { - tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); - } - } + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MkronChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + tgtcat[0].secfilt[n*Nsecfilt+Nsec].MapChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + } + } + set_average: /* Nm is updated, but not written out in -update mode (for existing entries) @@ -352,6 +371,23 @@ /* in UPDATE mode, this value is not saved; use relphot to recalculate */ if (Nsec > -1) { - tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); - } + if (isGPC1warp(tgtcat[0].measure[Nmeas].photcode)) { + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapWrp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + goto set_average_2; + } + if (isGPC1stack(tgtcat[0].measure[Nmeas].photcode)) { + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + // tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapStk = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + goto set_average_2; + } + if (isnan(tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp)) { + tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MkronChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_KRON); + tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MapChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_APER); + } + } + set_average_2: Nmeas ++; Index: /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/photcodes.c =================================================================== --- /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/photcodes.c (revision 42695) +++ /branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/photcodes.c (revision 42695) @@ -0,0 +1,37 @@ +# include "addstar.h" + +// for now (20140710) I need to identify gpc1 & gpc2 stacks explicitly. generalize in the future +int isGPC1stack (int photcode) { + + // 11000 : GPC1.g.SkyChip + // 11001 : VARIOUS.g.SkyChip (output of xcamera stacks) + // 14100 : SIMTEST.g.SkyChip + // 31000 : GPC2.g.SkyChip + + if ((photcode == 11000) || (photcode == 11001) || (photcode == 14100) || (photcode == 31000)) return TRUE; // g-band + if ((photcode == 11100) || (photcode == 11101) || (photcode == 14200) || (photcode == 31100)) return TRUE; // r-band + if ((photcode == 11200) || (photcode == 11201) || (photcode == 14300) || (photcode == 31200)) return TRUE; // i-band + if ((photcode == 11300) || (photcode == 11301) || (photcode == 14400) || (photcode == 31300)) return TRUE; // z-band + if ((photcode == 11400) || (photcode == 11401) || (photcode == 14500) || (photcode == 31400)) return TRUE; // y-band + if ((photcode == 11500) || (photcode == 11501) || (photcode == 14600) || (photcode == 31500)) return TRUE; // w-band + + return FALSE; +} + +// for now (20140710) I need to identify gpc1 & gpc2 warps explicitly. generalize in the future +int isGPC1warp (int photcode) { + + // 12000 : GPC1.g.ForcedWarp + // 15100 : SIMTEST.g.ForcedWarp + // 32000 : GPC2.g.ForcedWarp + + if ((photcode == 12000) || (photcode == 15100) || (photcode == 32000)) return TRUE; // g-band + if ((photcode == 12100) || (photcode == 15200) || (photcode == 32100)) return TRUE; // r-band + if ((photcode == 12200) || (photcode == 15300) || (photcode == 32200)) return TRUE; // i-band + if ((photcode == 12300) || (photcode == 15400) || (photcode == 32300)) return TRUE; // z-band + if ((photcode == 12400) || (photcode == 15500) || (photcode == 32400)) return TRUE; // y-band + if ((photcode == 12500) || (photcode == 15600) || (photcode == 32500)) return TRUE; // w-band + + return FALSE; +} +