Index: /trunk/Ohana/src/tools/src/mpcorb_predict.c =================================================================== --- /trunk/Ohana/src/tools/src/mpcorb_predict.c (revision 41710) +++ /trunk/Ohana/src/tools/src/mpcorb_predict.c (revision 41711) @@ -33,4 +33,13 @@ } Planets; +typedef struct { + char ID[8]; + double Robs; + double Dobs; + double Rvel; + double Dvel; + double dist; +} PlanetDatum; + int Shutdown (char *format, ...); int mpcorb_parseline (char *line, Planets *planet, int Nmax); @@ -38,7 +47,7 @@ Planets *mpcorb_read_fits (char *filename, int *nplanets); void mpcorb_save_fits (char *filename, Planets *planets, int Nplanets); -void mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets); - -void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc); +void mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets); + +PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc); void mpcorb_trange (int argc, char **argv); @@ -54,7 +63,41 @@ # define PS1_ALTITUDE 3067.7 /* meters */ +// some options (for testing) +int USE_AMP = FALSE; // otherwise use precess and/or nutation +int USE_PLANETARY_ABERRATION = FALSE; +int USE_MANUAL_PRECESS = FALSE; // otherwise use slaPreces +int USE_PRENUT = FALSE; // otherwise slaPrec + int main (int argc, char **argv) { + int N; + if (TESTING) mpcorb_testing (argc, argv); + + // -none option so I can always supply an option in dvo script + if ((N = get_argument (argc, argv, "-none"))) { + remove_argument (N, &argc, argv); + } + + // turn on / off some options + if ((N = get_argument (argc, argv, "-use-amp"))) { + USE_AMP = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-apply-planet-ab"))) { + USE_PLANETARY_ABERRATION = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-use-manual-precess"))) { + if (USE_AMP) Shutdown ("-use-amp and -use-manual-precess are incompatible"); + USE_MANUAL_PRECESS = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-use-prenut"))) { + if (USE_AMP) Shutdown ("-use-amp and -use-prenut are incompatible"); + if (!USE_MANUAL_PRECESS) Shutdown ("-use-prenut requires -use-manual-precess"); + USE_PRENUT = TRUE; + remove_argument (N, &argc, argv); + } if ((argc != 9) && (argc != 10)) goto usage; @@ -100,20 +143,18 @@ // is the object in the region for the start of the period? - double RatMin, DatMin; - mpcorb_predict (&planets[i], mjdMin, &RatMin, &DatMin, FALSE); + PlanetDatum plMin = mpcorb_predict (&planets[i], mjdMin, FALSE); int inRangeMin = TRUE; - inRangeMin = inRangeMin && (RatMin > Rmin); - inRangeMin = inRangeMin && (RatMin < Rmax); - inRangeMin = inRangeMin && (DatMin > Dmin); - inRangeMin = inRangeMin && (DatMin < Dmax); + inRangeMin = inRangeMin && (plMin.Robs > Rmin); + inRangeMin = inRangeMin && (plMin.Robs < Rmax); + inRangeMin = inRangeMin && (plMin.Dobs > Dmin); + inRangeMin = inRangeMin && (plMin.Dobs < Dmax); // is the object in the region for the end of the period? - double RatMax, DatMax; - mpcorb_predict (&planets[i], mjdMax, &RatMax, &DatMax, FALSE); + PlanetDatum plMax = mpcorb_predict (&planets[i], mjdMax, FALSE); int inRangeMax = TRUE; - inRangeMax = inRangeMax && (RatMax > Rmin); - inRangeMax = inRangeMax && (RatMax < Rmax); - inRangeMax = inRangeMax && (DatMax > Dmin); - inRangeMax = inRangeMax && (DatMax < Dmax); + inRangeMax = inRangeMax && (plMax.Robs > Rmin); + inRangeMax = inRangeMax && (plMax.Robs < Rmax); + inRangeMax = inRangeMax && (plMax.Dobs > Dmin); + inRangeMax = inRangeMax && (plMax.Dobs < Dmax); // NOTE: if an object is moving so fast that it traverses the entire region, we will @@ -123,9 +164,9 @@ planetsSave[NplanetsSave] = planets[i]; planetsSave[NplanetsSave].mjdMin = mjdMin; - planetsSave[NplanetsSave].RatMin = RatMin; - planetsSave[NplanetsSave].DatMin = DatMin; + planetsSave[NplanetsSave].RatMin = plMin.Robs; + planetsSave[NplanetsSave].DatMin = plMin.Dobs; planetsSave[NplanetsSave].mjdMax = mjdMax; - planetsSave[NplanetsSave].RatMax = RatMax; - planetsSave[NplanetsSave].DatMax = DatMax; + planetsSave[NplanetsSave].RatMax = plMax.Robs; + planetsSave[NplanetsSave].DatMax = plMax.Dobs; NplanetsSave++; @@ -134,8 +175,11 @@ fprintf (stderr, "\n"); + // NOTE : here we are saving the full orbital elements mpcorb_save_fits (output, planetsSave, NplanetsSave); exit (0); } + +# define VERBOSE_TEST FALSE // generate a table of asteroid positions which fall within the region at Tobs (in MJD) @@ -152,7 +196,4 @@ char *output = argv[8]; - // FILE *fout = fopen (output, "w"); - // if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output); - int Nplanets = 0; Planets *planets = mpcorb_read_fits (filename, &Nplanets); @@ -161,5 +202,13 @@ int NplanetsSave = 0; int NPLANETSSAVE = 1000; - ALLOCATE_PTR (planetsSave, Planets, NPLANETSSAVE); + ALLOCATE_PTR (planetsSave, PlanetDatum, NPLANETSSAVE); + + if (VERBOSE_TEST) { + // XXX for a test, I print out coords at steps in myAmpqk + for (int i = 0; i < 10; i++) { + mpcorb_predict (&planets[i], mjdObs, TRUE); + } + exit (0); + } // transform all objects and identify those in the target region: @@ -169,27 +218,21 @@ // is the object in the region for the start of the period? // first, use the fast, but inaccurate, calculation - double Robs, Dobs; - mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, FALSE); - if (Robs < Rmin) continue; - if (Robs > Rmax) continue; - if (Dobs < Dmin) continue; - if (Dobs > Dmax) continue; + PlanetDatum planetFast = mpcorb_predict (&planets[i], mjdObs, FALSE); + if (planetFast.Robs < Rmin) continue; + if (planetFast.Robs > Rmax) continue; + if (planetFast.Dobs < Dmin) continue; + if (planetFast.Dobs > Dmax) continue; // if it is in the region, use the more accurate calculation - mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE); + PlanetDatum planetSlow = mpcorb_predict (&planets[i], mjdObs, TRUE); // fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs); // save this prediction, with R,D in both @min & @max - planetsSave[NplanetsSave] = planets[i]; - planetsSave[NplanetsSave].mjdMin = mjdObs; - planetsSave[NplanetsSave].RatMin = Robs; - planetsSave[NplanetsSave].DatMin = Dobs; - planetsSave[NplanetsSave].mjdMax = mjdObs; - planetsSave[NplanetsSave].RatMax = Robs; - planetsSave[NplanetsSave].DatMax = Dobs; + planetsSave[NplanetsSave] = planetSlow; + strcpy (planetsSave[NplanetsSave].ID, planets[i].ID); // ID is not passed to mpcorb_predict NplanetsSave++; - CHECK_REALLOCATE (planetsSave, Planets, NPLANETSSAVE, NplanetsSave, 1000); + CHECK_REALLOCATE (planetsSave, PlanetDatum, NPLANETSSAVE, NplanetsSave, 1000); } fprintf (stderr, "\n"); @@ -204,5 +247,5 @@ // write minor planet positions (Rref,Dref,Mref) to a FITS table -void mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets) { +void mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets) { Header header; @@ -218,8 +261,11 @@ gfits_create_table_header (&theader, "BINTABLE", "DATA"); - gfits_define_bintable_column (&theader, "8A", "ID", "name", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "Rref", "RA", "degrees", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "Dref", "DEC", "degrees", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "Mref", "Mag", "magnitudes", 1.0, 0.0); + gfits_define_bintable_column (&theader, "8A", "ID", "name", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Rref", "RA", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Dref", "DEC", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Rvel", "RA_VEL", "arcsec/min", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Dvel", "DEC_VEL", "arcsec/min", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "dist", "DIST", "AU", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "Mref", "Mag", "magnitudes", 1.0, 0.0); // generate the output array that carries the data @@ -230,10 +276,16 @@ ALLOCATE_PTR (Rref, double, Nplanets); ALLOCATE_PTR (Dref, double, Nplanets); + ALLOCATE_PTR (Rvel, double, Nplanets); + ALLOCATE_PTR (Dvel, double, Nplanets); + ALLOCATE_PTR (dist, double, Nplanets); ALLOCATE_PTR (Mref, double, Nplanets); // set intermediate storage arrays for (int i = 0; i < Nplanets; i++) { - Rref[i] = planets[i].RatMin; - Dref[i] = planets[i].DatMin; + Rref[i] = planets[i].Robs; + Dref[i] = planets[i].Dobs; + Rvel[i] = planets[i].Rvel; + Dvel[i] = planets[i].Dvel; + dist[i] = planets[i].dist; Mref[i] = 16.0; // need to add this to prediction memcpy(&ID[i*8], planets[i].ID, 8); @@ -244,4 +296,7 @@ gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets); gfits_set_bintable_column (&theader, &ftable, "Dref", Dref, Nplanets); + gfits_set_bintable_column (&theader, &ftable, "Rvel", Rvel, Nplanets); + gfits_set_bintable_column (&theader, &ftable, "Dvel", Dvel, Nplanets); + gfits_set_bintable_column (&theader, &ftable, "dist", dist, Nplanets); gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets); @@ -250,4 +305,7 @@ free (Rref); free (Dref); + free (Rvel); + free (Dvel); + free (dist); free (Mref); @@ -530,11 +588,10 @@ if (TESTING == 2) { - double Robs, Dobs; for (int i = 0; i < 10; i++) { - mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE); + PlanetDatum myPlanet = mpcorb_predict (&planets[i], mjdObs, TRUE); char ra_string[64], de_string[64]; - hms_format (ra_string, 64, Robs/15.0); - hms_format (de_string, 64, Dobs); - fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", Robs, Dobs, ra_string, de_string); + hms_format (ra_string, 64, myPlanet.Robs/15.0); + hms_format (de_string, 64, myPlanet.Dobs); + fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", myPlanet.Robs, myPlanet.Dobs, ra_string, de_string); } exit (0); @@ -674,6 +731,8 @@ } +void myAmp ( double ra, double da, double date, double eq, double *rm, double *dm ); + // mjdObs is NOT in TT -void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc) { +PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc) { int jstat; @@ -683,6 +742,13 @@ Planets tmpPlanet; + PlanetDatum myPlanet; + myPlanet.Robs = NAN; + myPlanet.Dobs = NAN; + myPlanet.Rvel = NAN; + myPlanet.Dvel = NAN; + myPlanet.dist = NAN; // correct the orbital elements to the perturbed version + tmpPlanet = *planet; if (FullCalc) { slaPertel (2, planet->epoch, @@ -704,5 +770,12 @@ &jstat); } else { - tmpPlanet = *planet; + // the function calls below require elements in radians + tmpPlanet.epoch = planet->epoch; + tmpPlanet.inclination = planet->inclination *RAD_DEG; + tmpPlanet.ascend_node = planet->ascend_node *RAD_DEG; + tmpPlanet.perihelion = planet->perihelion *RAD_DEG; + tmpPlanet.semimajor_a = planet->semimajor_a ; + tmpPlanet.eccentricity = planet->eccentricity; + tmpPlanet.mean_anomaly = planet->mean_anomaly*RAD_DEG; } @@ -715,12 +788,47 @@ PS1_LATITUDE*RAD_DEG, 2, - planet->epoch, - planet->inclination*RAD_DEG, - planet->ascend_node*RAD_DEG, - planet->perihelion*RAD_DEG, - planet->semimajor_a, - planet->eccentricity, - planet->mean_anomaly*RAD_DEG, + tmpPlanet.epoch, + tmpPlanet.inclination, + tmpPlanet.ascend_node, + tmpPlanet.perihelion, + tmpPlanet.semimajor_a, + tmpPlanet.eccentricity, + tmpPlanet.mean_anomaly, 0.0, &Rtopo, &Dtopo, &dist, &jstat); + myPlanet.dist = dist; + + double dRobs = 0.0; // planetary aberration, if FullCalc + double dDobs = 0.0; // planetary aberration, if FullCalc + if (FullCalc) { + // calculate the velocity components + + // find position one minute later + double Rtop2, Dtop2, dist2; + slaPlante (mjdObsTT + 60/86400.0, + PS1_LONGITUDE*RAD_DEG, + PS1_LATITUDE*RAD_DEG, + 2, + tmpPlanet.epoch, + tmpPlanet.inclination, + tmpPlanet.ascend_node, + tmpPlanet.perihelion, + tmpPlanet.semimajor_a, + tmpPlanet.eccentricity, + tmpPlanet.mean_anomaly, + 0.0, &Rtop2, &Dtop2, &dist2, &jstat); + + myPlanet.Rvel = 3600.0*(Rtop2 - Rtopo)*DEG_RAD*cos(Dtopo); // this vector is wrong by the precession rotation + myPlanet.Dvel = 3600.0*(Dtop2 - Dtopo)*DEG_RAD; // this vector is wrong by the precession rotation + + // include the correction for planetary aberration + // XXX make these constants more accurate + double dist_km = dist * 1.5e8; + double time_min = dist_km / 299792.459/ 60.0; + double Rlag = myPlanet.Rvel * time_min; // velocity is in arcsec / min + double Dlag = myPlanet.Dvel * time_min; // velocity is in arcsec / min + + dRobs = Rlag / 3600.0 / cos(Dtopo); // RA degrees + dDobs = Dlag / 3600.0; // DEC degrees + } // NOTE: the atm effects seem to be large and incompatible with MPC outputs. @@ -749,9 +857,186 @@ double Rmean, Dmean; if (FullCalc) { - slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean); + if (USE_AMP) { + // slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean); + myAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean); + } else { + // use slaPreces and do not correct for stellar aberration + // convert mjdObsTT to year + double epochYear = 2000.0 + (mjdObsTT - 51544.5) / 365.25; // J2000 = 51544.5 + Rmean = Rap; + Dmean = Dap; + + // use slalib precession and nutation: + if (USE_MANUAL_PRECESS) { + double pm[3][3], v1[3], v2[3]; + + /* Convert RA,Dec to x,y,z */ + slaDcs2c ( Rmean, Dmean, v1 ); + + // calculate the precession & nutation combined matrix: + if (USE_PRENUT) { + // slaPrenut ( epochYear, 51544.5, pm ); + slaPrenut ( 2000.0, mjdObsTT, pm ); // prenut is always mean-to-apparent + slaDimxv ( pm, v1, v2 ); // apply the inverse transformation + } else { + slaPrec ( epochYear, 2000.0, pm ); // precess from apparent-to-mean + slaDmxv ( pm, v1, v2 ); // apply the forward transformatio + } + + /* Back to RA,Dec */ + slaDcc2s ( v2, &Rmean, &Dmean ); + + // renorm to 0-360 + Rmean = slaDranrm ( Rmean ); + } else { + slaPreces ("FK5", epochYear, 2000.000000, &Rmean, &Dmean); + } + // fprintf (stderr, "year: %f\n", epochYear); + } } else {Rmean = Rap; Dmean = Dap; } - *Robs = Rmean * DEG_RAD; - *Dobs = Dmean * DEG_RAD; + if (USE_PLANETARY_ABERRATION) { + myPlanet.Robs = Rmean * DEG_RAD + dRobs; + myPlanet.Dobs = Dmean * DEG_RAD + dDobs; + } else { + myPlanet.Robs = Rmean * DEG_RAD; + myPlanet.Dobs = Dmean * DEG_RAD; + } + return myPlanet; +} + +# define DISABLE_ABERRATION FALSE +# define DISABLE_DEFLECTION FALSE + +void myAmpqk ( double ra, double da, double amprms[21], double *rm, double *dm ) { + double gr2e; /* (grav rad Sun)*2/(Sun-Earth distance) */ + double ab1; /* sqrt(1-v*v) where v=modulus of Earth vel */ + double ehn[3]; /* Earth position wrt Sun (unit vector, FK5) */ + double abv[3]; /* Earth velocity wrt SSB (c, FK5) */ + double p[3], p1[3], p2[3], p3[3]; /* work vectors */ + double ab1p1, p1dv, p1dvp1, w, pde, pdep1; + int i, j; + +/* Unpack some of the parameters */ + gr2e = amprms[7]; + ab1 = amprms[11]; + for ( i = 0; i < 3; i++ ) { + ehn[i] = amprms[i + 4]; + abv[i] = amprms[i + 8]; + } + + // invert and print out: + if (VERBOSE_TEST) { + double Rmean, Dmean; + Rmean = ra * DEG_RAD; + Dmean = da * DEG_RAD; + fprintf (stderr, "%12.6f %12.6f ", Rmean, Dmean); + } + +/* Apparent RA,Dec to Cartesian */ + slaDcs2c ( ra, da, p3 ); + + // invert and print out: + if (VERBOSE_TEST) { + double Rtmp, Dtmp, Rmean, Dmean; + slaDcc2s ( p3, &Rtmp, &Dtmp ); + Rtmp = slaDranrm ( Rtmp ); + Rmean = Rtmp * DEG_RAD; + Dmean = Dtmp * DEG_RAD; + fprintf (stderr, "C: %12.6f %12.6f ", Rmean, Dmean); + } + +/* Precession and nutation */ + slaDimxv ( (double(*)[3]) &rms[12], p3, p2 ); + + // invert and print out: + if (VERBOSE_TEST) { + double Rtmp, Dtmp, Rmean, Dmean; + slaDcc2s ( p2, &Rtmp, &Dtmp ); + Rtmp = slaDranrm ( Rtmp ); + Rmean = Rtmp * DEG_RAD; + Dmean = Dtmp * DEG_RAD; + fprintf (stderr, "P: %12.6f %12.6f ", Rmean, Dmean); + } + + if (DISABLE_ABERRATION) { + // XXX : save result from p2 in p1 + for ( i = 0; i < 3; i++ ) { + p1[i] = p2[i]; + } + } else { + /* Aberration */ + ab1p1 = ab1 + 1.0; + for ( i = 0; i < 3; i++ ) { + p1[i] = p2[i]; + } + for ( j = 0; j < 2; j++ ) { + p1dv = slaDvdv ( p1, abv ); + p1dvp1 = 1.0 + p1dv; + w = 1.0 + p1dv / ab1p1; + for ( i = 0; i < 3; i++ ) { + p1[i] = ( p1dvp1 * p2[i] - w * abv[i] ) / ab1; + } + slaDvn ( p1, p3, &w ); + for ( i = 0; i < 3; i++ ) { + p1[i] = p3[i]; + } + } + } + + // invert and print out: + if (VERBOSE_TEST) { + double Rtmp, Dtmp, Rmean, Dmean; + slaDcc2s ( p1, &Rtmp, &Dtmp ); + Rtmp = slaDranrm ( Rtmp ); + Rmean = Rtmp * DEG_RAD; + Dmean = Dtmp * DEG_RAD; + fprintf (stderr, "A: %12.6f %12.6f ", Rmean, Dmean); + } + +/* Light deflection */ + for ( i = 0; i < 3; i++ ) { + p[i] = p1[i]; + } + +// disable deflection: + if (DISABLE_DEFLECTION) { + + } else { + for ( j = 0; j < 5; j++ ) { + pde = slaDvdv ( p, ehn ); + pdep1 = 1.0 + pde; + w = pdep1 - gr2e * pde; + for ( i = 0; i < 3; i++ ) { + p[i] = ( pdep1 * p1[i] - gr2e * ehn[i] ) / w; + } + slaDvn ( p, p2, &w ); + for ( i = 0; i < 3; i++ ) { + p[i] = p2[i]; + } + } + } + + // invert and print out: + if (VERBOSE_TEST) { + double Rtmp, Dtmp, Rmean, Dmean; + slaDcc2s ( p, &Rtmp, &Dtmp ); + Rtmp = slaDranrm ( Rtmp ); + Rmean = Rtmp * DEG_RAD; + Dmean = Dtmp * DEG_RAD; + fprintf (stderr, "D: %12.6f %12.6f\n", Rmean, Dmean); + } + +/* Mean RA,Dec */ + slaDcc2s ( p, rm, dm ); + *rm = slaDranrm ( *rm ); +} + +void myAmp ( double ra, double da, double date, double eq, double *rm, double *dm ) { + + double amprms[21]; /* Mean-to-apparent parameters */ + + slaMappa ( eq, date, amprms ); + myAmpqk ( ra, da, amprms, rm, dm ); } @@ -762,5 +1047,5 @@ // fields I need: - + // epoch : 21 - 25 // inclination : 60 - 68 @@ -912,9 +1197,8 @@ testPlanet.semimajor_a = 2.782049599999998; - double Robs, Dobs; - mpcorb_predict (&testPlanet, epochUT, &Robs, &Dobs, TRUE); - - fprintf (stderr, "%f vs %f : %f\n", 3.92105439135726600*DEG_RAD, Robs, 3600*(3.92105439135726600*DEG_RAD - Robs)); - fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, Dobs, 3600*(-0.33499850519808244*DEG_RAD - Dobs)); + PlanetDatum myPlanet = mpcorb_predict (&testPlanet, epochUT, TRUE); + + fprintf (stderr, "%f vs %f : %f\n", 3.92105439135726600*DEG_RAD, myPlanet.Robs, 3600*(3.92105439135726600*DEG_RAD - myPlanet.Robs)); + fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, myPlanet.Dobs, 3600*(-0.33499850519808244*DEG_RAD - myPlanet.Dobs)); exit (0); }