Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/Makefile (revision 37601) @@ -49,4 +49,5 @@ $(SRC)/version.$(ARCH).o \ $(SRC)/coordops.$(ARCH).o \ +$(SRC)/galaxy_model.$(ARCH).o \ $(SRC)/dvosorts.$(ARCH).o \ $(SRC)/dvo_photcode_ops.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvo.h (revision 37601) @@ -1038,4 +1038,13 @@ int RegionHostFindNeighbors (RegionHostTable *table, int Nhost); +// galaxy_model: +int TransformProperMotion_radians (double *uR, double *uD, double uL, double uB, double Rrad, double Drad, CoordTransform *transform); +int TransformProperMotion (double *uR, double *uD, double uL, double uB, double R, double D, CoordTransform *transform); +int SolarMotionModel_radians (double *uL_sol, double *uB_sol, double Lrad, double Brad, double distance); +int SolarMotionModel (double *uL_sol, double *uB_sol, double L, double B, double distance); +int GalaxyMotionModel_radians (double *uL_gal, double *uB_gal, double Lrad, double Brad); +int GalaxyMotionModel (double *uL_gal, double *uB_gal, double L, double B); +int InitGalaxyModel (char *version); + # define LENSFIELD(NAME) float LensValue_##NAME (PhotCode *code, Lensobj *lensobj); Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvodb.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvodb.h (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/include/dvodb.h (revision 37601) @@ -173,4 +173,10 @@ MEAS_XFIX, MEAS_YFIX, + MEAS_XOFF_KH, + MEAS_YOFF_KH, + MEAS_XOFF_DCR, + MEAS_YOFF_DCR, + MEAS_ROFF_GAL, + MEAS_DOFF_GAL, MEAS_POS_SYS_ERR, MEAS_XFIELD, Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37601) @@ -557,19 +557,16 @@ # endif break; - case MEAS_XFIX: /* OK */ - value.Flt = measure[0].Xfix; - break; - case MEAS_YFIX: /* OK */ - value.Flt = measure[0].Yfix; - break; - case MEAS_XCCD_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dXccd); - break; - case MEAS_YCCD_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dYccd); - break; - case MEAS_POS_SYS_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dRsys); - break; + case MEAS_XFIX: { value.Flt = measure[0].Xfix; break; } + case MEAS_YFIX: { value.Flt = measure[0].Yfix; break; } + case MEAS_XCCD_ERR: { value.Flt = FromShortPixels(measure[0].dXccd); break; } + case MEAS_YCCD_ERR: { value.Flt = FromShortPixels(measure[0].dYccd); break; } + case MEAS_POS_SYS_ERR: { value.Flt = FromShortPixels(measure[0].dRsys); } + + case MEAS_XOFF_KH: { value.Flt = measure[0].XoffKH; break; } + case MEAS_YOFF_KH: { value.Flt = measure[0].YoffKH; break; } + case MEAS_XOFF_DCR: { value.Flt = measure[0].XoffDCR; break; } + case MEAS_YOFF_DCR: { value.Flt = measure[0].YoffDCR; break; } + case MEAS_ROFF_GAL: { value.Flt = measure[0].RoffGAL; break; } + case MEAS_DOFF_GAL: { value.Flt = measure[0].DoffGAL; break; } case MEAS_XFIELD: /* offset relative to exposure center in ra,dec space */ Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbFields.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbFields.c (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dbFields.c (revision 37601) @@ -474,4 +474,10 @@ if (!strcasecmp (fieldName, "XFIX")) ESCAPE (MEAS_XFIX, OPIHI_FLT); if (!strcasecmp (fieldName, "YFIX")) ESCAPE (MEAS_YFIX, OPIHI_FLT); + if (!strcasecmp (fieldName, "XOFF_KH")) ESCAPE (MEAS_XOFF_KH, OPIHI_FLT); + if (!strcasecmp (fieldName, "YOFF_KH")) ESCAPE (MEAS_YOFF_KH, OPIHI_FLT); + if (!strcasecmp (fieldName, "XOFF_DCR")) ESCAPE (MEAS_XOFF_DCR, OPIHI_FLT); + if (!strcasecmp (fieldName, "YOFF_DCR")) ESCAPE (MEAS_YOFF_DCR, OPIHI_FLT); + if (!strcasecmp (fieldName, "ROFF_GAL")) ESCAPE (MEAS_ROFF_GAL, OPIHI_FLT); + if (!strcasecmp (fieldName, "DOFF_GAL")) ESCAPE (MEAS_DOFF_GAL, OPIHI_FLT); if (!strcasecmp (fieldName, "POS_SYS_ERR")) ESCAPE (MEAS_POS_SYS_ERR, OPIHI_FLT); if (!strcasecmp (fieldName, "XFIELD")) ESCAPE (MEAS_XFIELD, OPIHI_FLT); Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dvo_tiny_values.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 37600) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 37601) @@ -20,4 +20,6 @@ measureT[0].dM = measure[0].dM; measureT[0].airmass = measure[0].airmass; + measureT[0].RoffGAL = measure[0].RoffGAL; + measureT[0].DoffGAL = measure[0].DoffGAL; measureT[0].Xccd = measure[0].Xccd; measureT[0].Yccd = measure[0].Yccd; Index: /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/galaxy_model.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/galaxy_model.c (revision 37601) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/galaxy_model.c (revision 37601) @@ -0,0 +1,148 @@ +# include + +// velocity / distance / proper-motion scale: +static double iFkap = 0.001 / 4.74047; // (km/sec/kpc) / (arcsec/ year) + +static double A_oort = NAN; // km/sec/kpc +static double B_oort = NAN; // km/sec/kpc +static double U_sol = NAN; // km/sec +static double V_sol = NAN; // km/sec +static double W_sol = NAN; // km/sec + +int InitGalaxyModel (char *version) { + + // galaxy model parameters from Siegfreid Roeser: + if (!strcmp(version, "ROESER")) { + A_oort = +14.50; // km/sec/kpc + B_oort = -13.00; // km/sec/kpc + U_sol = 9.44; // km/sec + V_sol = 11.90; // km/sec + W_sol = 7.20; // km/sec + return TRUE; + } + if (!strcmp(version, "FEAST-HIPPARCOS")) { + // http://arxiv.org/pdf/astro-ph/9706293v1.pdf + A_oort = +14.82; // km/sec/kpc + B_oort = -12.37; // km/sec/kpc + U_sol = 9.32; // km/sec + V_sol = 11.18; // km/sec + W_sol = 7.61; // km/sec + return TRUE; + } + return FALSE; +} + +int GalaxyMotionModel (double *uL_gal, double *uB_gal, double L, double B) { + double Lrad = L * RAD_DEG; + double Brad = B * RAD_DEG; + *uL_gal = (A_oort * cos(2.0*Lrad) + B_oort) * cos(Brad) * iFkap; + *uB_gal = -0.5*A_oort * sin(2.0*Lrad) * sin(Brad*2.0) * iFkap; + return TRUE; +} + +int GalaxyMotionModel_radians (double *uL_gal, double *uB_gal, double Lrad, double Brad) { + *uL_gal = (A_oort * cos(2.0*Lrad) + B_oort) * cos(Brad) * iFkap; + *uB_gal = -0.5*A_oort * sin(2.0*Lrad) * sin(Brad*2.0) * iFkap; + return TRUE; +} + +int SolarMotionModel (double *uL_sol, double *uB_sol, double L, double B, double distance) { + double Lrad = L * RAD_DEG; + double Brad = B * RAD_DEG; + *uL_sol = (U_sol * sin(Lrad) - V_sol * cos(Lrad)) * iFkap / distance; + *uB_sol = ((U_sol * cos(Lrad) + V_sol * sin(Lrad))*sin(Brad) - W_sol*cos(Brad)) * iFkap / distance; + return TRUE; +} + +int SolarMotionModel_radians (double *uL_sol, double *uB_sol, double Lrad, double Brad, double distance) { + *uL_sol = (U_sol * sin(Lrad) - V_sol * cos(Lrad)) * iFkap / distance; + *uB_sol = ((U_sol * cos(Lrad) + V_sol * sin(Lrad))*sin(Brad) - W_sol*cos(Brad)) * iFkap / distance; + return TRUE; +} + +// Use the inverse transform: CELESTIAL->GALACTIC for uL,uB -> uR,uD +int TransformProperMotion (double *uR, double *uD, double uL, double uB, double R, double D, CoordTransform *transform) { + + double Rrad = R*RAD_DEG; + double Drad = D*RAD_DEG; + + // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf + double C1 = + cos(Drad)*transform->cos_phi + + sin(Drad)*cos(Rrad)*transform->sin_phi_sin_Xo - + sin(Drad)*sin(Rrad)*transform->sin_phi_cos_Xo; + + double C2 = + - cos(Rrad)*transform->sin_phi_cos_Xo + - sin(Rrad)*transform->sin_phi_sin_Xo; + + double cosBinv = 1.0 / sqrt(C1*C1 + C2*C2); + + // XXX add errors : I need to be able to choose the stars based on the error distribution + *uR = cosBinv * (C1 * uL - C2 * uB); + *uD = cosBinv * (C1 * uB + C2 * uL); + + return TRUE; +} +int TransformProperMotion_radians (double *uR, double *uD, double uL, double uB, double Rrad, double Drad, CoordTransform *transform) { + + // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf + double C1 = + cos(Drad)*transform->cos_phi + + sin(Drad)*cos(Rrad)*transform->sin_phi_sin_Xo - + sin(Drad)*sin(Rrad)*transform->sin_phi_cos_Xo; + + double C2 = + - cos(Rrad)*transform->sin_phi_cos_Xo + - sin(Rrad)*transform->sin_phi_sin_Xo; + + double cosBinv = 1.0 / sqrt(C1*C1 + C2*C2); + + // XXX add errors : I need to be able to choose the stars based on the error distribution + *uR = cosBinv * (C1 * uL - C2 * uB); + *uD = cosBinv * (C1 * uB + C2 * uL); + + return TRUE; +} + +// I am using the galactic rotation and solar motion model to predict the reflex proper +// motion of stars in the database. + +/***** + +From Siegfreid Roeser: + +Dear Gene, + + please find attached our approach for galactic rotation and + solar motion that we assumed in our approach last fall. +(Hope I understood correctly what you requested at the +telecon). If not, please come back. + + fkap = 4.74047d-0 {km/s/kpc}/{mas/y} + + gl,gb is galactic longitude resp. latitude + + + PARAMETER (distance= 1.d0,Aoort=14.5d0 ,Boort=-13.d0 ) + kpc km/s/kpc km/s/kpc + PARAMETER (Usol = 9.44d0, Vsol =11.90d0, Wsol = 7.20d0) + km/s +C +C Gal. Rotation +C + proper motions in longitude/latitude + + emulgal = (AOORT*dcos(2.d0*gl)+BOORT)*dcos(gb)/fkap + emubgal = -0.5d0*AOORT*dsin(2.d0*gl)*dsin(2.d0*gb)/fkap +c +c Solar Motion +cc + emulsol = Usol*dsin(gl) - Vsol*dcos(gl) + emulsol = emulsol/fkap/distance + emubsol = Usol*dsin(gb)*dcos(gl)+Vsol*dsin(gb)*dsin(gl) + emubsol = emubsol -Wsol*dcos(gb) + emubsol = emubsol/fkap/distance + + +******/