Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 39670) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 39926) @@ -114,30 +114,33 @@ /* Measure.flags values -- these values are 32 bit (as of PS1_V1) */ typedef enum { - ID_MEAS_NOCAL = 0x00000001, // detection ignored for this analysis (photcode, time range) -- internal only - ID_MEAS_POOR_PHOTOM = 0x00000002, // detection is photometry outlier - ID_MEAS_SKIP_PHOTOM = 0x00000004, // detection was ignored for photometry measurement - ID_MEAS_AREA = 0x00000008, // detection near image edge - ID_MEAS_POOR_ASTROM = 0x00000010, // detection is astrometry outlier - ID_MEAS_SKIP_ASTROM = 0x00000020, // detection was ignored for astrometry measurement - ID_MEAS_USED_OBJ = 0x00000040, // detection was used during update objects - ID_MEAS_USED_CHIP = 0x00000080, // detection was used during update chips (XXX this probably does not make it into the db) - ID_MEAS_BLEND_MEAS = 0x00000100, // detection is within radius of multiple objects - ID_MEAS_BLEND_OBJ = 0x00000200, // multiple detections within radius of object - ID_MEAS_WARP_USED = 0x00000400, // measurement used to find mean warp photometry - ID_MEAS_UNMASKED_ASTRO = 0x00000800, // measurement was not masked in final astrometry fit - ID_MEAS_BLEND_MEAS_X = 0x00001000, // detection is within radius of multiple objects across catalogs - ID_MEAS_ARTIFACT = 0x00002000, // detection is thought to be non-astronomical - ID_MEAS_SYNTH_MAG = 0x00004000, // magnitude is synthetic - ID_MEAS_PHOTOM_UBERCAL = 0x00008000, // externally-supplied zero point from ubercal analysis - ID_MEAS_STACK_PRIMARY = 0x00010000, // this stack measurement is in the primary skycell - ID_MEAS_STACK_PHOT_SRC = 0x00020000, // this measurement supplied the stack photometry - ID_MEAS_ICRF_QSO = 0x00040000, // this measurement is an ICRF reference position - ID_MEAS_IMAGE_EPOCH = 0x00080000, // this measurement is registered to the image epoch (not tied to ref catalog epoch) - ID_MEAS_PHOTOM_PSF = 0x00100000, // this measurement is used for the mean psf mag - ID_MEAS_PHOTOM_APER = 0x00200000, // this measurement is used for the mean ap mag - ID_MEAS_PHOTOM_KRON = 0x00400000, // this measurement is used for the mean kron mag - ID_MEAS_MASKED_PSF = 0x01000000, // this measurement is masked based on IRLS weights for mean psf mag - ID_MEAS_MASKED_APER = 0x02000000, // this measurement is masked based on IRLS weights for mean ap mag - ID_MEAS_MASKED_KRON = 0x04000000, // this measurement is masked based on IRLS weights for mean kron mag + ID_MEAS_NOCAL = 0x00000001, // detection ignored for this analysis (photcode, time range) -- internal only + ID_MEAS_POOR_PHOTOM = 0x00000002, // detection is photometry outlier + ID_MEAS_SKIP_PHOTOM = 0x00000004, // detection was ignored for photometry measurement + ID_MEAS_AREA = 0x00000008, // detection near image edge + ID_MEAS_POOR_ASTROM = 0x00000010, // detection is astrometry outlier + ID_MEAS_SKIP_ASTROM = 0x00000020, // detection was ignored for astrometry measurement + ID_MEAS_USED_OBJ = 0x00000040, // detection was used during update objects + ID_MEAS_USED_CHIP = 0x00000080, // detection was used during update chips (XXX this probably does not make it into the db) + ID_MEAS_BLEND_MEAS = 0x00000100, // detection is within radius of multiple objects + ID_MEAS_BLEND_OBJ = 0x00000200, // multiple detections within radius of object + ID_MEAS_WARP_USED = 0x00000400, // measurement used to find mean warp photometry + ID_MEAS_UNMASKED_ASTRO = 0x00000800, // measurement was not masked in final astrometry fit + ID_MEAS_BLEND_MEAS_X = 0x00001000, // detection is within radius of multiple objects across catalogs + ID_MEAS_ARTIFACT = 0x00002000, // detection is thought to be non-astronomical + ID_MEAS_SYNTH_MAG = 0x00004000, // magnitude is synthetic + ID_MEAS_PHOTOM_UBERCAL = 0x00008000, // externally-supplied zero point from ubercal analysis + ID_MEAS_STACK_PRIMARY = 0x00010000, // this stack measurement is in the primary skycell + ID_MEAS_STACK_PHOT_SRC = 0x00020000, // this measurement supplied the stack photometry + ID_MEAS_ICRF_QSO = 0x00040000, // this measurement is an ICRF reference position + ID_MEAS_IMAGE_EPOCH = 0x00080000, // this measurement is registered to the image epoch (not tied to ref catalog epoch) + ID_MEAS_PHOTOM_PSF = 0x00100000, // this measurement is used for the mean psf mag + ID_MEAS_PHOTOM_APER = 0x00200000, // this measurement is used for the mean ap mag + ID_MEAS_PHOTOM_KRON = 0x00400000, // this measurement is used for the mean kron mag + ID_MEAS_MASKED_PSF = 0x01000000, // this measurement is masked based on IRLS weights for mean psf mag + ID_MEAS_MASKED_APER = 0x02000000, // this measurement is masked based on IRLS weights for mean ap mag + ID_MEAS_MASKED_KRON = 0x04000000, // this measurement is masked based on IRLS weights for mean kron mag + ID_MEAS_OBJECT_HAS_2MASS = 0x10000000, // measurement comes from an object with 2mass data + ID_MEAS_OBJECT_HAS_GAIA = 0x20000000, // measurement comes from an object with gaia data + ID_MEAS_OBJECT_HAS_TYCHO = 0x40000000, // measurement comes from an object with tycho data } DVOMeasureFlags; @@ -233,4 +236,10 @@ ID_SECF_STACK_BESTDET = 0x00008000, // PS1 stack best measurement is a detection (not forced) ID_SECF_STACK_PRIMDET = 0x00010000, // PS1 stack primary measurement is a detection (not forced) + + ID_SECF_HAS_SDSS = 0x00100000, // this photcode has SDSS photometry + ID_SECF_HAS_HSC = 0x00200000, // this photcode has HSC photometry + ID_SECF_HAS_CFH = 0x00400000, // this photcode has CFH photometry (mostly Megacam) + ID_SECF_HAS_DES = 0x00800000, // this photcode has DES photometry + ID_SECF_OBJ_EXT = 0x01000000, // extended in this band Index: trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c =================================================================== --- trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c (revision 39670) +++ trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c (revision 39926) @@ -59,4 +59,6 @@ } +// note that this function is called by the dvo clients and uses the 'metadata' table +// generated by the dvo shell for the mosaic images only. Coords *MatchMosaicMetadata (unsigned int imageID) { @@ -73,4 +75,6 @@ mosaic.crval2 = image[m].crval2; + // note that image->theta is calculated based on pc1_1, pc1_2 when the metadata file is + // generated mosaic.pc1_1 = cos(RAD_DEG*image[m].theta); mosaic.pc1_2 = sin(RAD_DEG*image[m].theta); Index: trunk/Ohana/src/libdvo/src/coordops.c =================================================================== --- trunk/Ohana/src/libdvo/src/coordops.c (revision 39670) +++ trunk/Ohana/src/libdvo/src/coordops.c (revision 39926) @@ -84,4 +84,5 @@ /** extra polynomial terms **/ + // for ZPN, these are used to modify the radial distance and not the X,Y coords if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) { X2 = X*X; @@ -126,5 +127,7 @@ if (proj == PROJ_WRP) { if (!coords->mosaic) { - myAbort ("missing mosaic element"); + // myAbort ("missing mosaic element"); + *ra = L; + *dec = M; return (FALSE); } @@ -157,8 +160,17 @@ ctht = 0.0; } else { - T = DEG_RAD / R; - stht = T / sqrt ( 1.0 + T*T); - ctht = 1.0 / sqrt ( 1.0 + T*T); + // T = DEG_RAD / R; // T in 1/radians + // stht = T / sqrt ( 1.0 + T*T); + // ctht = 1.0 / sqrt ( 1.0 + T*T); + + T = RAD_DEG * R; + stht = 1.0 / sqrt ( 1.0 + T*T); + ctht = T / sqrt ( 1.0 + T*T); } + break; + case PROJ_SIN: + // R = (180/pi) cos (theta) + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); break; case PROJ_STG: @@ -167,9 +179,4 @@ ctht = sqrt (1 - stht*stht); break; - case PROJ_SIN: - // R = (180/pi) cos (theta) - ctht = RAD_DEG * R; - stht = sqrt (1 - ctht*ctht); - break; case PROJ_ARC: // R = 90 - theta (degrees) @@ -177,36 +184,4 @@ stht = cos (RAD_DEG * R); break; - - case PROJ_ZPN: - // R = 90 - theta (degrees) - // this is wrong because we are ignoring the distortion - // XXX For now, just solve for terms up to n = 3 - - // Ro = (pi/180)(90 - theta) - // R = (180/pi)sum (P_i R^i) - - if (UKIRT_ONLY) { - double Ro = RAD_DEG * R; - double P1 = coords[0].polyterms[0][1]; - double P3 = coords[0].polyterms[0][3]; - - // find the roots of f(gamma) = P1 gamma + P3 gamma^3 - Ro - // starting guess for gamma is Ro / P1 - double gamma = Ro / P1; - - int i; - for (i = 0; i < 5; i++) { - double F = P1*gamma + P3*gamma*gamma*gamma - Ro; - double dFdgamma = P1 + 3.0*P3*gamma*gamma; - - double gamma_new = gamma - F / dFdgamma; - gamma = gamma_new; - } - - double theta = 90.0 - gamma * DEG_RAD ; - ctht = cos (RAD_DEG * theta); - stht = sin (RAD_DEG * theta); - break; - } case PROJ_ZEA: @@ -220,4 +195,51 @@ ctht = sqrt (1 - stht*stht); break; + + case PROJ_ZPN: + + // the forward projection is: + // theta = atan2(stht, ctht) + // gamma = (pi/2 - theta) : theta in radians + // Ro = sum (P_i gamma^i) + // R = (180/pi) Ro + + // given R, we need to find theta: + // Ro = R * (pi / 180) = sum (P_i gamma^i) + // solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson + + // use Ro to get a guess for gamma and iterate + + { + double Ro = RAD_DEG * R; + + // find the roots of f(gamma) - Ro = 0 + // starting guess for gamma is (Ro - P0) / P1 + double gamma = (Ro - coords[0].polyterms[0][0]) / coords[0].polyterms[1][0]; + + int iter; + for (iter = 0; iter < 5; iter++) { + + double Rc = 0.0; // this will hold the ander + double dR = 0.0; + for (int i = coords[0].Npolyterms - 1; i > 1; i--) { + double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1]; + Rc = (Rc + Pi)*gamma; + dR = (dR + i*Pi)*gamma; + } + double P0 = coords[0].polyterms[0][0]; + double P1 = coords[0].polyterms[1][0]; + Rc = (Rc + P1)*gamma + P0; + dR = (dR + P1); + + double gamma_new = gamma - (Rc - Ro) / dR; + gamma = gamma_new; + } + + double theta = 90.0 - gamma * DEG_RAD ; + ctht = cos (RAD_DEG * theta); + stht = sin (RAD_DEG * theta); + break; + } + default: return (FALSE); @@ -354,32 +376,30 @@ case PROJ_ZPN: + // the forward projection is: + // theta = atan2(stht, ctht) + // gamma = (pi/2 - theta) : theta in radians + // Ro = sum (P_i gamma^i) + // R = (180/pi) Ro + // Ro = (pi/180)(90 - theta) // R = (180/pi)sum (P_i R^i) + + // is ZPN defined for Npolyterms = 0 or 1? + ctht = hypot(sphi, cphi); theta = atan2 (stht, ctht); - double Rc; - if (UKIRT_ONLY) { - double P1 = coords[0].polyterms[0][1]; - double P3 = coords[0].polyterms[0][3]; - double gamma = RAD_DEG * (90 - DEG_RAD * theta); - Rc = P1*gamma + P3*gamma*gamma*gamma; - } else { - double Ro = RAD_DEG * (90 - DEG_RAD * theta); - - // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21) - i = coords[0].Npolyterms - 1; - Rc = 0.0; - while (i > 0) { - if (i < 7) { - Rc = (Rc + coords[0].polyterms[0][i])*Ro; - } else { - Rc = (Rc + coords[0].polyterms[1][i-7])*Ro; - } - i --; - } - Rc += coords[0].polyterms[0][i]; + double Ro; + double gamma = M_PI_2 - theta; + + // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21) + Ro = 0.0; + for (i = coords[0].Npolyterms - 1; i > 0; i --) { + double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1]; + Ro = (Ro + Pi)*gamma; } - Rc = DEG_RAD * Rc; + Ro += coords[0].polyterms[0][0]; + + Rc = DEG_RAD * Ro; *L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ; @@ -609,7 +629,9 @@ enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN}; +int GetRadialZPN (Coords *coords, Header *header); + int GetCoords (Coords *coords, Header *header) { - int i, status, status1, status2, itmp, Polynomial, Polyterm; + int status, status1, status2, itmp, Polynomial, Polyterm; double Lambda, rotate, rotate1, rotate2, scale; double equinox; @@ -666,8 +688,16 @@ status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2); + ctype = &coords[0].ctype[4]; + + // read the ZPN coeffients PV2_i (i = 0 < 14) + // ZPN is inconsistent with the other Polynomial types + if (!strcmp (ctype, "-ZPN")) { + GetRadialZPN (coords, header); + break; + } + /* set NPLYTERM based on header. if NPLYTERM is missing, it should have a value of 0, unless the projection type is one of PLY, DIS, WRP, in which case it should be set to 3 */ - ctype = &coords[0].ctype[4]; Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP"); Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp); @@ -723,4 +753,7 @@ coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda; coords[0].pc2_2 = cos(rotate*RAD_DEG); + + // read the ZPN coeffients PV2_i (i = 0 < 14) + if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header); break; @@ -744,28 +777,6 @@ coords[0].pc2_2 /= scale; - if (!strcmp (&coords[0].ctype[4], "-ZPN")) { - int found; - for (i = 0; i < 14; i++) { - char name[64]; - snprintf (name, 64, "PV2_%d", i); - if (i < 7) { - found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[0][i]); - } else { - found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[1][i-7]); - } - if ((i == 0) && !found) { - coords[0].polyterms[0][0] = 0.0; - continue; - } - if ((i == 1) && !found) { - coords[0].polyterms[0][1] = 1.0; - continue; - } - if (!found) { - coords[0].Npolyterms = i; - break; - } - } - } + // read the ZPN coeffients PV2_i (i = 0 < 14) + if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header); break; @@ -819,4 +830,37 @@ } return (status); +} + +int GetRadialZPN (Coords *coords, Header *header) { + + // RA---ZPN can have up to 14 radial polynomial terms. these are stored in + // polyterms[0][0] - [6][0] for the first 7 and [0][1] - [6][1] for the rest + // these terms are coeffients of a polynomial of the radial distances + + // read the ZPN coeffients PV2_i (i = 0 < 14) + int found; + int Nmax = 0; + for (int i = 0; i < 14; i++) { + char name[64]; + snprintf (name, 64, "PV2_%d", i); + if (i < 7) { + coords[0].polyterms[i][0] = 0.0; + found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i][0]); + } else { + coords[0].polyterms[i-7][1] = 0.0; + found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i-7][1]); + } + // PV2_1 is implicit if not present + if ((i == 1) && !found) { + coords[0].polyterms[1][0] = 1.0; + continue; + } + // set Npolyterms based on the largest coefficient found + if (found) { + Nmax = i; + } + } + coords[0].Npolyterms = Nmax + 1; + return TRUE; } Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 39670) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 39926) @@ -427,8 +427,8 @@ break; case MEAS_TMEAN: /* OK */ + value.Flt = TimeValue (average[0].Tmean, TimeReference, TimeFormat); + break; + case MEAS_TRANGE: /* OK */ value.Flt = GetTimeRange (average[0].Trange, TimeFormat); - break; - case MEAS_TRANGE: /* OK */ - value.Flt = TimeValue (average[0].Trange, 0, TimeFormat); break; case MEAS_NMEAS: /* OK */ @@ -461,12 +461,15 @@ break; case MEAS_RA_FIT_OFFSET: /* OK */ + // RA_epoch_fit = RA_mean + uR*(t - Tmean)/cos(dec) + plx*parR + // note that this extraction ignores parallax dT = (measure[0].t - average[0].Tmean) / (86400*365.25); - dR = dvoOffsetR (measure, average); - value.Flt = average[0].uR * dT + dR; + dR = dvoOffsetR (measure, average); // RA_epoch - RA_mean (** NOT local linear distance **) + value.Flt = dR*cos(RAD_DEG*measure[0].D) - average[0].uR * dT; + // this is the local linear distance of the measurement from the fit break; case MEAS_DEC_FIT_OFFSET: /* OK */ dT = (measure[0].t - average[0].Tmean) / (86400*365.25); dD = dvoOffsetD (measure, average); - value.Flt = average[0].uD * dT + dD; + value.Flt = dD - average[0].uD * dT; break; case MEAS_RA_OFFSET_ERR: /* OK */ Index: trunk/Ohana/src/libdvo/src/galaxy_model.c =================================================================== --- trunk/Ohana/src/libdvo/src/galaxy_model.c (revision 39670) +++ trunk/Ohana/src/libdvo/src/galaxy_model.c (revision 39926) @@ -28,4 +28,22 @@ V_sol = 11.18; // km/sec W_sol = 7.61; // km/sec + return TRUE; + } + if (!strcmp(version, "TEST-CONSTANT")) { + // use for testing + A_oort = +47.40; // km/sec/kpc + B_oort = -47.40; // km/sec/kpc + U_sol = 0.00; // km/sec + V_sol = 0.00; // km/sec + W_sol = 0.00; // km/sec + return TRUE; + } + if (!strcmp(version, "TEST-ZERO")) { + // use for testing + A_oort = 0.00; // km/sec/kpc + B_oort = 0.00; // km/sec/kpc + U_sol = 0.00; // km/sec + V_sol = 0.00; // km/sec + W_sol = 0.00; // km/sec return TRUE; } Index: trunk/Ohana/src/libdvo/test/coords.sh =================================================================== --- trunk/Ohana/src/libdvo/test/coords.sh (revision 39926) +++ trunk/Ohana/src/libdvo/test/coords.sh (revision 39926) @@ -0,0 +1,62 @@ + +macro fakeimage-tan + + mcreate test 1000 1000 + + keyword test CTYPE1 -w RA---TAN ; keyword test CTYPE1 -wc "Algorithm type for axis 1" + keyword test CTYPE2 -w DEC--TAN ; keyword test CTYPE2 -wc "Algorithm type for axis 2" + keyword test CRPIX1 -wf 2997.03245628555 ; keyword test CRPIX1 -wc "Dither offset Y" + keyword test CRPIX2 -wf -944.140242055343 ; keyword test CRPIX2 -wc "Dither offset Y" + keyword test CRVAL1 -wf 207.605373039722 ; keyword test CRVAL1 -wc "[deg] Right ascension at the reference pixel" + keyword test CRVAL2 -wf 0.113109076969833 ; keyword test CRVAL2 -wc "[deg] Declination at the reference pixel" + keyword test CRUNIT1 -w deg ; keyword test CRUNIT1 -wc "Unit of right ascension co-ordinates" + keyword test CRUNIT2 -w deg ; keyword test CRUNIT2 -wc "Unit of declination co-ordinates" + keyword test CD1_1 -wf -3.34355087531346E-07 ; keyword test CD1_1 -wc "Transformation matrix element" + keyword test CD1_2 -wf -0.000111484161257407 ; keyword test CD1_2 -wc "Transformation matrix element" + keyword test CD2_1 -wf 0.000111383424202905 ; keyword test CD2_1 -wc "Transformation matrix element" + keyword test CD2_2 -wf -3.65044007533155E-07 ; keyword test CD2_2 -wc "Transformation matrix element" + # keyword test PV2_1 -wf 1.000000E+00 ; keyword test -c Pol.coeff. for pixel -> celestial coord + # keyword test PV2_2 -wf 0.000000E+00 ; keyword test -c Pol.coeff. for pixel -> celestial coord + # keyword test PV2_3 -wf -50. + +end + +macro fakeimage-zpn1 + + mcreate test 1000 1000 + + keyword test CTYPE1 -w RA---ZPN ; keyword test CTYPE1 -wc "Algorithm type for axis 1" + keyword test CTYPE2 -w DEC--ZPN ; keyword test CTYPE2 -wc "Algorithm type for axis 2" + keyword test CRPIX1 -wf 2997.03245628555 ; keyword test CRPIX1 -wc "Dither offset Y" + keyword test CRPIX2 -wf -944.140242055343 ; keyword test CRPIX2 -wc "Dither offset Y" + keyword test CRVAL1 -wf 207.605373039722 ; keyword test CRVAL1 -wc "[deg] Right ascension at the reference pixel" + keyword test CRVAL2 -wf 0.113109076969833 ; keyword test CRVAL2 -wc "[deg] Declination at the reference pixel" + keyword test CRUNIT1 -w deg ; keyword test CRUNIT1 -wc "Unit of right ascension co-ordinates" + keyword test CRUNIT2 -w deg ; keyword test CRUNIT2 -wc "Unit of declination co-ordinates" + keyword test CD1_1 -wf -3.34355087531346E-07 ; keyword test CD1_1 -wc "Transformation matrix element" + keyword test CD1_2 -wf -0.000111484161257407 ; keyword test CD1_2 -wc "Transformation matrix element" + keyword test CD2_1 -wf 0.000111383424202905 ; keyword test CD2_1 -wc "Transformation matrix element" + keyword test CD2_2 -wf -3.65044007533155E-07 ; keyword test CD2_2 -wc "Transformation matrix element" + keyword test PV2_1 -wf 1.000000E+00 ; keyword test PV2_1 -wc "Pol.coeff. for pixel -> celestial coord" + keyword test PV2_2 -wf 0.000000E+00 ; keyword test PV2_2 -wc "Pol.coeff. for pixel -> celestial coord" + keyword test PV2_3 -wf -50. + +end + +macro test.coords + for ix 50 5000 500 + for iy 50 5000 500 + coords test -p $ix $iy + coords test -c $RA $DEC + if (abs($ix - $Xc) > 0.01) + echo bad + break + end + if (abs($iy - $Yc) > 0.01) + echo bad + break + end + end + end +end +