IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41713


Ignore:
Timestamp:
Jul 12, 2021, 1:40:08 PM (5 years ago)
Author:
eugene
Message:

cleanup mpcorb_predict code; move test code to src/mpcorb_predict_tests.c

Location:
trunk/Ohana/src/tools/src
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/tools/src/mpcorb_predict.c

    r41711 r41713  
    1616# define TESTING 0
    1717
     18// structure to represent the model of a planet
    1819typedef struct {
    1920  char   ID[8];
    … …  
    3334} Planets;
    3435
     36// structure to represent properties of an observed planet
    3537typedef struct {
    3638  char   ID[8];
    3739  double Robs;
    3840  double Dobs;
    39   double Rvel;
    40   double Dvel;
    4141  double dist;
    4242} PlanetDatum;
    … …  
    5454void mpcorb_moment (int argc, char **argv);
    5555
    56 void mpcorb_testing (int argc, char **argv);
    57 
    58 void mpcorb_predict_test ();
    59 void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs);
    60 
    6156# define PS1_LONGITUDE -156.2559041668965 /* degrees East */
    6257# define PS1_LATITUDE    20.7070994446013 /* degrees North */
    6358# define PS1_ALTITUDE  3067.7 /* meters */
    6459
    65 // some options (for testing)
    66 int USE_AMP                  = FALSE; // otherwise use precess and/or nutation
    67 int USE_PLANETARY_ABERRATION = FALSE;
    68 int USE_MANUAL_PRECESS       = FALSE; // otherwise use slaPreces
    69 int USE_PRENUT               = FALSE; // otherwise slaPrec
    70 
    7160int main (int argc, char **argv) {
    72 
    73   int N;
    74 
    75   if (TESTING) mpcorb_testing (argc, argv);
    76 
    77   // -none option so I can always supply an option in dvo script
    78   if ((N = get_argument (argc, argv, "-none"))) {
    79     remove_argument (N, &argc, argv);
    80   }
    81 
    82   // turn on / off some options
    83   if ((N = get_argument (argc, argv, "-use-amp"))) {
    84     USE_AMP = TRUE;
    85     remove_argument (N, &argc, argv);
    86   }
    87   if ((N = get_argument (argc, argv, "-apply-planet-ab"))) {
    88     USE_PLANETARY_ABERRATION = TRUE;
    89     remove_argument (N, &argc, argv);
    90   }
    91   if ((N = get_argument (argc, argv, "-use-manual-precess"))) {
    92     if (USE_AMP) Shutdown ("-use-amp and -use-manual-precess are incompatible");
    93     USE_MANUAL_PRECESS = TRUE;
    94     remove_argument (N, &argc, argv);
    95   }
    96   if ((N = get_argument (argc, argv, "-use-prenut"))) {
    97     if (USE_AMP) Shutdown ("-use-amp and -use-prenut are incompatible");
    98     if (!USE_MANUAL_PRECESS) Shutdown ("-use-prenut requires -use-manual-precess");
    99     USE_PRENUT = TRUE;
    100     remove_argument (N, &argc, argv);
    101   }
    10261
    10362  if ((argc != 9) && (argc != 10)) goto usage;
    … …  
    181140}
    182141
    183 # define VERBOSE_TEST FALSE
    184 
    185142// generate a table of asteroid positions which fall within the region at Tobs (in MJD)
    186143void mpcorb_moment (int argc, char **argv) {
    … …  
    203160  int NPLANETSSAVE = 1000;
    204161  ALLOCATE_PTR (planetsSave, PlanetDatum, NPLANETSSAVE);
    205 
    206   if (VERBOSE_TEST) {
    207     // XXX for a test, I print out coords at steps in myAmpqk
    208     for (int i = 0; i < 10; i++) {
    209       mpcorb_predict (&planets[i], mjdObs, TRUE);
    210     }
    211     exit (0);
    212   }
    213162
    214163  // transform all objects and identify those in the target region:
    … …  
    227176    PlanetDatum planetSlow = mpcorb_predict (&planets[i], mjdObs, TRUE);
    228177
    229     // fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
    230 
    231178    // save this prediction, with R,D in both @min & @max
    232179    planetsSave[NplanetsSave] = planetSlow;
    … …  
    264211  gfits_define_bintable_column (&theader,  "D", "Rref", "RA",      "degrees",     1.0, 0.0);
    265212  gfits_define_bintable_column (&theader,  "D", "Dref", "DEC",     "degrees",     1.0, 0.0);
    266   gfits_define_bintable_column (&theader,  "D", "Rvel", "RA_VEL",  "arcsec/min",  1.0, 0.0);
    267   gfits_define_bintable_column (&theader,  "D", "Dvel", "DEC_VEL", "arcsec/min",  1.0, 0.0);
    268213  gfits_define_bintable_column (&theader,  "D", "dist", "DIST",    "AU",          1.0, 0.0);
    269214  gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",     "magnitudes",  1.0, 0.0);
    … …  
    276221  ALLOCATE_PTR (Rref,                 double, Nplanets);
    277222  ALLOCATE_PTR (Dref,                 double, Nplanets);
    278   ALLOCATE_PTR (Rvel,                 double, Nplanets);
    279   ALLOCATE_PTR (Dvel,                 double, Nplanets);
    280223  ALLOCATE_PTR (dist,                 double, Nplanets);
    281224  ALLOCATE_PTR (Mref,                 double, Nplanets);
    … …  
    285228    Rref[i]        = planets[i].Robs;
    286229    Dref[i]        = planets[i].Dobs;
    287     Rvel[i]        = planets[i].Rvel;
    288     Dvel[i]        = planets[i].Dvel;
    289230    dist[i]        = planets[i].dist;
    290231    Mref[i]        = 16.0; // need to add this to prediction
    … …  
    296237  gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets);
    297238  gfits_set_bintable_column (&theader, &ftable, "Dref", Dref, Nplanets);
    298   gfits_set_bintable_column (&theader, &ftable, "Rvel", Rvel, Nplanets);
    299   gfits_set_bintable_column (&theader, &ftable, "Dvel", Dvel, Nplanets);
    300239  gfits_set_bintable_column (&theader, &ftable, "dist", dist, Nplanets);
    301240  gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets);
    … …  
    305244  free (Rref);
    306245  free (Dref);
    307   free (Rvel);
    308   free (Dvel);
    309246  free (dist);
    310247  free (Mref);
    … …  
    574511}
    575512
    576 void mpcorb_testing (int argc, char **argv) {
    577 
    578   if (argc != 3) Shutdown ("USAGE: %s (MPCORB.DAT) (epoch)\n", argv[0]);
    579 
    580   char  *filename = argv[1];
    581   double mjdObs   = atof (argv[2]);
    582 
    583   int Nplanets = 0;
    584   Planets *planets = mpcorb_read_text (filename, &Nplanets);
    585   if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
    586 
    587   if (TESTING == 1) { mpcorb_predict_test(planets, Nplanets); exit (0); }
    588 
    589   if (TESTING == 2) {
    590     for (int i = 0; i < 10; i++) {
    591       PlanetDatum myPlanet = mpcorb_predict (&planets[i], mjdObs, TRUE);
    592       char ra_string[64], de_string[64];
    593       hms_format (ra_string, 64, myPlanet.Robs/15.0);
    594       hms_format (de_string, 64, myPlanet.Dobs);
    595       fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", myPlanet.Robs, myPlanet.Dobs, ra_string, de_string);
    596     }
    597     exit (0);
    598   }
    599   fprintf (stderr, "ERROR: invalid value of TESTING : %d\n", TESTING);
    600   exit (2);
    601 }
    602    
    603 // Supply this test with the full list of minor planets from MPCORB.DAT.
    604 // It will compare the positions for the first 10 minor planets with the
    605 // reference data downloaded from the MPC for:
    606 // MJD = 59396.0 @ (Long,Lat,Alt) = (0.0, 0.0, 0.0)
    607 
    608 void mpcorb_predict_test (Planets *planets, int Nplanets) {
    609  
    610   int    Nref = 10;
    611   double Rref[10] = { 49.14708333, 352.90125000, 251.95416667, 175.74666667, 266.25541667, 298.38000000,  52.84458333, 121.94916667, 192.29500000, 144.20583333 };
    612   double Dref[10] = { 11.48944444,   8.68722222,  -4.19666667,   9.37111111, -16.97583333,  -8.30500000,  22.13472222,  21.50888889,  -1.52000000,  11.57138889 };
    613 
    614   for (int i = 0; (i < Nref) && (i < Nplanets); i++) {
    615 
    616     if (VERBOSE) fprintf (stderr, "%8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.7f %12.7f\n",
    617              planets[i].ID,
    618              planets[i].epoch,
    619              planets[i].mean_anomaly,
    620              planets[i].perihelion,
    621              planets[i].ascend_node,
    622              planets[i].inclination,
    623              planets[i].eccentricity,
    624              planets[i].semimajor_a);
    625    
    626     // generate predictions for (Long,Lat,Alt) = (0.0, 0.0, 0.0)
    627     double Rpre, Dpre;
    628     mpcorb_predict_test_calc (&planets[i], 59396.0, &Rpre, &Dpre);
    629 
    630     fprintf (stderr, "result: %8s : %12.6f %12.6f : %6.2f %6.2f\n", planets[i].ID, Rpre, Dpre, 3600*(Rref[i]-Rpre)*cos(Dref[i]*RAD_DEG), 3600*(Dref[i]-Dpre));
    631   }
    632 }
    633 
    634 // attempt to predict as a test for Long,Lat = 0,0, with all other corrections
    635 // this attempt uses slaPertel to calculate perturbed versions of the orbital elements
    636 // NOTE this test seems to have worked: the calculated positions of the first 10 asteroids reported
    637 // by MPC match to within 2 arcsec.  This is much larger than the Planetary abberation,
    638 // so slaPlante must already be accounting for that effect.
    639 
    640 void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
     513// mjdObs is NOT in TT
     514PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc) {
    641515
    642516  int jstat;
    … …  
    646520
    647521  Planets tmpPlanet;
    648 
    649   // correct the orbital elements to the perturbed version
    650   slaPertel (2,
    651              planet->epoch,
    652              mjdObsTT,
    653              planet->epoch,
    654              planet->inclination*RAD_DEG,
    655              planet->ascend_node*RAD_DEG,
    656              planet->perihelion*RAD_DEG,
    657              planet->semimajor_a,
    658              planet->eccentricity,
    659              planet->mean_anomaly*RAD_DEG,
    660              &tmpPlanet.epoch,
    661              &tmpPlanet.inclination,  // returned as radians and used as radians below
    662              &tmpPlanet.ascend_node,  // returned as radians and used as radians below
    663              &tmpPlanet.perihelion,   // returned as radians and used as radians below
    664              &tmpPlanet.semimajor_a,
    665              &tmpPlanet.eccentricity,
    666              &tmpPlanet.mean_anomaly, // returned as radians and used as radians below
    667              &jstat);
    668 
    669   // this function returns the topocentric apparent R,D.  I assume it
    670   // has no atm corrections.
    671 
    672   double Rtopo, Dtopo, dist;
    673   slaPlante (mjdObsTT,
    674              0.0, // PS1_LONGITUDE*RAD_DEG,
    675              0.0, // PS1_LATITUDE*RAD_DEG,
    676              2,
    677              tmpPlanet.epoch,
    678              tmpPlanet.inclination,
    679              tmpPlanet.ascend_node,
    680              tmpPlanet.perihelion,
    681              tmpPlanet.semimajor_a,
    682              tmpPlanet.eccentricity,
    683              tmpPlanet.mean_anomaly,
    684              0.0, &Rtopo, &Dtopo, &dist, &jstat);
    685 
    686   if (VERBOSE) fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rtopo*DEG_RAD, Dtopo*DEG_RAD);
    687 
    688   // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
    689   // Setting the pressure to 0.0 seems to avoid those corrections.
    690 
    691   double Rap, Dap;
    692   slaOap ("R",
    693           Rtopo,
    694           Dtopo,
    695           mjdObs,
    696           -0.4,
    697           0.0, // LONGITUDE
    698           0.0, // LATITUDE
    699           0.0, // ALTITUDE,
    700           0.0, // xp
    701           0.0, // yp
    702           273.15, // T(K)
    703           0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
    704           0.0, // RH(%)
    705           0.7, // Wavelength (um)
    706           0.0, // lapse rate
    707           &Rap, &Dap);
    708 
    709   if (VERBOSE) fprintf (stderr, "Oap   : %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rtopo)*DEG_RAD*cos(Dtopo), 3600*(Dap-Dtopo)*DEG_RAD);
    710 
    711   double Rmean, Dmean;
    712   slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
    713 
    714   if (VERBOSE) fprintf (stderr, "Amp   : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
    715 
    716   if (0) {
    717     // slaAmp applies both precession/nutation as well as abberation
    718     // but, I think abberation is not appropriate to solar system bodies.
    719 
    720     // here, I instead calculate just the precession and nutation with slalib
    721 
    722     double Rpre, Dpre;
    723     Rpre = Rap;
    724     Dpre = Dap;
    725     slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
    726     if (VERBOSE) fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
    727   }
    728 
    729   *Robs = Rmean * DEG_RAD;
    730   *Dobs = Dmean * DEG_RAD;
    731 }
    732 
    733 void myAmp ( double ra, double da, double date, double eq, double *rm, double *dm );
    734 
    735 // mjdObs is NOT in TT
    736 PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc) {
    737 
    738   int jstat;
    739  
    740   double dTT = slaDtt (mjdObs);
    741   double mjdObsTT = mjdObs + dTT / 86400.0;
    742 
    743   Planets tmpPlanet;
    744522  PlanetDatum myPlanet;
    745523  myPlanet.Robs = NAN;
    746524  myPlanet.Dobs = NAN;
    747   myPlanet.Rvel = NAN;
    748   myPlanet.Dvel = NAN;
    749525  myPlanet.dist = NAN;
    750526
    751527  // correct the orbital elements to the perturbed version
    752   tmpPlanet = *planet;
    753528  if (FullCalc) { slaPertel (2,
    754529             planet->epoch,
    … …  
    780555  }
    781556
    782   // this function returns the topocentric apparent R,D.  I assume it
    783   // has no atm corrections.
     557  // this function returns the topocentric apparent R,D.  It does no atm corrections.
    784558
    785559  double Rtopo, Dtopo, dist;
    … …  
    798572  myPlanet.dist = dist;
    799573
    800   double dRobs = 0.0; // planetary aberration, if FullCalc
    801   double dDobs = 0.0; // planetary aberration, if FullCalc
    802   if (FullCalc) {
    803     // calculate the velocity components
    804    
    805     // find position one minute later
    806     double Rtop2, Dtop2, dist2;
    807     slaPlante (mjdObsTT + 60/86400.0,
    808                PS1_LONGITUDE*RAD_DEG,
    809                PS1_LATITUDE*RAD_DEG,
    810                2,
    811                tmpPlanet.epoch,
    812                tmpPlanet.inclination,
    813                tmpPlanet.ascend_node,
    814                tmpPlanet.perihelion,
    815                tmpPlanet.semimajor_a,
    816                tmpPlanet.eccentricity,
    817                tmpPlanet.mean_anomaly,
    818                0.0, &Rtop2, &Dtop2, &dist2, &jstat);
    819 
    820     myPlanet.Rvel = 3600.0*(Rtop2 - Rtopo)*DEG_RAD*cos(Dtopo); // this vector is wrong by the precession rotation
    821     myPlanet.Dvel = 3600.0*(Dtop2 - Dtopo)*DEG_RAD;            // this vector is wrong by the precession rotation
    822 
    823     // include the correction for planetary aberration
    824     // XXX make these constants more accurate
    825     double dist_km = dist * 1.5e8;
    826     double time_min = dist_km / 299792.459/ 60.0;
    827     double Rlag = myPlanet.Rvel * time_min; // velocity is in arcsec / min
    828     double Dlag = myPlanet.Dvel * time_min; // velocity is in arcsec / min
    829    
    830     dRobs = Rlag / 3600.0 / cos(Dtopo); // RA degrees
    831     dDobs = Dlag / 3600.0;              // DEC degrees
    832   }
    833 
    834574  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
    835575  // Setting the pressure to 0.0 seems to avoid those corrections.
     576  // Should I set altitude to zero?
    836577
    837578  double Rap, Dap;
    … …  
    857598  double Rmean, Dmean;
    858599  if (FullCalc) {
    859     if (USE_AMP) {
    860       // slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
    861       myAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
    862     } else {
    863       // use slaPreces and do not correct for stellar aberration
    864       // convert mjdObsTT to year
    865       double epochYear = 2000.0 + (mjdObsTT - 51544.5) / 365.25; // J2000 = 51544.5
    866       Rmean = Rap;
    867       Dmean = Dap;
    868 
    869       // use slalib precession and nutation:
    870       if (USE_MANUAL_PRECESS) {
    871         double pm[3][3], v1[3], v2[3];
    872 
    873         /* Convert RA,Dec to x,y,z */
    874         slaDcs2c ( Rmean, Dmean, v1 );
    875 
    876         // calculate the precession & nutation combined matrix:
    877         if (USE_PRENUT) {
    878           // slaPrenut ( epochYear, 51544.5, pm );
    879           slaPrenut ( 2000.0, mjdObsTT, pm ); // prenut is always mean-to-apparent
    880           slaDimxv ( pm, v1, v2 );            // apply the inverse transformation
    881         } else {
    882           slaPrec ( epochYear, 2000.0, pm );  // precess from apparent-to-mean
    883           slaDmxv ( pm, v1, v2 );             // apply the forward transformatio
    884         }
    885 
    886         /* Back to RA,Dec */
    887         slaDcc2s ( v2, &Rmean, &Dmean );
    888 
    889         // renorm to 0-360
    890         Rmean = slaDranrm ( Rmean );
    891       } else {
    892         slaPreces ("FK5", epochYear, 2000.000000, &Rmean, &Dmean);
    893       }
    894       // fprintf (stderr, "year: %f\n", epochYear);
    895     }
     600    slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
    896601  } else {Rmean = Rap; Dmean = Dap; }
    897602
    898   if (USE_PLANETARY_ABERRATION) {
    899     myPlanet.Robs = Rmean * DEG_RAD + dRobs;
    900     myPlanet.Dobs = Dmean * DEG_RAD + dDobs;
    901   } else {
    902     myPlanet.Robs = Rmean * DEG_RAD;
    903     myPlanet.Dobs = Dmean * DEG_RAD;
    904   }
     603  myPlanet.Robs = Rmean * DEG_RAD;
     604  myPlanet.Dobs = Dmean * DEG_RAD;
    905605  return myPlanet;
    906 }
    907 
    908 # define DISABLE_ABERRATION FALSE
    909 # define DISABLE_DEFLECTION FALSE
    910 
    911 void myAmpqk ( double ra, double da, double amprms[21], double *rm, double *dm ) {
    912    double gr2e;    /* (grav rad Sun)*2/(Sun-Earth distance) */
    913    double ab1;     /* sqrt(1-v*v) where v=modulus of Earth vel */
    914    double ehn[3];  /* Earth position wrt Sun (unit vector, FK5) */
    915    double abv[3];  /* Earth velocity wrt SSB (c, FK5) */
    916    double p[3], p1[3], p2[3], p3[3];  /* work vectors */
    917    double ab1p1, p1dv, p1dvp1, w, pde, pdep1;
    918    int i, j;
    919 
    920 /* Unpack some of the parameters */
    921    gr2e = amprms[7];
    922    ab1  = amprms[11];
    923    for ( i = 0; i < 3; i++ ) {
    924       ehn[i] = amprms[i + 4];
    925       abv[i] = amprms[i + 8];
    926    }
    927 
    928    // invert and print out:
    929    if (VERBOSE_TEST) {
    930      double Rmean, Dmean;
    931      Rmean = ra * DEG_RAD;
    932      Dmean = da * DEG_RAD;
    933      fprintf (stderr, "%12.6f %12.6f ", Rmean, Dmean);
    934    }
    935 
    936 /* Apparent RA,Dec to Cartesian */
    937    slaDcs2c ( ra, da, p3 );
    938 
    939    // invert and print out:
    940    if (VERBOSE_TEST) {
    941      double Rtmp, Dtmp, Rmean, Dmean;
    942      slaDcc2s ( p3, &Rtmp, &Dtmp );
    943      Rtmp = slaDranrm ( Rtmp );
    944      Rmean = Rtmp * DEG_RAD;
    945      Dmean = Dtmp * DEG_RAD;
    946      fprintf (stderr, "C: %12.6f %12.6f ", Rmean, Dmean);
    947    }
    948 
    949 /* Precession and nutation */
    950    slaDimxv ( (double(*)[3]) &amprms[12], p3, p2 );
    951 
    952    // invert and print out:
    953    if (VERBOSE_TEST) {
    954      double Rtmp, Dtmp, Rmean, Dmean;
    955      slaDcc2s ( p2, &Rtmp, &Dtmp );
    956      Rtmp = slaDranrm ( Rtmp );
    957      Rmean = Rtmp * DEG_RAD;
    958      Dmean = Dtmp * DEG_RAD;
    959      fprintf (stderr, "P: %12.6f %12.6f ", Rmean, Dmean);
    960    }
    961 
    962    if (DISABLE_ABERRATION) {
    963      // XXX : save result from p2 in p1
    964      for ( i = 0; i < 3; i++ ) {
    965        p1[i] = p2[i];
    966      }
    967    } else {
    968      /* Aberration */
    969      ab1p1 = ab1 + 1.0;
    970      for ( i = 0; i < 3; i++ ) {
    971        p1[i] = p2[i];
    972      }
    973      for ( j = 0; j < 2; j++ ) {
    974        p1dv = slaDvdv ( p1, abv );
    975        p1dvp1 = 1.0 + p1dv;
    976        w = 1.0 + p1dv / ab1p1;
    977        for ( i = 0; i < 3; i++ ) {
    978          p1[i] = ( p1dvp1 * p2[i] - w * abv[i] ) / ab1;
    979        }
    980        slaDvn ( p1, p3, &w );
    981        for ( i = 0; i < 3; i++ ) {
    982          p1[i] = p3[i];
    983        }
    984      }
    985    }
    986 
    987    // invert and print out:
    988    if (VERBOSE_TEST) {
    989      double Rtmp, Dtmp, Rmean, Dmean;
    990      slaDcc2s ( p1, &Rtmp, &Dtmp );
    991      Rtmp = slaDranrm ( Rtmp );
    992      Rmean = Rtmp * DEG_RAD;
    993      Dmean = Dtmp * DEG_RAD;
    994      fprintf (stderr, "A: %12.6f %12.6f ", Rmean, Dmean);
    995    }
    996 
    997 /* Light deflection */
    998    for ( i = 0; i < 3; i++ ) {
    999       p[i] = p1[i];
    1000    }
    1001 
    1002 // disable deflection:
    1003    if (DISABLE_DEFLECTION) {
    1004      
    1005    } else {
    1006      for ( j = 0; j < 5; j++ ) {
    1007        pde = slaDvdv ( p, ehn );
    1008        pdep1 = 1.0 + pde;
    1009        w = pdep1 - gr2e * pde;
    1010        for ( i = 0; i < 3; i++ ) {
    1011          p[i] = ( pdep1 * p1[i] - gr2e * ehn[i] ) / w;
    1012        }
    1013        slaDvn ( p, p2, &w );
    1014        for ( i = 0; i < 3; i++ ) {
    1015          p[i] = p2[i];
    1016        }
    1017      }
    1018    }
    1019 
    1020    // invert and print out:
    1021    if (VERBOSE_TEST) {
    1022      double Rtmp, Dtmp, Rmean, Dmean;
    1023      slaDcc2s ( p, &Rtmp, &Dtmp );
    1024      Rtmp = slaDranrm ( Rtmp );
    1025      Rmean = Rtmp * DEG_RAD;
    1026      Dmean = Dtmp * DEG_RAD;
    1027      fprintf (stderr, "D: %12.6f %12.6f\n", Rmean, Dmean);
    1028    }
    1029 
    1030 /* Mean RA,Dec */
    1031    slaDcc2s ( p, rm, dm );
    1032    *rm = slaDranrm ( *rm );
    1033 }
    1034 
    1035 void myAmp ( double ra, double da, double date, double eq, double *rm, double *dm ) {
    1036 
    1037    double amprms[21];    /* Mean-to-apparent parameters */
    1038 
    1039    slaMappa ( eq, date, amprms );
    1040     myAmpqk ( ra, da, amprms, rm, dm );
    1041606}
    1042607
    … …  
    1181746}
    1182747
    1183 // The following test case provided by Serge can be perfectly reproduced, but only if I
    1184 // use Dtt for the epoch, and this is clearly wrong.
    1185 void mpcorb_predict_test_from_serge () {
    1186 
    1187   double epochUT = 57435.5;
    1188 
    1189   Planets testPlanet;
    1190   // testPlanet.epoch        = epochUT + slaDtt(epochUT) / 86400.0; // XXX this is what the epoch should be
    1191   testPlanet.epoch        = slaDtt(epochUT); // XXX this is the value that works with the test values
    1192   testPlanet.inclination  = 0.043306854729735354 * DEG_RAD;
    1193   testPlanet.ascend_node  = 4.728915154005972    * DEG_RAD;
    1194   testPlanet.perihelion   = 2.2756608097151414   * DEG_RAD;
    1195   testPlanet.mean_anomaly = 2.888163341284418    * DEG_RAD;
    1196   testPlanet.eccentricity = 0.09528339999999975;
    1197   testPlanet.semimajor_a  = 2.782049599999998;
    1198                
    1199   PlanetDatum myPlanet = mpcorb_predict (&testPlanet, epochUT, TRUE);
    1200 
    1201   fprintf (stderr, "%f vs %f : %f\n",  3.92105439135726600*DEG_RAD, myPlanet.Robs, 3600*(3.92105439135726600*DEG_RAD - myPlanet.Robs));
    1202   fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, myPlanet.Dobs, 3600*(-0.33499850519808244*DEG_RAD - myPlanet.Dobs));
    1203   exit (0);
    1204 }
    1205 
    1206 // The following test code uses unperturbed orbital elements and gets the
    1207 // wrong answer by 20 - 60 arcsec for the top 10 minor planets.
    1208 void mpcorb_predict_t0 (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
    1209 
    1210   int jstat;
    1211  
    1212   double dTT = slaDtt (mjdObs);
    1213   double mjdObsTT = mjdObs + dTT / 86400.0;
    1214 
    1215   // this function returns the topocentric apparent R,D.  I assume it
    1216   // has no atm corrections.
    1217 
    1218   double Rrad, Drad, dist;
    1219   slaPlante (mjdObsTT,
    1220              0.0, // PS1_LONGITUDE*RAD_DEG,
    1221              0.0, // PS1_LATITUDE*RAD_DEG,
    1222              2,
    1223              planet->epoch,
    1224              planet->inclination*RAD_DEG,
    1225              planet->ascend_node*RAD_DEG,
    1226              planet->perihelion*RAD_DEG,
    1227              planet->semimajor_a,
    1228              planet->eccentricity,
    1229              planet->mean_anomaly*RAD_DEG,
    1230              0.0, &Rrad, &Drad, &dist, &jstat);
    1231 
    1232   fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rrad*DEG_RAD, Drad*DEG_RAD);
    1233 
    1234   // I'm not certain if I should use realistic (STD) atm parameters or not.
    1235   // If I set them to zero, is the correction not applied?
    1236  
    1237   double Rap, Dap;
    1238   slaOap ("R",
    1239           Rrad,
    1240           Drad,
    1241           mjdObs,
    1242           -0.4,
    1243           0.0, // LONGITUDE
    1244           0.0, // LATITUDE
    1245           0.0, // ALTITUDE,
    1246           0.0, // xp
    1247           0.0, // yp
    1248           273.15, // T(K)
    1249           1013.25, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
    1250           0.0, // RH(%)
    1251           0.7, // Wavelength (um)
    1252           0.0065, // lapse rate
    1253           &Rap, &Dap);
    1254 
    1255   // XXX the atm effects seem to be large and incompatible with MPC outputs
    1256   // fprintf (stderr, "Oap (STD): %12.6f %12.6f\n", Rap*DEG_RAD, Dap*DEG_RAD);
    1257 
    1258   slaOap ("R",
    1259           Rrad,
    1260           Drad,
    1261           mjdObs,
    1262           -0.4,
    1263           0.0, // LONGITUDE
    1264           0.0, // LATITUDE
    1265           0.0, // ALTITUDE,
    1266           0.0, // xp
    1267           0.0, // yp
    1268           273.15, // T(K)
    1269           0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
    1270           0.0, // RH(%)
    1271           0.7, // Wavelength (um)
    1272           0.0, // lapse rate
    1273           &Rap, &Dap);
    1274 
    1275   fprintf (stderr, "Oap (0.0): %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rrad)*DEG_RAD*cos(Drad), 3600*(Dap-Drad)*DEG_RAD);
    1276 
    1277   double Rmean, Dmean;
    1278   slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
    1279 
    1280   fprintf (stderr, "Amp1     : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
    1281 
    1282   // apply Amp to output from Plante
    1283   double Ralt, Dalt;
    1284   slaAmp (Rrad, Drad, mjdObsTT, 2000.0, &Ralt, &Dalt);
    1285 
    1286   fprintf (stderr, "Amp2     : %12.6f %12.6f : %6.2f %6.2f\n", Ralt*DEG_RAD, Dalt*DEG_RAD, 3600*(Ralt-Rrad)*DEG_RAD*cos(Drad), 3600*(Dalt-Drad)*DEG_RAD);
    1287 
    1288 
    1289   // slaAmp applies both precession/nutation as well as abberation
    1290   // but, I think abberation is not appropriate to solar system bodies.
    1291 
    1292   // here, I instead calculate just the precession and nutation with slalib
    1293 
    1294   double Rpre, Dpre;
    1295   Rpre = Rap;
    1296   Dpre = Dap;
    1297   slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
    1298   fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
    1299 
    1300   *Robs = Rmean * DEG_RAD;
    1301   *Dobs = Dmean * DEG_RAD;
    1302 }
    1303 
Note: See TracChangeset for help on using the changeset viewer.