Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile (revision 32853) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile (revision 32854) @@ -18,5 +18,5 @@ relastro: $(BIN)/relastro.$(ARCH) -install: $(DESTBIN)/relastro +install: $(DESTBIN)/relastro $(DESTBIN)/testparallax RELASTRO = \ @@ -68,2 +68,14 @@ $(RELASTRO): $(INC)/relastro.h $(KAPA_INCS) $(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS) + +testparallax: $(BIN)/testparallax.$(ARCH) + +TESTPAR = \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/testparallax.$(ARCH).o + +$(TESTPAR): $(INC)/relastro.h $(KAPA_INCS) +$(BIN)/testparallax.$(ARCH): $(TESTPAR) $(KAPA_LIBS) Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/plx_factor.pro =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/plx_factor.pro (revision 32854) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/plx_factor.pro (revision 32854) @@ -0,0 +1,116 @@ +;+ +; +; PLX_FACTOR +; +; Given a Julian date (scalar or vector), this returns the parallax +; factor for a set of equatorial coordinates (RA,Dec). User must have +; an appropriate Earth ephemeris savefile in the path specified at the +; top of the procedure. The savefile should return the variables: +; +; jd, x, y, z +; +; that specify the position of the Earth in the coordinate system of +; solar system's barycenter. For my savefile, I computed this from +; the JPL ephemeris DE405. +; +; INPUT: +; +; epoch (JD) -- julian date(s) for which to compute the parallax +; factor +; +; ra (deg) -- "catalog" RA of object (unperturbed by parallax, +; aberration, etc.) +; +; de (deg) -- "catalog" Dec of object (unperturbed by parallax, +; aberration, etc.) +; +; NOTE: Epoch input can be an array, with RA and Dec input as +; scalars. However, if RA and Dec are input as arrays (i.e., the +; parallax factor is to be computed for multiple objects), then they +; must have the same dimensions as the epoch array. +; +; OUTPUT: +; +; The output produced by this program are two arrays of parallax +; factors (f_ra, f_de) such that: +; +; RA = RA_0 + f_ra * parallax +; Dec = Dec_0 + f_de * parallax +; +; Note that f_ra and f_de have dimensionless units. Also note that +; the RA correction factor includes a factor of cos(Dec) such that: +; f_ra * parallax = delta(RA) * cos(Dec) +; +; Written by Trent J. Dupuy -- 2010 Jul 14 +; +; REVISIONS +; --------- +; 2011 Jan 25: Added /SPITZER keyword that uses the Spitzer ephemeris +; instead of the Earth ephemeris. +;- +pro plx_factor, epoch, ra_in, de_in, f_ra, f_de, spitzer=spitzer + +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;; SPECIFY PATH FOR EPHEMERIS SAVEFILE HERE ;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +savefile = '/Users/tdupuy/data/jpl/'+$ + 'JPL_DE405_EARTH-SOLARBARYCENTER.sav' +if keyword_set(spitzer) then begin + savefile = '/Users/tdupuy/data/jpl/'+$ + 'JPL_SPITZER-SOLARBARYCENTER.sav' + message,'Using spitzer ephemeris.',/continue +endif +restore, savefile +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + +; +; check inputs +; +if n_params() lt 5 then begin + print, 'PLX_FACTOR, epoch (JD), ra (deg), dec (deg), f_ra, f_de' + return +endif +if total(size(ra_in,/dim) ne size(de_in,/dim)) ne 0 then $ + message, 'input RA and Dec dimensions must match each other' +if (size(ra_in,/dim))[0] ne 0 and $ + total(size(ra_in,/dim) ne size(epoch,/dim)) ne 0 then $ + message, '(RA, Dec) not specified as scalars, so their '+$ + 'dimensions must match the input epoch array' + +; +; convert to radians and compute trig functions +; +radeg = 180d/!dpi +ra = ra_in/radeg +de = de_in/radeg +sina = sin(ra) +cosa = cos(ra) +sind = sin(de) +cosd = cos(de) + +; +; interpolate the Earth ephemeris at the input epochs +; +x_in = interpol( x, jd, double(epoch) ) +y_in = interpol( y, jd, double(epoch) ) +z_in = interpol( z, jd, double(epoch) ) +if min(epoch) lt min(jd) or max(epoch) gt max(jd) then $ + message, /continue, 'inputs epochs extend beyond the range '+$ + 'of the ephemeris being used.' + +; +; could also apply an offset here due to the position of the observer +; on the earth being slightly different than the geocenter. as +; described in Sec 6.2.3 of Kovalesky & Seidelmann this effect is +; about 4 uas for a parallax of 0.1" -- i.e., negligible +; + +; +; compute parallax factor (from pg. B28 of Astronomical Almanac 2010, +; also see Sec 6.2.1 of Kovalesky & Seidelmann's Fundamentals +; of Astrometry for a derivation) +; +f_ra = x_in*sina - y_in*cosa +f_de = x_in*cosa*sind + y_in*sina*sind - z_in*cosd + +end Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ParFactor.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ParFactor.c (revision 32853) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ParFactor.c (revision 32854) @@ -1,3 +1,4 @@ # include "relastro.h" +# define J2000 2451545. /* Julian date at standard epoch */ # if (0) @@ -31,5 +32,4 @@ double n, L, g; -# define J2000 2451545. /* Julian date at standard epoch */ n = jd - J2000; @@ -46,6 +46,4 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) { -# define J2000 2451545. /* Julian date at standard epoch */ - double n = jd - J2000; // day number double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) @@ -60,4 +58,5 @@ /* given RA, DEC, Time, calculate the parallax factor */ +// Time is relative to Tmean, Tmean is years relative to J2000 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) { @@ -66,5 +65,8 @@ /* given a time T in UNIX seconds, determine the solar longitude S */ - jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); + // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); + jd = 365.25*(Time + Tmean) + J2000; + // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd); + sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 32853) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 32854) @@ -133,4 +133,7 @@ R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + + // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects + // to get Time in years relative to UNIX Tzero (1970/01/01) T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/testparallax.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/testparallax.c (revision 32853) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/testparallax.c (revision 32854) @@ -1,31 +1,114 @@ # include "relastro.h" +# define NPTS 64 +# define J2000 2451545. /* Julian date at standard epoch */ +# define T1970 2440587.500000 /* JD at UNIX ref time */ int main (int argc, char **argv) { + + int i, Npts; + char line[1024]; + double R[NPTS], D[NPTS], Time[NPTS], JD[NPTS]; + double X[NPTS], Y[NPTS], dX[NPTS], dY[NPTS], pX[NPTS], pY[NPTS]; + double Tref[NPTS], Tjyrs[NPTS], TrefS, TjyrsS, TrefMean, TjyrsMean, Ro, Do; + Coords coords; + PMFit fitPAR; - if (argc != N) { - fprintf (stderr, "USAGE: testparallax (file.dat)\n"); - exit (1); + if (argc != 2) { + fprintf (stderr, "USAGE: testparallax (file.dat)\n"); + exit (1); + } + + // test parallax program + FILE *f = fopen (argv[1], "r"); + if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); } + + // XXX uncomment to skip first line scan_line(f, line); + + Npts = 0; + for (i = 0; scan_line(f, line) != EOF; i++) { + if (Npts == NPTS) { + fprintf (stderr, "too many point: use dynamic alloc\n"); + exit (2); } + if (line[0] == '#') continue; + dparse (&R[Npts], 1, line); + dparse (&D[Npts], 2, line); + dparse (&JD[Npts], 3, line); + dparse (&Time[Npts], 3, line); - // test parallax program - FILE *f = fopen (argv[1], "r"); - if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); } + dparse (&dX[Npts], 4, line); + dparse (&dY[Npts], 4, line); - char buffer[1024]; - scan_line(f, line); - - while (scan_line(f, line) != EOF) { - dparse (&jd, 1, line); - dparse (&time, 2, line); - dparse (&ra, 3, line); - dparse (&dec, 4, line); + // trent's file + // dparse (&JD[Npts], 1, line); + // dparse (&Time[Npts], 2, line); + // dparse (&R[Npts], 3, line); + // dparse (&D[Npts], 4, line); + Npts ++; + } - } + /* project coordinates to a plane centered on the object with units of arcsec */ + coords.crval1 = 0; + coords.crval2 = 0; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "RA---SIN"); - // project to local coords + // use one point as a local reference + coords.crval1 = R[0]; + coords.crval2 = D[0]; - // run fitter - FitPMandPar + TrefS = TjyrsS = 0.0; + // project to local coords + for (i = 0; i < Npts; i++) { + RD_to_XY (&X[i], &Y[i], R[i], D[i], &coords); + Tjyrs[i] = Time[i] / 365.25; + Tref[i] = Time[i] / 365.25; + // Tunix[i] = (JD[i] - T1970) / 365.25; // time relative to T1970 in years + // Tjyrs[i] = (JD[i] - J2000) / 365.25; // time relative to J2000 in years + + TrefS += Tref[i]; + TjyrsS += Tjyrs[i]; + + // nominal for PS1 + // dX[i] = 0.020; + // dY[i] = 0.020; + } + TrefMean = TrefS / Npts; + TjyrsMean = TjyrsS / Npts; + + for (i = 0; i < Npts; i++) { + Tref[i] -= TrefMean; + ParFactor (&pX[i], &pY[i], R[i], D[i], Tjyrs[i], 0.0); + fprintf (stderr, "%f %f : %f %f : %f %f : %f %f\n", R[i], D[i], X[i], Y[i], Tref[i], Tjyrs[i], pX[i], pY[i]); + } + + // run fitter + // FitPMandPar (&fitPAR, X, dX, Y, dY, Tref, pX, pY, Npts, TRUE); + FitPMandPar (&fitPAR, X, dX, Y, dY, Tjyrs, pX, pY, Npts, TRUE); + XY_to_RD (&Ro, &Do, fitPAR.Ro, fitPAR.Do, &coords); + + fprintf (stderr, "Rx : %f\n", fitPAR.Ro); + fprintf (stderr, "Dx : %f\n", fitPAR.Do); + fprintf (stderr, "Ro : %f\n", Ro); + fprintf (stderr, "Do : %f\n", Do); + fprintf (stderr, "dRo : %f\n", fitPAR.dRo); + fprintf (stderr, "dDo : %f\n", fitPAR.dDo); + fprintf (stderr, "uR : %f\n", fitPAR.uR); + fprintf (stderr, "uD : %f\n", fitPAR.uD); + fprintf (stderr, "duR : %f\n", fitPAR.duR); + fprintf (stderr, "duD : %f\n", fitPAR.duD); + fprintf (stderr, "p : %f\n", fitPAR.p); + fprintf (stderr, "dp : %f\n", fitPAR.dp); + fprintf (stderr, "uTot : %f\n", hypot(fitPAR.uR,fitPAR.uD)); + fprintf (stderr, "PA : %f\n", DEG_RAD*atan2(fitPAR.uR,fitPAR.uD) + 360.0); + + exit (0); } +