IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41711


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

adding distance to calculate apparent mags, fix perturbed orbits

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r41710 r41711  
    3333} Planets;
    3434
     35typedef struct {
     36  char   ID[8];
     37  double Robs;
     38  double Dobs;
     39  double Rvel;
     40  double Dvel;
     41  double dist;
     42} PlanetDatum;
     43
    3544int Shutdown (char *format, ...);
    3645int mpcorb_parseline (char *line, Planets *planet, int Nmax);
    … …  
    3847Planets *mpcorb_read_fits (char *filename, int *nplanets);
    3948void     mpcorb_save_fits (char *filename, Planets *planets, int Nplanets);
    40 void     mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets);
    41 
    42 void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc);
     49void     mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets);
     50
     51PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc);
    4352
    4453void mpcorb_trange (int argc, char **argv);
    … …  
    5463# define PS1_ALTITUDE  3067.7 /* meters */
    5564
     65// some options (for testing)
     66int USE_AMP                  = FALSE; // otherwise use precess and/or nutation
     67int USE_PLANETARY_ABERRATION = FALSE;
     68int USE_MANUAL_PRECESS       = FALSE; // otherwise use slaPreces
     69int USE_PRENUT               = FALSE; // otherwise slaPrec
     70
    5671int main (int argc, char **argv) {
    5772
     73  int N;
     74
    5875  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  }
    59102
    60103  if ((argc != 9) && (argc != 10)) goto usage;
    … …  
    100143
    101144    // is the object in the region for the start of the period?
    102     double RatMin, DatMin;
    103     mpcorb_predict (&planets[i], mjdMin, &RatMin, &DatMin, FALSE);
     145    PlanetDatum plMin = mpcorb_predict (&planets[i], mjdMin, FALSE);
    104146    int inRangeMin = TRUE;
    105     inRangeMin = inRangeMin && (RatMin > Rmin);
    106     inRangeMin = inRangeMin && (RatMin < Rmax);
    107     inRangeMin = inRangeMin && (DatMin > Dmin);
    108     inRangeMin = inRangeMin && (DatMin < Dmax);
     147    inRangeMin = inRangeMin && (plMin.Robs > Rmin);
     148    inRangeMin = inRangeMin && (plMin.Robs < Rmax);
     149    inRangeMin = inRangeMin && (plMin.Dobs > Dmin);
     150    inRangeMin = inRangeMin && (plMin.Dobs < Dmax);
    109151   
    110152    // is the object in the region for the end of the period?
    111     double RatMax, DatMax;
    112     mpcorb_predict (&planets[i], mjdMax, &RatMax, &DatMax, FALSE);
     153    PlanetDatum plMax = mpcorb_predict (&planets[i], mjdMax, FALSE);
    113154    int inRangeMax = TRUE;
    114     inRangeMax = inRangeMax && (RatMax > Rmin);
    115     inRangeMax = inRangeMax && (RatMax < Rmax);
    116     inRangeMax = inRangeMax && (DatMax > Dmin);
    117     inRangeMax = inRangeMax && (DatMax < Dmax);
     155    inRangeMax = inRangeMax && (plMax.Robs > Rmin);
     156    inRangeMax = inRangeMax && (plMax.Robs < Rmax);
     157    inRangeMax = inRangeMax && (plMax.Dobs > Dmin);
     158    inRangeMax = inRangeMax && (plMax.Dobs < Dmax);
    118159
    119160    // NOTE: if an object is moving so fast that it traverses the entire region, we will
    … …  
    123164    planetsSave[NplanetsSave] = planets[i];
    124165    planetsSave[NplanetsSave].mjdMin = mjdMin;
    125     planetsSave[NplanetsSave].RatMin = RatMin;
    126     planetsSave[NplanetsSave].DatMin = DatMin;
     166    planetsSave[NplanetsSave].RatMin = plMin.Robs;
     167    planetsSave[NplanetsSave].DatMin = plMin.Dobs;
    127168    planetsSave[NplanetsSave].mjdMax = mjdMax;
    128     planetsSave[NplanetsSave].RatMax = RatMax;
    129     planetsSave[NplanetsSave].DatMax = DatMax;
     169    planetsSave[NplanetsSave].RatMax = plMax.Robs;
     170    planetsSave[NplanetsSave].DatMax = plMax.Dobs;
    130171   
    131172    NplanetsSave++;
    … …  
    134175  fprintf (stderr, "\n");
    135176
     177  // NOTE : here we are saving the full orbital elements
    136178  mpcorb_save_fits (output, planetsSave, NplanetsSave);
    137179
    138180  exit (0);
    139181}
     182
     183# define VERBOSE_TEST FALSE
    140184
    141185// generate a table of asteroid positions which fall within the region at Tobs (in MJD)
    … …  
    152196  char  *output   = argv[8];
    153197
    154   // FILE *fout = fopen (output, "w");
    155   // if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output);
    156 
    157198  int Nplanets = 0;
    158199  Planets *planets = mpcorb_read_fits (filename, &Nplanets);
    … …  
    161202  int NplanetsSave = 0;
    162203  int NPLANETSSAVE = 1000;
    163   ALLOCATE_PTR (planetsSave, Planets, NPLANETSSAVE);
     204  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  }
    164213
    165214  // transform all objects and identify those in the target region:
    … …  
    169218    // is the object in the region for the start of the period?
    170219    // first, use the fast, but inaccurate, calculation
    171     double Robs, Dobs;
    172     mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, FALSE);
    173     if (Robs < Rmin) continue;
    174     if (Robs > Rmax) continue;
    175     if (Dobs < Dmin) continue;
    176     if (Dobs > Dmax) continue;
     220    PlanetDatum planetFast = mpcorb_predict (&planets[i], mjdObs, FALSE);
     221    if (planetFast.Robs < Rmin) continue;
     222    if (planetFast.Robs > Rmax) continue;
     223    if (planetFast.Dobs < Dmin) continue;
     224    if (planetFast.Dobs > Dmax) continue;
    177225   
    178226    // if it is in the region, use the more accurate calculation
    179     mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
     227    PlanetDatum planetSlow = mpcorb_predict (&planets[i], mjdObs, TRUE);
    180228
    181229    // fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
    182230
    183231    // save this prediction, with R,D in both @min & @max
    184     planetsSave[NplanetsSave] = planets[i];
    185     planetsSave[NplanetsSave].mjdMin = mjdObs;
    186     planetsSave[NplanetsSave].RatMin = Robs;
    187     planetsSave[NplanetsSave].DatMin = Dobs;
    188     planetsSave[NplanetsSave].mjdMax = mjdObs;
    189     planetsSave[NplanetsSave].RatMax = Robs;
    190     planetsSave[NplanetsSave].DatMax = Dobs;
     232    planetsSave[NplanetsSave] = planetSlow;
     233    strcpy (planetsSave[NplanetsSave].ID, planets[i].ID); // ID is not passed to mpcorb_predict
    191234   
    192235    NplanetsSave++;
    193     CHECK_REALLOCATE (planetsSave, Planets, NPLANETSSAVE, NplanetsSave, 1000);
     236    CHECK_REALLOCATE (planetsSave, PlanetDatum, NPLANETSSAVE, NplanetsSave, 1000);
    194237  }
    195238  fprintf (stderr, "\n");
    … …  
    204247
    205248// write minor planet positions (Rref,Dref,Mref) to a FITS table
    206 void mpcorb_refs_fits (char *filename, Planets *planets, int Nplanets) {
     249void mpcorb_refs_fits (char *filename, PlanetDatum *planets, int Nplanets) {
    207250
    208251  Header header;
    … …  
    218261  gfits_create_table_header (&theader, "BINTABLE", "DATA");
    219262
    220   gfits_define_bintable_column (&theader, "8A", "ID",   "name", "none",        1.0, 0.0);
    221   gfits_define_bintable_column (&theader,  "D", "Rref", "RA",   "degrees",     1.0, 0.0);
    222   gfits_define_bintable_column (&theader,  "D", "Dref", "DEC",  "degrees",     1.0, 0.0);
    223   gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",  "magnitudes",  1.0, 0.0);
     263  gfits_define_bintable_column (&theader, "8A", "ID",   "name",    "none",        1.0, 0.0);
     264  gfits_define_bintable_column (&theader,  "D", "Rref", "RA",      "degrees",     1.0, 0.0);
     265  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);
     268  gfits_define_bintable_column (&theader,  "D", "dist", "DIST",    "AU",          1.0, 0.0);
     269  gfits_define_bintable_column (&theader,  "D", "Mref", "Mag",     "magnitudes",  1.0, 0.0);
    224270
    225271  // generate the output array that carries the data
    … …  
    230276  ALLOCATE_PTR (Rref,                 double, Nplanets);
    231277  ALLOCATE_PTR (Dref,                 double, Nplanets);
     278  ALLOCATE_PTR (Rvel,                 double, Nplanets);
     279  ALLOCATE_PTR (Dvel,                 double, Nplanets);
     280  ALLOCATE_PTR (dist,                 double, Nplanets);
    232281  ALLOCATE_PTR (Mref,                 double, Nplanets);
    233282
    234283  // set intermediate storage arrays
    235284  for (int i = 0; i < Nplanets; i++) {
    236     Rref[i]        = planets[i].RatMin;
    237     Dref[i]        = planets[i].DatMin;
     285    Rref[i]        = planets[i].Robs;
     286    Dref[i]        = planets[i].Dobs;
     287    Rvel[i]        = planets[i].Rvel;
     288    Dvel[i]        = planets[i].Dvel;
     289    dist[i]        = planets[i].dist;
    238290    Mref[i]        = 16.0; // need to add this to prediction
    239291    memcpy(&ID[i*8], planets[i].ID, 8);
    … …  
    244296  gfits_set_bintable_column (&theader, &ftable, "Rref", Rref, Nplanets);
    245297  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);
     300  gfits_set_bintable_column (&theader, &ftable, "dist", dist, Nplanets);
    246301  gfits_set_bintable_column (&theader, &ftable, "Mref", Mref, Nplanets);
    247302
    … …  
    250305  free (Rref);
    251306  free (Dref);
     307  free (Rvel);
     308  free (Dvel);
     309  free (dist);
    252310  free (Mref);
    253311
    … …  
    530588
    531589  if (TESTING == 2) {
    532     double Robs, Dobs;
    533590    for (int i = 0; i < 10; i++) {
    534       mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
     591      PlanetDatum myPlanet = mpcorb_predict (&planets[i], mjdObs, TRUE);
    535592      char ra_string[64], de_string[64];
    536       hms_format (ra_string, 64, Robs/15.0);
    537       hms_format (de_string, 64, Dobs);
    538       fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", Robs, Dobs, ra_string, de_string);
     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);
    539596    }
    540597    exit (0);
    … …  
    674731}
    675732
     733void myAmp ( double ra, double da, double date, double eq, double *rm, double *dm );
     734
    676735// mjdObs is NOT in TT
    677 void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc) {
     736PlanetDatum mpcorb_predict (Planets *planet, double mjdObs, int FullCalc) {
    678737
    679738  int jstat;
    … …  
    683742
    684743  Planets tmpPlanet;
     744  PlanetDatum myPlanet;
     745  myPlanet.Robs = NAN;
     746  myPlanet.Dobs = NAN;
     747  myPlanet.Rvel = NAN;
     748  myPlanet.Dvel = NAN;
     749  myPlanet.dist = NAN;
    685750
    686751  // correct the orbital elements to the perturbed version
     752  tmpPlanet = *planet;
    687753  if (FullCalc) { slaPertel (2,
    688754             planet->epoch,
    … …  
    704770             &jstat);
    705771  } else {
    706     tmpPlanet = *planet;
     772    // the function calls below require elements in radians
     773    tmpPlanet.epoch        = planet->epoch;
     774    tmpPlanet.inclination  = planet->inclination *RAD_DEG;
     775    tmpPlanet.ascend_node  = planet->ascend_node *RAD_DEG;
     776    tmpPlanet.perihelion   = planet->perihelion  *RAD_DEG;
     777    tmpPlanet.semimajor_a  = planet->semimajor_a ;
     778    tmpPlanet.eccentricity = planet->eccentricity;
     779    tmpPlanet.mean_anomaly = planet->mean_anomaly*RAD_DEG;
    707780  }
    708781
    … …  
    715788             PS1_LATITUDE*RAD_DEG,
    716789             2,
    717              planet->epoch,
    718              planet->inclination*RAD_DEG,
    719              planet->ascend_node*RAD_DEG,
    720              planet->perihelion*RAD_DEG,
    721              planet->semimajor_a,
    722              planet->eccentricity,
    723              planet->mean_anomaly*RAD_DEG,
     790             tmpPlanet.epoch,
     791             tmpPlanet.inclination,
     792             tmpPlanet.ascend_node,
     793             tmpPlanet.perihelion,
     794             tmpPlanet.semimajor_a,
     795             tmpPlanet.eccentricity,
     796             tmpPlanet.mean_anomaly,
    724797             0.0, &Rtopo, &Dtopo, &dist, &jstat);
     798  myPlanet.dist = dist;
     799
     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  }
    725833
    726834  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
    … …  
    749857  double Rmean, Dmean;
    750858  if (FullCalc) {
    751     slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
     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    }
    752896  } else {Rmean = Rap; Dmean = Dap; }
    753897
    754   *Robs = Rmean * DEG_RAD;
    755   *Dobs = Dmean * DEG_RAD;
     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  }
     905  return myPlanet;
     906}
     907
     908# define DISABLE_ABERRATION FALSE
     909# define DISABLE_DEFLECTION FALSE
     910
     911void 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
     1035void 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 );
    7561041}
    7571042
    … …  
    7621047
    7631048  // fields I need:
    764 
     1049 
    7651050  // epoch       : 21 - 25
    7661051  // inclination : 60 - 68
    … …  
    9121197  testPlanet.semimajor_a  = 2.782049599999998;
    9131198               
    914   double Robs, Dobs;
    915   mpcorb_predict (&testPlanet, epochUT, &Robs, &Dobs, TRUE);
    916 
    917   fprintf (stderr, "%f vs %f : %f\n",  3.92105439135726600*DEG_RAD, Robs, 3600*(3.92105439135726600*DEG_RAD - Robs));
    918   fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, Dobs, 3600*(-0.33499850519808244*DEG_RAD - Dobs));
     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));
    9191203  exit (0);
    9201204}
Note: See TracChangeset for help on using the changeset viewer.