#ifndef SLAAINCLD #define SLAAINCLD #include #include #include using std::string; using std::ostream; using std::endl; extern"C" { void sla_pv2el_ ( const double[6],const double*,const double*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*); } extern"C" { void sla_geoc_ ( const double*,const double*,double*,double*); } extern"C" { float sla_sep_ ( const float*,const float*,const float*,const float*); } extern"C" { void sla_ampqk_ ( const double*,const double*,const double[21],double*,double*); } extern"C" { void sla_aoppat_ ( const double*,double[14]); } extern"C" { double sla_epb2d_ ( const double*); } extern"C" { void sla_dafin_ ( const char*,const int*,int*,double*,int*); } extern"C" { void sla_ecmat_ ( const double*,double[3][3]); } extern"C" { void sla_amp_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_tp2v_ ( const float*,const float*,const float[3],float[3]); } extern"C" { void sla_refz_ ( const double*,const double*,const double*,double*); } extern"C" { void sla_caldj_ ( const int*,const int*,const int*,double*,int*); } extern"C" { double sla_gmst_ ( const double*); } extern"C" { void sla_calyd_ ( const int*,const int*,const int*,int*,int*,int*); } extern"C" { float sla_rvgalc_ ( const float*,const float*); } extern"C" { void sla_mxv_ ( const float[3][3],const float[3],float[3]); } extern"C" { void sla_epv_ ( const double*,double[3],double[3],double[3],double[3]); } extern"C" { double sla_dranrm_ ( const double*); } extern"C" { void sla_dtf2d_ ( const int*,const int*,const double*,double*,int*); } extern"C" { void sla_galsup_ ( const double*,const double*,double*,double*); } extern"C" { void sla_preces_ ( const char*,const double*,const double*,double*,double*); } extern"C" { void sla_euler_ ( const char*,const int*,const float*,const float*,const float*,float[3][3]); } extern"C" { void sla_nutc_ ( const double*,double*,double*,double*); } extern"C" { void sla_refco_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { float sla_range_ ( const double*); } extern"C" { void sla_dvn_ ( const double[3],double[3],double*); } extern"C" { void sla_pertel_ ( const int*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,double*,double*,double*,double*,int*); } extern"C" { void sla_mapqkz_ ( const double*,const double*,const double[21],double*,double*); } extern"C" { void sla_map_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_aoppa_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double[14]); } extern"C" { void sla_polmo_ ( const double*,const double*,double*,const double*,double*,double*,double*); } extern"C" { void sla_dmxv_ ( const double[3][3],const double[3],double[3]); } extern"C" { void sla_deuler_ ( const char*,const int*,const double*,const double*,const double*,double[3][3]); } extern"C" { void sla_cs2c_ ( const float*,const float*,float[3]); } extern"C" { void sla_imxv_ ( const float[3][3],const float[3],float[3]); } extern"C" { void sla_fk5hz_ ( const double*,const double*,const double*,double*,double*); } extern"C" { float sla_vdv_ ( const float[3],const float[3]); } extern"C" { void sla_dcc2s_ ( const double[3],double*,double*); } extern"C" { void sla_refro_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*); } extern"C" { void sla_av2m_ ( const float[3],float[3][3]); } extern"C" { void sla_el2ue_ ( const double*,const int*,double*,double*,double*,double*,double*,double*,double*,double*,double[13],int*); } extern"C" { void sla_cr2af_ ( const int*,const float*,char*,int[4]); } extern"C" { void sla_ecleq_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_mappa_ ( const double*,const double*,double[21]); } extern"C" { void sla_ctf2d_ ( const int*,const int*,const float*,float*,int*); } extern"C" { void sla_unpcd_ ( const double*,double*,double*); } extern"C" { void sla_cd2tf_ ( const int*,const float*,char*,int[4]); } extern"C" { void sla_mxm_ ( const float[3][3],const float[3][3],float[3][3]); } extern"C" { void sla_dfltin_ ( const char*,int*,double*,int*); } extern"C" { void sla_plantu_ ( const double*,const double*,const double*,const double[13],double*,double*,double*,int*); } extern"C" { float sla_pav_ ( const float[3],const float[3]); } extern"C" { void sla_dvxv_ ( const double[3],const double[3],double[3]); } extern"C" { void sla_altaz_ ( const double*,const double*,const double*,double*,double*,double*,double*,double*,double*,double*,double*,double*); } extern"C" { void sla_dav2m_ ( const double[3],double[3][3]); } extern"C" { void sla_de2h_ ( const double*,const double*,const double*,double*,double*); } extern"C" { double sla_airmas_ ( const double*); } extern"C" { void sla_moon_ ( const int*,const int*,const float*,float[6]); } extern"C" { double sla_dvdv_ ( const double[3],const double[3]); } extern"C" { void sla_dd2tf_ ( const int*,const double*,char*,int[4]); } extern"C" { void sla_galeq_ ( const double*,const double*,double*,double*); } extern"C" { void sla_oap_ ( const char*,const int*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_nutc80_ ( const double*,double*,double*,double*); } extern"C" { void sla_planel_ ( const double*,const int*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double[6],int*); } extern"C" { void sla_h2fk5_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { void sla_dmxm_ ( const double[3][3],const double[3][3],double[3][3]); } extern"C" { void sla_vn_ ( const float[3],float[3],float*); } extern"C" { void sla_clyd_ ( const int*,const int*,const int*,int*,int*,int*); } extern"C" { void sla_m2av_ ( const float[3][3],float[3]); } extern"C" { void sla_dr2tf_ ( const int*,const double*,char*,int[4]); } extern"C" { void sla__atms_ ( const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { double sla_pa_ ( const double*,const double*,const double*); } extern"C" { void sla_dcmpf_ ( const double[6],double*,double*,double*,double*,double*,double*); } extern"C" { double sla_dbear_ ( const double*,const double*,const double*,const double*); } extern"C" { void sla_afin_ ( const char*,const int*,int*,float*,int*); } extern"C" { void sla_ds2tp_ ( const double*,const double*,const double*,const double*,double*,double*,int*); } extern"C" { void sla_dtp2s_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_fk52h_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { void sla_mapqk_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double[21],double*,double*); } extern"C" { float sla_rvlsrd_ ( const float*,const float*); } extern"C" { double sla_dt_ ( const double*); } extern"C" { void sla_ctf2r_ ( const int*,const int*,const float*,float*,int*); } extern"C" { double sla_dpav_ ( const double[3],const double[3]); } extern"C" { void sla_pertue_ ( const double*,const double[13],int*); } extern"C" { double sla_epj_ ( const double*); } extern"C" { void sla_atmdsp_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_prebn_ ( const double*,const double*,double[3][3]); } extern"C" { void sla_h2e_ ( const float*,const float*,const float*,float*,float*); } extern"C" { double sla_rcc_ ( const double*,const double*,const double*,const double*,const double*); } extern"C" { void sla_ue2el_ ( const double[13],int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*); } extern"C" { void sla_eg50_ ( const double*,const double*,double*,double*); } extern"C" { void sla_invf_ ( const double[6],double[6],int*); } extern"C" { void sla_fk45z_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_dimxv_ ( const double[3][3],const double[3],double[3]); } extern"C" { void sla_daf2r_ ( const int*,const int*,const double*,double*,int*); } extern"C" { void sla_pdq2h_ ( const double*,const double*,const double*,double*,int*,double*,int*); } extern"C" { void sla_dh2e_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_dtps2c_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*,int*); } extern"C" { void sla_dm2av_ ( const double[3][3],double[3]); } extern"C" { void sla_prec_ ( const double*,const double*,double[3][3]); } extern"C" { float sla_rvlsrk_ ( const float*,const float*); } extern"C" { void sla_plante_ ( const double*,const double*,const double*,const int*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,int*); } extern"C" { double sla_epb_ ( const double*); } extern"C" { void sla_precl_ ( const double*,const double*,double[3][3]); } extern"C" { void sla_kbj_ ( const int*,const double*,char*,int*); } extern"C" { void sla_pv2ue_ ( const double[6],const double*,const double*,double[13],int*); } extern"C" { void sla_dcs2c_ ( const double*,const double*,double[3]); } extern"C" { void sla_fk54z_ ( const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { void sla__idchf_ ( const char*,int*,int*,int*,double*); } extern"C" { void sla_djcl_ ( const double*,int*,int*,int*,double*,int*); } extern"C" { double sla_eqeqx_ ( const double*); } extern"C" { void sla_dc62s_ ( const double[6],double*,double*,double*,double*,double*,double*); } extern"C" { void sla_eqgal_ ( const double*,const double*,double*,double*); } extern"C" { void sla_ge50_ ( const double*,const double*,double*,double*); } extern"C" { double sla_epco_ ( const char*,const char*,const double*); } extern"C" { float sla_rvlg_ ( const float*,const float*); } extern"C" { void sla_cs2c6_ ( const float*,const float*,const float*,const float*,const float*,const float*,float[6]); } extern"C" { void sla_s2tp_ ( const float*,const float*,const float*,const float*,float*,float*,int*); } extern"C" { void sla_cc2s_ ( const float[3],float*,float*); } extern"C" { float sla_ranorm_ ( const double*); } extern"C" { void sla_nut_ ( const double*,double[3][3]); } extern"C" { void sla_v2tp_ ( const float[3],const float[3],float*,float*,int*); } extern"C" { void sla_djcal_ ( int*,const double*,int[4],int*); } extern"C" { void sla_dtpv2c_ ( const double*,const double*,const double[3],double[3],double[3],int*); } extern"C" { void sla__idchi_ ( const char*,int*,int*,double*); } extern"C" { double sla_dtt_ ( const double*); } extern"C" { void sla_ds2c6_ ( const double*,const double*,const double*,const double*,const double*,const double*,double[6]); } extern"C" { void sla_etrms_ ( const double*,double[3]); } extern"C" { void sla_refcoq_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_eqecl_ ( const double*,const double*,const double*,double*,double*); } extern"C" { double sla_gmsta_ ( const double*,const double*); } extern"C" { void sla_tps2c_ ( const float*,const float*,const float*,const float*,float*,float*,float*,float*,int*); } extern"C" { void sla_fk425_ ( const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,double*,double*,double*); } extern"C" { void sla_dv2tp_ ( const double[3],const double[3],double*,double*,int*); } extern"C" { void sla_tp2s_ ( const float*,const float*,const float*,const float*,float*,float*); } extern"C" { void sla_dmoon_ ( const double*,double[6]); } extern"C" { void sla_caf2r_ ( const int*,const int*,const float*,float*,int*); } extern"C" { void sla_earth_ ( const int*,const int*,const float*,float[6]); } extern"C" { void sla_prenut_ ( const double*,const double*,double[3][3]); } extern"C" { double sla_drange_ ( const double*); } extern"C" { void sla_xy2xy_ ( const double*,const double*,const double[6],double*,double*); } extern"C" { float sla_bear_ ( const float*,const float*,const float*,const float*); } extern"C" { void sla_ue2pv_ ( const double*,const double[13],double[6],int*); } extern"C" { void sla_supgal_ ( const double*,const double*,double*,double*); } extern"C" { double sla_epj2d_ ( const double*); } extern"C" { double sla_dat_ ( const double*); } extern"C" { void sla_dr2af_ ( const int*,const double*,char*,int[4]); } extern"C" { double sla_dsepv_ ( const double[3],const double[3]); } extern"C" { void sla_addet_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_refv_ ( const double*,const double*,const double*,double*); } extern"C" { double sla_dsep_ ( const double*,const double*,const double*,const double*); } extern"C" { void sla__atmt_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*); } extern"C" { void sla_cc62s_ ( const float[6],float*,float*,float*,float*,float*,float*); } extern"C" { void sla_subet_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_cldj_ ( const int*,const int*,const int*,double*,int*); } extern"C" { void sla_dtf2r_ ( const int*,const int*,const double*,double*,int*); } extern"C" { void sla_aop_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,double*,double*); } extern"C" { float sla_rverot_ ( const float*,const float*,const float*,const float*); } extern"C" { void sla_rdplan_ ( const double*,const int*,const double*,const double*,double*,double*,double*); } extern"C" { void sla_dbjin_ ( const char*,int*,double*,int*,int*); } extern"C" { void sla_tpv2c_ ( const float*,const float*,const float[3],float[3],float[3],int*); } extern"C" { void sla_oapqk_ ( const char*,const int*,const double*,const double*,const double[14],double*,double*); } extern"C" { void sla_pda2h_ ( const double*,const double*,const double*,double*,int*,double*,int*); } extern"C" { void sla_aopqk_ ( const double*,const double*,const double[14],double*,double*,double*,double*,double*); } extern"C" { void sla_pm_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { float sla_sepv_ ( const float[3],const float[3]); } extern"C" { void sla_intin_ ( const char*,int*,int*,int*); } extern"C" { void sla_cr2tf_ ( const int*,const float*,char*,int[4]); } extern"C" { void sla_planet_ ( const double*,const int*,double[6],int*); } extern"C" { void sla_vxv_ ( const float[3],const float[3],float[3]); } extern"C" { double sla_zd_ ( const double*,const double*,const double*); } extern"C" { void sla_pcd_ ( const double*,double*,double*); } extern"C" { void sla_ecor_ ( const float*,const float*,const int*,const int*,const float*,float*,float*); } extern"C" { void sla_flotin_ ( const char*,int*,float*,int*); } extern"C" { void sla_dtp2v_ ( const double*,const double*,const double[3],double[3]); } extern"C" { void sla_e2h_ ( const float*,const float*,const float*,float*,float*); } extern"C" { void sla_fk524_ ( const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,double*,double*,double*); } namespace sla { // helper to make a pointer to an rvalue #ifndef POINTABLEDEFD #define POINTABLEDEFD template struct _pointable { T x; _pointable(const T& x): x(x) {} operator T*() { return &x; } }; #endif // - - - - - - // P V 2 E L // - - - - - - // // Heliocentric osculating elements obtained from instantaneous position // and velocity. // // Given: // PV d(6) heliocentric x,y,z,xdot,ydot,zdot of date, // J2000 equatorial triad (AU,AU/s; Note 1) // DATE d date (TT Modified Julian Date = JD-2400000.5) // PMASS d mass of the planet (Sun=1; Note 2) // JFORMR i requested element set (1-3; Note 3) // // Returned: // JFORM d element set actually returned (1-3; Note 4) // EPOCH d epoch of elements (TT MJD) // ORBINC d inclination (radians) // ANODE d longitude of the ascending node (radians) // PERIH d longitude or argument of perihelion (radians) // AORQ d mean distance or perihelion distance (AU) // E d eccentricity // AORL d mean anomaly or longitude (radians, JFORM=1,2 only) // DM d daily motion (radians, JFORM=1 only) // JSTAT i status: 0 = OK // -1 = illegal PMASS // -2 = illegal JFORMR // -3 = position/velocity out of range // // Notes // // 1 The PV 6-vector is with respect to the mean equator and equinox of // epoch J2000. The orbital elements produced are with respect to // the J2000 ecliptic and mean equinox. // // 2 The mass, PMASS, is important only for the larger planets. For // most purposes (e.g. asteroids) use 0D0. Values less than zero // are illegal. // // 3 Three different element-format options are supported: // // Option JFORM=1, suitable for the major planets: // // EPOCH = epoch of elements (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = longitude of perihelion, curly pi (radians) // AORQ = mean distance, a (AU) // E = eccentricity, e // AORL = mean longitude L (radians) // DM = daily motion (radians) // // Option JFORM=2, suitable for minor planets: // // EPOCH = epoch of elements (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = mean distance, a (AU) // E = eccentricity, e // AORL = mean anomaly M (radians) // // Option JFORM=3, suitable for comets: // // EPOCH = epoch of perihelion (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = perihelion distance, q (AU) // E = eccentricity, e // // 4 It may not be possible to generate elements in the form // requested through JFORMR. The caller is notified of the form // of elements actually returned by means of the JFORM argument: // // JFORMR JFORM meaning // // 1 1 OK - elements are in the requested format // 1 2 never happens // 1 3 orbit not elliptical // // 2 1 never happens // 2 2 OK - elements are in the requested format // 2 3 orbit not elliptical // // 3 1 never happens // 3 2 never happens // 3 3 OK - elements are in the requested format // // 5 The arguments returned for each value of JFORM (cf Note 5: JFORM // may not be the same as JFORMR) are as follows: // // JFORM 1 2 3 // EPOCH t0 t0 T // ORBINC i i i // ANODE Omega Omega Omega // PERIH curly pi omega omega // AORQ a a q // E e e e // AORL L M - // DM n - - // // where: // // t0 is the epoch of the elements (MJD, TT) // T " epoch of perihelion (MJD, TT) // i " inclination (radians) // Omega " longitude of the ascending node (radians) // curly pi " longitude of perihelion (radians) // omega " argument of perihelion (radians) // a " mean distance (AU) // q " perihelion distance (AU) // e " eccentricity // L " longitude (radians, 0-2pi) // M " mean anomaly (radians, 0-2pi) // n " daily motion (radians) // - means no value is set // // 6 At very small inclinations, the longitude of the ascending node // ANODE becomes indeterminate and under some circumstances may be // set arbitrarily to zero. Similarly, if the orbit is close to // circular, the true anomaly becomes indeterminate and under some // circumstances may be set arbitrarily to zero. In such cases, // the other elements are automatically adjusted to compensate, // and so the elements remain a valid description of the orbit. // // 7 The osculating epoch for the returned elements is the argument // DATE. // // Reference: Sterne, Theodore E., "An Introduction to Celestial // Mechanics", Interscience Publishers, 1960 // // Called: sla_DRANRM // // Last revision: 8 September 2005 // // Copyright P.T.Wallace. All rights reserved. // struct pv2el { double JFORM;double EPOCH;double ORBINC;double ANODE;double PERIH;double AORQ;double E;double AORL;double DM;int JSTAT; pv2el( const double PV[6],const double DATE,const double PMASS,int JFORMR ) { sla_pv2el_( PV,&DATE,&PMASS,&JFORMR,&JFORM,&EPOCH,&ORBINC,&ANODE,&PERIH,&AORQ,&E,&AORL,&DM,&JSTAT ); } std::string __str__() const { std::ostringstream s; s << "JFORM,EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,JSTAT : " << JFORM<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const aoppat& o ) { return s << o.__str__(); } // - - - - - - // E P B 2 D // - - - - - - // // Conversion of Besselian Epoch to Modified Julian Date // (double precision) // // Given: // EPB dp Besselian Epoch // // The result is the Modified Julian Date (JD - 2400000.5). // // Reference: // Lieske,J.H., 1979. Astron.Astrophys.,73,282. // // P.T.Wallace Starlink February 1984 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline double epb2d( const double EPB ) { return sla_epb2d_( &EPB ); }; // - - - - - - // D A F I N // - - - - - - // // Sexagesimal character string to angle (double precision) // // Given: // STRING c*(*) string containing deg, arcmin, arcsec fields // IPTR i pointer to start of decode (1st = 1) // // Returned: // IPTR i advanced past the decoded angle // A d angle in radians // J i status: 0 = OK // +1 = default, A unchanged // -1 = bad degrees ) // -2 = bad arcminutes ) (note 3) // -3 = bad arcseconds ) // // Example: // // argument before after // // STRING '-57 17 44.806 12 34 56.7' unchanged // IPTR 1 16 (points to 12...) // A ? -1.00000D0 // J ? 0 // // A further call to sla_DAFIN, without adjustment of IPTR, will // decode the second angle, 12deg 34min 56.7sec. // // Notes: // // 1) The first three "fields" in STRING are degrees, arcminutes, // arcseconds, separated by spaces or commas. The degrees field // may be signed, but not the others. The decoding is carried // out by the DFLTIN routine and is free-format. // // 2) Successive fields may be absent, defaulting to zero. For // zero status, the only combinations allowed are degrees alone, // degrees and arcminutes, and all three fields present. If all // three fields are omitted, a status of +1 is returned and A is // unchanged. In all other cases A is changed. // // 3) Range checking: // // The degrees field is not range checked. However, it is // expected to be integral unless the other two fields are absent. // // The arcminutes field is expected to be 0-59, and integral if // the arcseconds field is present. If the arcseconds field // is absent, the arcminutes is expected to be 0-59.9999... // // The arcseconds field is expected to be 0-59.9999... // // 4) Decoding continues even when a check has failed. Under these // circumstances the field takes the supplied value, defaulting // to zero, and the result A is computed and returned. // // 5) Further fields after the three expected ones are not treated // as an error. The pointer IPTR is left in the correct state // for further decoding with the present routine or with DFLTIN // etc. See the example, above. // // 6) If STRING contains hours, minutes, seconds instead of degrees // etc, or if the required units are turns (or days) instead of // radians, the result A should be multiplied as follows: // // for to obtain multiply // STRING A in A by // // d ' " radians 1 = 1D0 // d ' " turns 1/2pi = 0.1591549430918953358D0 // h m s radians 15 = 15D0 // h m s days 15/2pi = 2.3873241463784300365D0 // // Called: sla_DFLTIN // // P.T.Wallace Starlink 1 August 1996 // // Copyright (C) 1996 Rutherford Appleton Laboratory // struct dafin { double A;int J; dafin( const string& STRING,int IPTR ) { sla_dafin_( STRING.c_str(),_pointable(STRING.length()),&IPTR,&A,&J ); } std::string __str__() const { std::ostringstream s; s << "A,J : " << A<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const preces& o ) { return s << o.__str__(); } // - - - - - - // E U L E R // - - - - - - // // Form a rotation matrix from the Euler angles - three successive // rotations about specified Cartesian axes (single precision) // // Given: // ORDER c*(*) specifies about which axes the rotations occur // PHI r 1st rotation (radians) // THETA r 2nd rotation ( " ) // PSI r 3rd rotation ( " ) // // Returned: // RMAT r(3,3) rotation matrix // // A rotation is positive when the reference frame rotates // anticlockwise as seen looking towards the origin from the // positive region of the specified axis. // // The characters of ORDER define which axes the three successive // rotations are about. A typical value is 'ZXZ', indicating that // RMAT is to become the direction cosine matrix corresponding to // rotations of the reference frame through PHI radians about the // old Z-axis, followed by THETA radians about the resulting X-axis, // then PSI radians about the resulting Z-axis. // // The axis names can be any of the following, in any order or // combination: X, Y, Z, uppercase or lowercase, 1, 2, 3. Normal // axis labelling/numbering conventions apply; the xyz (=123) // triad is right-handed. Thus, the 'ZXZ' example given above // could be written 'zxz' or '313' (or even 'ZxZ' or '3xZ'). ORDER // is terminated by length or by the first unrecognized character. // // Fewer than three rotations are acceptable, in which case the later // angle arguments are ignored. If all rotations are zero, the // identity matrix is produced. // // Called: sla_DEULER // // P.T.Wallace Starlink 23 May 1997 // // Copyright (C) 1997 Rutherford Appleton Laboratory // struct euler { float RMAT[3][3]; euler( const string& ORDER,const float PHI,const float THETA,const float PSI ) { sla_euler_( ORDER.c_str(),_pointable(ORDER.length()),&PHI,&THETA,&PSI,RMAT ); } std::string __str__() const { std::ostringstream s; s << "RMAT : " << RMAT << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const euler& o ) { return s << o.__str__(); } // - - - - - // N U T C // - - - - - // // Nutation: longitude & obliquity components and mean obliquity, // using the Shirai & Fukushima (2001) theory. // // Given: // DATE d TDB (loosely ET) as Modified Julian Date // (JD-2400000.5) // Returned: // DPSI,DEPS d nutation in longitude,obliquity // EPS0 d mean obliquity // // Notes: // // 1 The routine predicts forced nutation (but not free core nutation) // plus corrections to the IAU 1976 precession model. // // 2 Earth attitude predictions made by combining the present nutation // model with IAU 1976 precession are accurate to 1 mas (with respect // to the ICRF) for a few decades around 2000. // // 3 The sla_NUTC80 routine is the equivalent of the present routine // but using the IAU 1980 nutation theory. The older theory is less // accurate, leading to errors as large as 350 mas over the interval // 1900-2100, mainly because of the error in the IAU 1976 precession. // // References: // // Shirai, T. & Fukushima, T., Astron.J. 121, 3270-3283 (2001). // // Fukushima, T., Astron.Astrophys. 244, L11 (1991). // // Simon, J. L., Bretagnon, P., Chapront, J., Chapront-Touze, M., // Francou, G. & Laskar, J., Astron.Astrophys. 282, 663 (1994). // // This revision: 24 November 2005 // // Copyright P.T.Wallace. All rights reserved. // struct nutc { double DPSI;double DEPS;double EPS0; nutc( const double DATE ) { sla_nutc_( &DATE,&DPSI,&DEPS,&EPS0 ); } std::string __str__() const { std::ostringstream s; s << "DPSI,DEPS,EPS0 : " << DPSI<<" "< 100 micrometres. // // 3 The routine is a slower but more accurate alternative to the // sla_REFCOQ routine. The constants it produces give perfect // agreement with sla_REFRO at zenith distances arctan(1) (45 deg) // and arctan(4) (about 76 deg). It achieves 0.5 arcsec accuracy // for ZD < 80 deg, 0.01 arcsec accuracy for ZD < 60 deg, and // 0.001 arcsec accuracy for ZD < 45 deg. // // P.T.Wallace Starlink 22 May 2004 // // Copyright (C) 2004 Rutherford Appleton Laboratory // struct refco { double REFA;double REFB; refco( const double HM,const double TDK,const double PMB,const double RH,const double WL,const double PHI,const double TLR,const double EPS ) { sla_refco_( &HM,&TDK,&PMB,&RH,&WL,&PHI,&TLR,&EPS,&REFA,&REFB ); } std::string __str__() const { std::ostringstream s; s << "REFA,REFB : " << REFA<<" "< 100 years) // +1 to +10 = coincident with planet (Note 6) // 0 = OK // -1 = illegal JFORM // -2 = illegal E0 // -3 = illegal AORQ0 // -4 = internal error // -5 = numerical error // // Notes: // // 1 Two different element-format options are available: // // Option JFORM=2, suitable for minor planets: // // EPOCH = epoch of elements (TT MJD) // ORBI = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = mean distance, a (AU) // E = eccentricity, e // AM = mean anomaly M (radians) // // Option JFORM=3, suitable for comets: // // EPOCH = epoch of perihelion (TT MJD) // ORBI = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = perihelion distance, q (AU) // E = eccentricity, e // // 2 DATE0, DATE1, EPOCH0 and EPOCH1 are all instants of time in // the TT timescale (formerly Ephemeris Time, ET), expressed // as Modified Julian Dates (JD-2400000.5). // // DATE0 is the instant at which the given (i.e. unperturbed) // osculating elements are correct. // // DATE1 is the specified instant at which the updated osculating // elements are correct. // // EPOCH0 and EPOCH1 will be the same as DATE0 and DATE1 // (respectively) for the JFORM=2 case, normally used for minor // planets. For the JFORM=3 case, the two epochs will refer to // perihelion passage and so will not, in general, be the same as // DATE0 and/or DATE1 though they may be similar to one another. // // 3 The elements are with respect to the J2000 ecliptic and equinox. // // 4 Unused elements (AM0 and AM1 for JFORM=3) are not accessed. // // 5 See the sla_PERTUE routine for details of the algorithm used. // // 6 This routine is not intended to be used for major planets, which // is why JFORM=1 is not available and why there is no opportunity // to specify either the longitude of perihelion or the daily // motion. However, if JFORM=2 elements are somehow obtained for a // major planet and supplied to the routine, sensible results will, // in fact, be produced. This happens because the sla_PERTUE routine // that is called to perform the calculations checks the separation // between the body and each of the planets and interprets a // suspiciously small value (0.001 AU) as an attempt to apply it to // the planet concerned. If this condition is detected, the // contribution from that planet is ignored, and the status is set to // the planet number (1-10 = Mercury, Venus, EMB, Mars, Jupiter, // Saturn, Uranus, Neptune, Earth, Moon) as a warning. // // Reference: // // Sterne, Theodore E., "An Introduction to Celestial Mechanics", // Interscience Publishers Inc., 1960. Section 6.7, p199. // // Called: sla_EL2UE, sla_PERTUE, sla_UE2EL // // This revision: 19 June 2004 // // Copyright (C) 2004 P.T.Wallace. All rights reserved. // struct pertel { double EPOCH1;double ORBI1;double ANODE1;double PERIH1;double AORQ1;double E1;double AM1;int JSTAT; pertel( const int JFORM,const double DATE0,const double DATE1,const double EPOCH0,const double ORBI0,const double ANODE0,const double PERIH0,const double AORQ0,const double E0,const double AM0 ) { sla_pertel_( &JFORM,&DATE0,&DATE1,&EPOCH0,&ORBI0,&ANODE0,&PERIH0,&AORQ0,&E0,&AM0,&EPOCH1,&ORBI1,&ANODE1,&PERIH1,&AORQ1,&E1,&AM1,&JSTAT ); } std::string __str__() const { std::ostringstream s; s << "EPOCH1,ORBI1,ANODE1,PERIH1,AORQ1,E1,AM1,JSTAT : " << EPOCH1<<" "<(ORDER.length()),&PHI,&THETA,&PSI,RMAT ); } std::string __str__() const { std::ostringstream s; s << "RMAT : " << RMAT << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const deuler& o ) { return s << o.__str__(); } // - - - - - // C S 2 C // - - - - - // // Spherical coordinates to direction cosines (single precision) // // Given: // A,B real spherical coordinates in radians // (RA,Dec), (long,lat) etc. // // Returned: // V real(3) x,y,z unit vector // // The spherical coordinates are longitude (+ve anticlockwise looking // from the +ve latitude pole) and latitude. The Cartesian coordinates // are right handed, with the x axis at zero longitude and latitude, and // the z axis at the +ve latitude pole. // // Last revision: 22 July 2004 // // Copyright P.T.Wallace. All rights reserved. // struct cs2c { float V[3]; cs2c( const float A,const float B ) { sla_cs2c_( &A,&B,V ); } std::string __str__() const { std::ostringstream s; s << "V : " << V << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const cs2c& o ) { return s << o.__str__(); } // - - - - - // I M X V // - - - - - // // Performs the 3-D backward unitary transformation: // // vector VB = (inverse of matrix RM) * vector VA // // (single precision) // // (n.b. the matrix must be unitary, as this routine assumes that // the inverse and transpose are identical) // // Given: // RM real(3,3) matrix // VA real(3) vector // // Returned: // VB real(3) result vector // // P.T.Wallace Starlink November 1984 // // Copyright (C) 1995 Rutherford Appleton Laboratory // struct imxv { float VB[3]; imxv( const float RM[3][3],const float VA[3] ) { sla_imxv_( RM,VA,VB ); } std::string __str__() const { std::ostringstream s; s << "VB : " << VB << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const imxv& o ) { return s << o.__str__(); } // - - - - - - // F K 5 H Z // - - - - - - // // Transform an FK5 (J2000) star position into the frame of the // Hipparcos catalogue, assuming zero Hipparcos proper motion. // // (double precision) // // This routine converts a star position from the FK5 system to // the Hipparcos system, in such a way that the Hipparcos proper // motion is zero. Because such a star has, in general, a non-zero // proper motion in the FK5 system, the routine requires the epoch // at which the position in the FK5 system was determined. // // Given: // R5 d FK5 RA (radians), equinox J2000, epoch EPOCH // D5 d FK5 Dec (radians), equinox J2000, epoch EPOCH // EPOCH d Julian epoch (TDB) // // Returned (all Hipparcos): // RH d RA (radians) // DH d Dec (radians) // // Called: sla_DCS2C, sla_DAV2M, sla_DIMXV, sla_DMXV, sla_DCC2S, // sla_DRANRM // // Notes: // // 1) The FK5 to Hipparcos transformation consists of a pure // rotation and spin; zonal errors in the FK5 catalogue are // not taken into account. // // 2) The published orientation and spin components are interpreted // as "axial vectors". An axial vector points at the pole of the // rotation and its length is the amount of rotation in radians. // // 3) See also sla_FK52H, sla_H2FK5, sla_HFK5Z. // // Reference: // // M.Feissel & F.Mignard, Astron. Astrophys. 331, L33-L36 (1998). // // P.T.Wallace Starlink 22 June 1999 // // Copyright (C) 1999 Rutherford Appleton Laboratory // struct fk5hz { double RH;double DH; fk5hz( const double R5,const double D5,const double EPOCH ) { sla_fk5hz_( &R5,&D5,&EPOCH,&RH,&DH ); } std::string __str__() const { std::ostringstream s; s << "RH,DH : " << RH<<" "< 100 micrometres. // Because the algorithm takes no account of the ionosphere, the // accuracy deteriorates at low frequencies, below about 30 MHz. // // 6 Before use, the value of ZOBS is expressed in the range +/- pi. // If this ranged ZOBS is -ve, the result REF is computed from its // absolute value before being made -ve to match. In addition, if // it has an absolute value greater than 93 deg, a fixed REF value // equal to the result for ZOBS = 93 deg is returned, appropriately // signed. // // 7 As in the original Hohenkerk and Sinclair algorithm, fixed values // of the water vapour polytrope exponent, the height of the // tropopause, and the height at which refraction is negligible are // used. // // 8 The radio refraction has been tested against work done by // Iain Coulson, JACH, (private communication 1995) for the // James Clerk Maxwell Telescope, Mauna Kea. For typical conditions, // agreement at the 0.1 arcsec level is achieved for moderate ZD, // worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg. At hot and // humid sea-level sites the accuracy will not be as good. // // 9 It should be noted that the relative humidity RH is formally // defined in terms of "mixing ratio" rather than pressures or // densities as is often stated. It is the mass of water per unit // mass of dry air divided by that for saturated air at the same // temperature and pressure (see Gill 1982). // // 10 The algorithm is designed for observers in the troposphere. The // supplied temperature, pressure and lapse rate are assumed to be // for a point in the troposphere and are used to define a model // atmosphere with the tropopause at 11km altitude and a constant // temperature above that. However, in practice, the refraction // values returned for stratospheric observers, at altitudes up to // 25km, are quite usable. // // Called: sla_DRANGE, sla__ATMT, sla__ATMS // // Last revision: 5 December 2005 // // Copyright P.T.Wallace. All rights reserved. // struct refro { double REF; refro( const double ZOBS,const double HM,const double TDK,const double PMB,const double RH,const double WL,const double PHI,const double TLR,const double EPS ) { sla_refro_( &ZOBS,&HM,&TDK,&PMB,&RH,&WL,&PHI,&TLR,&EPS,&REF ); } std::string __str__() const { std::ostringstream s; s << "REF : " << REF << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const refro& o ) { return s << o.__str__(); } // - - - - - // A V 2 M // - - - - - // // Form the rotation matrix corresponding to a given axial vector. // // (single precision) // // A rotation matrix describes a rotation about some arbitrary axis, // called the Euler axis. The "axial vector" supplied to this routine // has the same direction as the Euler axis, and its magnitude is the // amount of rotation in radians. // // Given: // AXVEC r(3) axial vector (radians) // // Returned: // RMAT r(3,3) rotation matrix // // If AXVEC is null, the unit matrix is returned. // // The reference frame rotates clockwise as seen looking along // the axial vector from the origin. // // Last revision: 26 November 2005 // // Copyright P.T.Wallace. All rights reserved. // struct av2m { float RMAT[3][3]; av2m( const float AXVEC[3] ) { sla_av2m_( AXVEC,RMAT ); } std::string __str__() const { std::ostringstream s; s << "RMAT : " << RMAT << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const av2m& o ) { return s << o.__str__(); } // - - - - - - // E L 2 U E // - - - - - - // // Transform conventional osculating orbital elements into "universal" // form. // // Given: // DATE d epoch (TT MJD) of osculation (Note 3) // JFORM i choice of element set (1-3, Note 6) // EPOCH d epoch (TT MJD) of the elements // ORBINC d inclination (radians) // ANODE d longitude of the ascending node (radians) // PERIH d longitude or argument of perihelion (radians) // AORQ d mean distance or perihelion distance (AU) // E d eccentricity // AORL d mean anomaly or longitude (radians, JFORM=1,2 only) // DM d daily motion (radians, JFORM=1 only) // // Returned: // U d(13) universal orbital elements (Note 1) // // (1) combined mass (M+m) // (2) total energy of the orbit (alpha) // (3) reference (osculating) epoch (t0) // (4-6) position at reference epoch (r0) // (7-9) velocity at reference epoch (v0) // (10) heliocentric distance at reference epoch // (11) r0.v0 // (12) date (t) // (13) universal eccentric anomaly (psi) of date, approx // // JSTAT i status: 0 = OK // -1 = illegal JFORM // -2 = illegal E // -3 = illegal AORQ // -4 = illegal DM // -5 = numerical error // // Called: sla_UE2PV, sla_PV2UE // // Notes // // 1 The "universal" elements are those which define the orbit for the // purposes of the method of universal variables (see reference). // They consist of the combined mass of the two bodies, an epoch, // and the position and velocity vectors (arbitrary reference frame) // at that epoch. The parameter set used here includes also various // quantities that can, in fact, be derived from the other // information. This approach is taken to avoiding unnecessary // computation and loss of accuracy. The supplementary quantities // are (i) alpha, which is proportional to the total energy of the // orbit, (ii) the heliocentric distance at epoch, (iii) the // outwards component of the velocity at the given epoch, (iv) an // estimate of psi, the "universal eccentric anomaly" at a given // date and (v) that date. // // 2 The companion routine is sla_UE2PV. This takes the set of numbers // that the present routine outputs and uses them to derive the // object's position and velocity. A single prediction requires one // call to the present routine followed by one call to sla_UE2PV; // for convenience, the two calls are packaged as the routine // sla_PLANEL. Multiple predictions may be made by again calling the // present routine once, but then calling sla_UE2PV multiple times, // which is faster than multiple calls to sla_PLANEL. // // 3 DATE is the epoch of osculation. It is in the TT timescale // (formerly Ephemeris Time, ET) and is a Modified Julian Date // (JD-2400000.5). // // 4 The supplied orbital elements are with respect to the J2000 // ecliptic and equinox. The position and velocity parameters // returned in the array U are with respect to the mean equator and // equinox of epoch J2000, and are for the perihelion prior to the // specified epoch. // // 5 The universal elements returned in the array U are in canonical // units (solar masses, AU and canonical days). // // 6 Three different element-format options are available: // // Option JFORM=1, suitable for the major planets: // // EPOCH = epoch of elements (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = longitude of perihelion, curly pi (radians) // AORQ = mean distance, a (AU) // E = eccentricity, e (range 0 to <1) // AORL = mean longitude L (radians) // DM = daily motion (radians) // // Option JFORM=2, suitable for minor planets: // // EPOCH = epoch of elements (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = mean distance, a (AU) // E = eccentricity, e (range 0 to <1) // AORL = mean anomaly M (radians) // // Option JFORM=3, suitable for comets: // // EPOCH = epoch of perihelion (TT MJD) // ORBINC = inclination i (radians) // ANODE = longitude of the ascending node, big omega (radians) // PERIH = argument of perihelion, little omega (radians) // AORQ = perihelion distance, q (AU) // E = eccentricity, e (range 0 to 10) // // 7 Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are // not accessed. // // 8 The algorithm was originally adapted from the EPHSLA program of // D.H.P.Jones (private communication, 1996). The method is based // on Stumpff's Universal Variables. // // Reference: Everhart & Pitkin, Am.J.Phys. 51, 712 (1983). // // Last revision: 8 September 2005 // // Copyright P.T.Wallace. All rights reserved. // struct el2ue { double U[13];int JSTAT; el2ue( const double DATE,const int JFORM,double EPOCH,double ORBINC,double ANODE,double PERIH,double AORQ,double E,double AORL,double DM ) { sla_el2ue_( &DATE,&JFORM,&EPOCH,&ORBINC,&ANODE,&PERIH,&AORQ,&E,&AORL,&DM,U,&JSTAT ); } std::string __str__() const { std::ostringstream s; s << "U,JSTAT : " << U<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const unpcd& o ) { return s << o.__str__(); } // - - - - - - // C D 2 T F // - - - - - - // // Convert an interval in days into hours, minutes, seconds // // (single precision) // // Given: // NDP int number of decimal places of seconds // DAYS real interval in days // // Returned: // SIGN char '+' or '-' // IHMSF int(4) hours, minutes, seconds, fraction // // Notes: // // 1) NDP less than zero is interpreted as zero. // // 2) The largest useful value for NDP is determined by the size of // DAYS, the format of REAL floating-point numbers on the target // machine, and the risk of overflowing IHMSF(4). On some // architectures, for DAYS up to 1.0, the available floating- // point precision corresponds roughly to NDP=3. This is well // below the ultimate limit of NDP=9 set by the capacity of a // typical 32-bit IHMSF(4). // // 3) The absolute value of DAYS may exceed 1.0. In cases where it // does not, it is up to the caller to test for and handle the // case where DAYS is very nearly 1.0 and rounds up to 24 hours, // by testing for IHMSF(1)=24 and setting IHMSF(1-4) to zero. // // Called: sla_DD2TF // // Last revision: 26 December 2004 // // Copyright P.T.Wallace. All rights reserved. // struct cd2tf { char SIGN;int IHMSF[4]; cd2tf( const int NDP,const float DAYS ) { sla_cd2tf_( &NDP,&DAYS,&SIGN,IHMSF ); } std::string __str__() const { std::ostringstream s; s << "SIGN,IHMSF : " << SIGN<<" "< degrees // velocities * (2pi/86400)*(360/2pi) -> degree/sec // accelerations * ((2pi/86400)**2)*(360/2pi) -> degree/sec/sec // // Note that the seconds here are sidereal rather than SI. One // sidereal second is about 0.99727 SI seconds. // // The velocity and acceleration factors assume the sidereal // tracking case. Their respective numerical values are (exactly) // 1/240 and (approximately) 1/3300236.9. // // 2) Azimuth is returned in the range 0-2pi; north is zero, // and east is +pi/2. Elevation and parallactic angle are // returned in the range +/-pi. Parallactic angle is +ve for // a star west of the meridian and is the angle NP-star-zenith. // // 3) The latitude is geodetic as opposed to geocentric. The // hour angle and declination are topocentric. Refraction and // deficiencies in the telescope mounting are ignored. The // purpose of the routine is to give the general form of the // quantities. The details of a real telescope could profoundly // change the results, especially close to the zenith. // // 4) No range checking of arguments is carried out. // // 5) In applications which involve many such calculations, rather // than calling the present routine it will be more efficient to // use inline code, having previously computed fixed terms such // as sine and cosine of latitude, and (for tracking a star) // sine and cosine of declination. // // This revision: 29 October 2004 // // Copyright P.T.Wallace. All rights reserved. // struct altaz { double AZ;double AZD;double AZDD;double EL;double ELD;double ELDD;double PA;double PAD;double PADD; altaz( const double HA,const double DEC,const double PHI ) { sla_altaz_( &HA,&DEC,&PHI,&AZ,&AZD,&AZDD,&EL,&ELD,&ELDD,&PA,&PAD,&PADD ); } std::string __str__() const { std::ostringstream s; s << "AZ,AZD,AZDD,EL,ELD,ELDD,PA,PAD,PADD : " << AZ<<" "<(TYPE.length()),&OB1,&OB2,&DATE,&DUT,&ELONGM,&PHIM,&HM,&XP,&YP,&TDK,&PMB,&RH,&WL,&TLR,&RAP,&DAP ); } std::string __str__() const { std::ostringstream s; s << "RAP,DAP : " << RAP<<" "<(STRING.length()),&IPTR,&A,&J ); } std::string __str__() const { std::ostringstream s; s << "A,J : " << A<<" "< 100 years) // +1 to +10 = coincident with major planet (Note 5) // 0 = OK // -1 = numerical error // // Called: sla_EPJ, sla_PLANET, sla_PV2UE, sla_UE2PV, sla_EPV, // sla_PREC, sla_DMOON, sla_DMXV // // Notes: // // 1 The "universal" elements are those which define the orbit for the // purposes of the method of universal variables (see reference 2). // They consist of the combined mass of the two bodies, an epoch, // and the position and velocity vectors (arbitrary reference frame) // at that epoch. The parameter set used here includes also various // quantities that can, in fact, be derived from the other // information. This approach is taken to avoiding unnecessary // computation and loss of accuracy. The supplementary quantities // are (i) alpha, which is proportional to the total energy of the // orbit, (ii) the heliocentric distance at epoch, (iii) the // outwards component of the velocity at the given epoch, (iv) an // estimate of psi, the "universal eccentric anomaly" at a given // date and (v) that date. // // 2 The universal elements are with respect to the J2000 equator and // equinox. // // 3 The epochs DATE, U(3) and U(12) are all Modified Julian Dates // (JD-2400000.5). // // 4 The algorithm is a simplified form of Encke's method. It takes as // a basis the unperturbed motion of the body, and numerically // integrates the perturbing accelerations from the major planets. // The expression used is essentially Sterne's 6.7-2 (reference 1). // Everhart and Pitkin (reference 2) suggest rectifying the orbit at // each integration step by propagating the new perturbed position // and velocity as the new universal variables. In the present // routine the orbit is rectified less frequently than this, in order // to gain a slight speed advantage. However, the rectification is // done directly in terms of position and velocity, as suggested by // Everhart and Pitkin, bypassing the use of conventional orbital // elements. // // The f(q) part of the full Encke method is not used. The purpose // of this part is to avoid subtracting two nearly equal quantities // when calculating the "indirect member", which takes account of the // small change in the Sun's attraction due to the slightly displaced // position of the perturbed body. A simpler, direct calculation in // double precision proves to be faster and not significantly less // accurate. // // Apart from employing a variable timestep, and occasionally // "rectifying the orbit" to keep the indirect member small, the // integration is done in a fairly straightforward way. The // acceleration estimated for the middle of the timestep is assumed // to apply throughout that timestep; it is also used in the // extrapolation of the perturbations to the middle of the next // timestep, to predict the new disturbed position. There is no // iteration within a timestep. // // Measures are taken to reach a compromise between execution time // and accuracy. The starting-point is the goal of achieving // arcsecond accuracy for ordinary minor planets over a ten-year // timespan. This goal dictates how large the timesteps can be, // which in turn dictates how frequently the unperturbed motion has // to be recalculated from the osculating elements. // // Within predetermined limits, the timestep for the numerical // integration is varied in length in inverse proportion to the // magnitude of the net acceleration on the body from the major // planets. // // The numerical integration requires estimates of the major-planet // motions. Approximate positions for the major planets (Pluto // alone is omitted) are obtained from the routine sla_PLANET. Two // levels of interpolation are used, to enhance speed without // significantly degrading accuracy. At a low frequency, the routine // sla_PLANET is called to generate updated position+velocity "state // vectors". The only task remaining to be carried out at the full // frequency (i.e. at each integration step) is to use the state // vectors to extrapolate the planetary positions. In place of a // strictly linear extrapolation, some allowance is made for the // curvature of the orbit by scaling back the radius vector as the // linear extrapolation goes off at a tangent. // // Various other approximations are made. For example, perturbations // by Pluto and the minor planets are neglected and relativistic // effects are not taken into account. // // In the interests of simplicity, the background calculations for // the major planets are carried out en masse. The mean elements and // state vectors for all the planets are refreshed at the same time, // without regard for orbit curvature, mass or proximity. // // The Earth-Moon system is treated as a single body when the body is // distant but as separate bodies when closer to the EMB than the // parameter RNE, which incurs a time penalty but improves accuracy // for near-Earth objects. // // 5 This routine is not intended to be used for major planets. // However, if major-planet elements are supplied, sensible results // will, in fact, be produced. This happens because the routine // checks the separation between the body and each of the planets and // interprets a suspiciously small value (0.001 AU) as an attempt to // apply the routine to the planet concerned. If this condition is // detected, the contribution from that planet is ignored, and the // status is set to the planet number (1-10 = Mercury, Venus, EMB, // Mars, Jupiter, Saturn, Uranus, Neptune, Earth, Moon) as a warning. // // References: // // 1 Sterne, Theodore E., "An Introduction to Celestial Mechanics", // Interscience Publishers Inc., 1960. Section 6.7, p199. // // 2 Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983. // // Last revision: 27 December 2004 // // Copyright P.T.Wallace. All rights reserved. // struct pertue { int JSTAT; pertue( const double DATE,const double U[13] ) { sla_pertue_( &DATE,U,&JSTAT ); } std::string __str__() const { std::ostringstream s; s << "JSTAT : " << JSTAT << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const pertue& o ) { return s << o.__str__(); } // - - - - // E P J // - - - - // // Conversion of Modified Julian Date to Julian Epoch (double precision) // // Given: // DATE dp Modified Julian Date (JD - 2400000.5) // // The result is the Julian Epoch. // // Reference: // Lieske,J.H., 1979. Astron.Astrophys.,73,282. // // P.T.Wallace Starlink February 1984 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline double epj( const double DATE ) { return sla_epj_( &DATE ); }; // - - - - - - - // A T M D S P // - - - - - - - // // Apply atmospheric-dispersion adjustments to refraction coefficients. // // Given: // TDK d ambient temperature, K // PMB d ambient pressure, millibars // RH d ambient relative humidity, 0-1 // WL1 d reference wavelength, micrometre (0.4D0 recommended) // A1 d refraction coefficient A for wavelength WL1 (radians) // B1 d refraction coefficient B for wavelength WL1 (radians) // WL2 d wavelength for which adjusted A,B required // // Returned: // A2 d refraction coefficient A for wavelength WL2 (radians) // B2 d refraction coefficient B for wavelength WL2 (radians) // // Notes: // // 1 To use this routine, first call sla_REFCO specifying WL1 as the // wavelength. This yields refraction coefficients A1,B1, correct // for that wavelength. Subsequently, calls to sla_ATMDSP specifying // different wavelengths will produce new, slightly adjusted // refraction coefficients which apply to the specified wavelength. // // 2 Most of the atmospheric dispersion happens between 0.7 micrometre // and the UV atmospheric cutoff, and the effect increases strongly // towards the UV end. For this reason a blue reference wavelength // is recommended, for example 0.4 micrometres. // // 3 The accuracy, for this set of conditions: // // height above sea level 2000 m // latitude 29 deg // pressure 793 mb // temperature 17 degC // humidity 50% // lapse rate 0.0065 degC/m // reference wavelength 0.4 micrometre // star elevation 15 deg // // is about 2.5 mas RMS between 0.3 and 1.0 micrometres, and stays // within 4 mas for the whole range longward of 0.3 micrometres // (compared with a total dispersion from 0.3 to 20.0 micrometres // of about 11 arcsec). These errors are typical for ordinary // conditions and the given elevation; in extreme conditions values // a few times this size may occur, while at higher elevations the // errors become much smaller. // // 4 If either wavelength exceeds 100 micrometres, the radio case // is assumed and the returned refraction coefficients are the // same as the given ones. Note that radio refraction coefficients // cannot be turned into optical values using this routine, nor // vice versa. // // 5 The algorithm consists of calculation of the refractivity of the // air at the observer for the two wavelengths, using the methods // of the sla_REFRO routine, and then scaling of the two refraction // coefficients according to classical refraction theory. This // amounts to scaling the A coefficient in proportion to (n-1) and // the B coefficient almost in the same ratio (see R.M.Green, // "Spherical Astronomy", Cambridge University Press, 1985). // // Last revision 2 December 2005 // // Copyright P.T.Wallace. All rights reserved. // struct atmdsp { double A2;double B2; atmdsp( const double TDK,const double PMB,const double RH,const double WL1,const double A1,const double B1,const double WL2 ) { sla_atmdsp_( &TDK,&PMB,&RH,&WL1,&A1,&B1,&WL2,&A2,&B2 ); } std::string __str__() const { std::ostringstream s; s << "A2,B2 : " << A2<<" "< 100 micrometres. // // Notes: // // 1 The model is an approximation, for moderate zenith distances, // to the predictions of the sla_REFRO routine. The approximation // is maintained across a range of conditions, and applies to // both optical/IR and radio. // // 2 The algorithm is a fast alternative to the sla_REFCO routine. // The latter calls the sla_REFRO routine itself: this involves // integrations through a model atmosphere, and is costly in // processor time. However, the model which is produced is precisely // correct for two zenith distance (45 degrees and about 76 degrees) // and at other zenith distances is limited in accuracy only by the // A tan Z + B tan**3 Z formulation itself. The present routine // is not as accurate, though it satisfies most practical // requirements. // // 3 The model omits the effects of (i) height above sea level (apart // from the reduced pressure itself), (ii) latitude (i.e. the // flattening of the Earth) and (iii) variations in tropospheric // lapse rate. // // The model was tested using the following range of conditions: // // lapse rates 0.0055, 0.0065, 0.0075 K/metre // latitudes 0, 25, 50, 75 degrees // heights 0, 2500, 5000 metres ASL // pressures mean for height -10% to +5% in steps of 5% // temperatures -10 deg to +20 deg with respect to 280 deg at SL // relative humidity 0, 0.5, 1 // wavelengths 0.4, 0.6, ... 2 micron, + radio // zenith distances 15, 45, 75 degrees // // The accuracy with respect to direct use of the sla_REFRO routine // was as follows: // // worst RMS // // optical/IR 62 mas 8 mas // radio 319 mas 49 mas // // For this particular set of conditions: // // lapse rate 0.0065 K/metre // latitude 50 degrees // sea level // pressure 1005 mb // temperature 280.15 K // humidity 80% // wavelength 5740 Angstroms // // the results were as follows: // // ZD sla_REFRO sla_REFCOQ Saastamoinen // // 10 10.27 10.27 10.27 // 20 21.19 21.20 21.19 // 30 33.61 33.61 33.60 // 40 48.82 48.83 48.81 // 45 58.16 58.18 58.16 // 50 69.28 69.30 69.27 // 55 82.97 82.99 82.95 // 60 100.51 100.54 100.50 // 65 124.23 124.26 124.20 // 70 158.63 158.68 158.61 // 72 177.32 177.37 177.31 // 74 200.35 200.38 200.32 // 76 229.45 229.43 229.42 // 78 267.44 267.29 267.41 // 80 319.13 318.55 319.10 // // deg arcsec arcsec arcsec // // The values for Saastamoinen's formula (which includes terms // up to tan^5) are taken from Hohenkerk and Sinclair (1985). // // The results from the much slower but more accurate sla_REFCO // routine have not been included in the tabulation as they are // identical to those in the sla_REFRO column to the 0.01 arcsec // resolution used. // // 4 Outlandish input parameters are silently limited to mathematically // safe values. Zero pressure is permissible, and causes zeroes to // be returned. // // 5 The algorithm draws on several sources, as follows: // // a) The formula for the saturation vapour pressure of water as // a function of temperature and temperature is taken from // expressions A4.5-A4.7 of Gill (1982). // // b) The formula for the water vapour pressure, given the // saturation pressure and the relative humidity, is from // Crane (1976), expression 2.5.5. // // c) The refractivity of air is a function of temperature, // total pressure, water-vapour pressure and, in the case // of optical/IR but not radio, wavelength. The formulae // for the two cases are developed from Hohenkerk & Sinclair // (1985) and Rueger (2002). // // The above three items are as used in the sla_REFRO routine. // // d) The formula for beta, the ratio of the scale height of the // atmosphere to the geocentric distance of the observer, is // an adaption of expression 9 from Stone (1996). The // adaptations, arrived at empirically, consist of (i) a // small adjustment to the coefficient and (ii) a humidity // term for the radio case only. // // e) The formulae for the refraction constants as a function of // n-1 and beta are from Green (1987), expression 4.31. // // References: // // Crane, R.K., Meeks, M.L. (ed), "Refraction Effects in the Neutral // Atmosphere", Methods of Experimental Physics: Astrophysics 12B, // Academic Press, 1976. // // Gill, Adrian E., "Atmosphere-Ocean Dynamics", Academic Press, 1982. // // Green, R.M., "Spherical Astronomy", Cambridge University Press, 1987. // // Hohenkerk, C.Y., & Sinclair, A.T., NAO Technical Note No. 63, 1985. // // Rueger, J.M., "Refractive Index Formulae for Electronic Distance // Measurement with Radio and Millimetre Waves", in Unisurv Report // S-68, School of Surveying and Spatial Information Systems, // University of New South Wales, Sydney, Australia, 2002. // // Stone, Ronald C., P.A.S.P. 108 1051-1058, 1996. // // Last revision: 2 December 2005 // // Copyright P.T.Wallace. All rights reserved. // struct refcoq { double REFA;double REFB; refcoq( const double TDK,const double PMB,const double RH,const double WL ) { sla_refcoq_( &TDK,&PMB,&RH,&WL,&REFA,&REFB ); } std::string __str__() const { std::ostringstream s; s << "REFA,REFB : " << REFA<<" "<(TYPE.length()),&OB1,&OB2,AOPRMS,&RAP,&DAP ); } std::string __str__() const { std::ostringstream s; s << "RAP,DAP : " << RAP<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const pcd& o ) { return s << o.__str__(); } // - - - - - // E C O R // - - - - - // // Component of Earth orbit velocity and heliocentric // light time in a given direction (single precision) // // Given: // RM,DM real mean RA, Dec of date (radians) // IY int year // ID int day in year (1 = Jan 1st) // FD real fraction of day // // Returned: // RV real component of Earth orbital velocity (km/sec) // TL real component of heliocentric light time (sec) // // Notes: // // 1 The date and time is TDB (loosely ET) in a Julian calendar // which has been aligned to the ordinary Gregorian // calendar for the interval 1900 March 1 to 2100 February 28. // The year and day can be obtained by calling sla_CALYD or // sla_CLYD. // // 2 Sign convention: // // The velocity component is +ve when the Earth is receding from // the given point on the sky. The light time component is +ve // when the Earth lies between the Sun and the given point on // the sky. // // 3 Accuracy: // // The velocity component is usually within 0.004 km/s of the // correct value and is never in error by more than 0.007 km/s. // The error in light time correction is about 0.03s at worst, // but is usually better than 0.01s. For applications requiring // higher accuracy, see the sla_EVP and sla_EPV routines. // // Called: sla_EARTH, sla_CS2C, sla_VDV // // Last revision: 5 April 2005 // // Copyright P.T.Wallace. All rights reserved. // struct ecor { float RV;float TL; ecor( const float RM,const float DM,const int IY,const int ID,const float FD ) { sla_ecor_( &RM,&DM,&IY,&ID,&FD,&RV,&TL ); } std::string __str__() const { std::ostringstream s; s << "RV,TL : " << RV<<" "<