#ifndef SLAAINCLD #define SLAAINCLD #include #include #include using std::string; using std::ostream; using std::endl; extern"C" { void sla_dmoon_ ( const double*,double[6]); } extern"C" { double sla_dt_ ( const double*); } extern"C" { void sla_fk524_ ( const double*,const double*,const double*,const double*,const double*,const double*,double*,double*,double*,double*,double*,double*); } extern"C" { void sla_kbj_ ( const int*,const double*,char*,int*); } extern"C" { void sla_epv_ ( const double*,double[3],double[3],double[3],double[3]); } extern"C" { void sla_ds2c6_ ( const double*,const double*,const double*,const double*,const double*,const double*,double[6]); } extern"C" { void sla_dtpv2c_ ( const double*,const double*,const double[3],double[3],double[3],int*); } extern"C" { void sla_rdplan_ ( const double*,const int*,const double*,const double*,double*,double*,double*); } extern"C" { void sla_planet_ ( const double*,const int*,double[6],int*); } extern"C" { void sla_dr2tf_ ( const int*,const double*,char*,int[4]); } extern"C" { void sla_cldj_ ( const int*,const int*,const int*,double*,int*); } extern"C" { void sla_dm2av_ ( const double[3][3],double[3]); } 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_aopqk_ ( const double*,const double*,const double[14],double*,double*,double*,double*,double*); } extern"C" { void sla_el2ue_ ( const double*,const int*,double*,double*,double*,double*,double*,double*,double*,double*,double[13],int*); } extern"C" { void sla_dtp2v_ ( const double*,const double*,const double[3],double[3]); } extern"C" { float sla_bear_ ( const float*,const float*,const float*,const float*); } extern"C" { void sla_daf2r_ ( const int*,const int*,const double*,double*,int*); } extern"C" { void sla_dtf2d_ ( const int*,const int*,const double*,double*,int*); } 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" { void sla_flotin_ ( const char*,int*,float*,int*); } 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" { void sla_dcc2s_ ( const double[3],double*,double*); } extern"C" { void sla_dtf2r_ ( const int*,const int*,const double*,double*,int*); } extern"C" { double sla_pa_ ( const double*,const double*,const double*); } extern"C" { double sla_zd_ ( const double*,const double*,const double*); } extern"C" { void sla_m2av_ ( const float[3][3],float[3]); } extern"C" { void sla_dimxv_ ( const double[3][3],const double[3],double[3]); } extern"C" { void sla_tpv2c_ ( const float*,const float*,const float[3],float[3],float[3],int*); } extern"C" { void sla_dd2tf_ ( 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" { void sla_deuler_ ( const char*,const int*,const double*,const double*,const double*,double[3][3]); } 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_galeq_ ( const double*,const double*,double*,double*); } extern"C" { void sla_oapqk_ ( const char*,const int*,const double*,const double*,const double[14],double*,double*); } 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_map_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_etrms_ ( const double*,double[3]); } extern"C" { void sla_cs2c_ ( const float*,const float*,float[3]); } extern"C" { void sla_ctf2d_ ( const int*,const int*,const float*,float*,int*); } extern"C" { void sla_dvxv_ ( const double[3],const double[3],double[3]); } extern"C" { void sla_plantu_ ( const double*,const double*,const double*,const double[13],double*,double*,double*,int*); } extern"C" { void sla_dr2af_ ( const int*,const double*,char*,int[4]); } extern"C" { void sla_mxv_ ( const float[3][3],const float[3],float[3]); } extern"C" { void sla_aoppat_ ( const double*,double[14]); } extern"C" { void sla_caf2r_ ( const int*,const int*,const float*,float*,int*); } extern"C" { double sla_dranrm_ ( const double*); } extern"C" { void sla_de2h_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_h2fk5_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { void sla_ampqk_ ( const double*,const double*,const double[21],double*,double*); } extern"C" { void sla_dmxm_ ( const double[3][3],const double[3][3],double[3][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_dtp2s_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_pv2ue_ ( const double[6],const double*,const double*,double[13],int*); } extern"C" { double sla_dbear_ ( const double*,const double*,const double*,const double*); } extern"C" { void sla_intin_ ( const char*,int*,int*,int*); } extern"C" { void sla_dvn_ ( const double[3],double[3],double*); } extern"C" { void sla_pertue_ ( const double*,const double[13],int*); } extern"C" { void sla_s2tp_ ( const float*,const float*,const float*,const float*,float*,float*,int*); } extern"C" { void sla_fk54z_ ( const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { float sla_ranorm_ ( const double*); } extern"C" { void sla_afin_ ( const char*,const int*,int*,float*,int*); } extern"C" { void sla_dtps2c_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*,int*); } extern"C" { double sla_dsepv_ ( const double[3],const double[3]); } extern"C" { double sla_gmst_ ( const double*); } extern"C" { double sla_epco_ ( const char*,const char*,const double*); } extern"C" { double sla_epb2d_ ( const double*); } extern"C" { float sla_pav_ ( const float[3],const float[3]); } extern"C" { void sla_earth_ ( const int*,const int*,const float*,float[6]); } extern"C" { void sla_fk5hz_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_moon_ ( const int*,const int*,const float*,float[6]); } extern"C" { void sla_cr2af_ ( const int*,const float*,char*,int[4]); } extern"C" { void sla_geoc_ ( const double*,const double*,double*,double*); } extern"C" { void sla_djcal_ ( int*,const double*,int[4],int*); } extern"C" { void sla_cd2tf_ ( const int*,const float*,char*,int[4]); } extern"C" { double sla_eqeqx_ ( const double*); } extern"C" { double sla_epj2d_ ( const double*); } extern"C" { void sla_tp2v_ ( const float*,const float*,const float[3],float[3]); } extern"C" { void sla_cc62s_ ( const float[6],float*,float*,float*,float*,float*,float*); } extern"C" { void sla_tps2c_ ( const float*,const float*,const float*,const float*,float*,float*,float*,float*,int*); } extern"C" { float sla_rvlsrd_ ( const float*,const float*); } extern"C" { float sla_sepv_ ( const float[3],const float[3]); } extern"C" { void sla_ecleq_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla__idchi_ ( const char*,int*,int*,double*); } extern"C" { void sla_dafin_ ( const char*,const int*,int*,double*,int*); } extern"C" { void sla_ue2el_ ( const double[13],int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*); } extern"C" { void sla_imxv_ ( const float[3][3],const float[3],float[3]); } 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_dh2e_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_pda2h_ ( const double*,const double*,const double*,double*,int*,double*,int*); } 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_dv2tp_ ( const double[3],const double[3],double*,double*,int*); } extern"C" { void sla_nutc_ ( const double*,double*,double*,double*); } extern"C" { void sla_refcoq_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { void sla_dc62s_ ( const double[6],double*,double*,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" { void sla_dcs2c_ ( const double*,const double*,double[3]); } extern"C" { void sla_pcd_ ( const double*,double*,double*); } extern"C" { double sla_dtt_ ( const double*); } extern"C" { double sla_epb_ ( const double*); } extern"C" { void sla_dav2m_ ( const double[3],double[3][3]); } extern"C" { void sla_ge50_ ( const double*,const double*,double*,double*); } extern"C" { void sla_cc2s_ ( const float[3],float*,float*); } extern"C" { void sla_addet_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_precl_ ( const double*,const double*,double[3][3]); } extern"C" { double sla_gmsta_ ( const double*,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_cr2tf_ ( const int*,const float*,char*,int[4]); } extern"C" { void sla_ds2tp_ ( const double*,const double*,const double*,const double*,double*,double*,int*); } extern"C" { void sla_ecor_ ( const float*,const float*,const int*,const int*,const float*,float*,float*); } extern"C" { void sla_ue2pv_ ( const double*,const double[13],double[6],int*); } extern"C" { void sla_xy2xy_ ( const double*,const double*,const double[6],double*,double*); } extern"C" { float sla_sep_ ( const float*,const float*,const float*,const float*); } extern"C" { float sla_range_ ( const double*); } extern"C" { double sla_dsep_ ( const double*,const double*,const double*,const double*); } extern"C" { void sla_tp2s_ ( const float*,const float*,const float*,const float*,float*,float*); } extern"C" { void sla_pdq2h_ ( const double*,const double*,const double*,double*,int*,double*,int*); } extern"C" { double sla_epj_ ( const double*); } extern"C" { void sla_mapqkz_ ( const double*,const double*,const double[21],double*,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_euler_ ( const char*,const int*,const float*,const float*,const float*,float[3][3]); } extern"C" { float sla_rvlg_ ( const float*,const float*); } extern"C" { void sla__idchf_ ( const char*,int*,int*,int*,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_prenut_ ( const double*,const double*,double[3][3]); } extern"C" { void sla_preces_ ( const char*,const double*,const double*,double*,double*); } extern"C" { void sla_djcl_ ( const double*,int*,int*,int*,double*,int*); } extern"C" { void sla_invf_ ( const double[6],double[6],int*); } extern"C" { void sla_dbjin_ ( const char*,int*,double*,int*,int*); } extern"C" { void sla_mappa_ ( const double*,const double*,double[21]); } extern"C" { void sla_nutc80_ ( const double*,double*,double*,double*); } extern"C" { void sla_eg50_ ( const double*,const double*,double*,double*); } extern"C" { void sla_ctf2r_ ( const int*,const int*,const float*,float*,int*); } extern"C" { void sla_vn_ ( const float[3],float[3],float*); } extern"C" { double sla_drange_ ( const double*); } extern"C" { float sla_rvgalc_ ( const float*,const float*); } extern"C" { void sla_polmo_ ( const double*,const double*,double*,const double*,double*,double*,double*); } extern"C" { void sla_caldj_ ( const int*,const int*,const int*,double*,int*); } extern"C" { void sla_e2h_ ( const float*,const float*,const float*,float*,float*); } extern"C" { void sla_subet_ ( const double*,const double*,const double*,double*,double*); } extern"C" { double sla_dvdv_ ( const double[3],const double[3]); } extern"C" { void sla_calyd_ ( const int*,const int*,const int*,int*,int*,int*); } extern"C" { void sla_galsup_ ( const double*,const double*,double*,double*); } extern"C" { void sla_ecmat_ ( const double*,double[3][3]); } extern"C" { void sla_mapqk_ ( const double*,const double*,const double*,const double*,const double*,const double*,const double[21],double*,double*); } extern"C" { void sla_vxv_ ( const float[3],const float[3],float[3]); } extern"C" { void sla_dmxv_ ( const double[3][3],const double[3],double[3]); } 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_refv_ ( const double*,const double*,const double*,double*); } extern"C" { void sla_supgal_ ( const double*,const double*,double*,double*); } extern"C" { void sla_h2e_ ( const float*,const float*,const float*,float*,float*); } extern"C" { void sla_amp_ ( const double*,const double*,const double*,const double*,double*,double*); } extern"C" { float sla_rvlsrk_ ( const float*,const float*); } extern"C" { void sla_dfltin_ ( const char*,int*,double*,int*); } extern"C" { void sla_unpcd_ ( const double*,double*,double*); } extern"C" { void sla_prebn_ ( const double*,const double*,double[3][3]); } extern"C" { double sla_dat_ ( const double*); } extern"C" { void sla_fk45z_ ( const double*,const double*,const double*,double*,double*); } extern"C" { float sla_vdv_ ( const float[3],const float[3]); } extern"C" { void sla_dcmpf_ ( const double[6],double*,double*,double*,double*,double*,double*); } extern"C" { void sla_av2m_ ( const float[3],float[3][3]); } extern"C" { void sla_clyd_ ( const int*,const int*,const int*,int*,int*,int*); } extern"C" { void sla_fk52h_ ( const double*,const double*,const double*,const double*,double*,double*,double*,double*); } extern"C" { double sla_rcc_ ( const double*,const double*,const double*,const double*,const double*); } extern"C" { void sla_cs2c6_ ( const float*,const float*,const float*,const float*,const float*,const float*,float[6]); } extern"C" { void sla_prec_ ( const double*,const double*,double[3][3]); } extern"C" { void sla_eqecl_ ( const double*,const double*,const double*,double*,double*); } extern"C" { void sla_eqgal_ ( const double*,const double*,double*,double*); } extern"C" { void sla_mxm_ ( const float[3][3],const float[3][3],float[3][3]); } extern"C" { double sla_dpav_ ( const double[3],const double[3]); } extern"C" { void sla_refz_ ( const double*,const double*,const 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" { double sla_airmas_ ( 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*); } 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 // - - - - - - // D M O O N // - - - - - - // // Approximate geocentric position and velocity of the Moon // (double precision) // // Given: // DATE D TDB (loosely ET) as a Modified Julian Date // (JD-2400000.5) // // Returned: // PV D(6) Moon x,y,z,xdot,ydot,zdot, mean equator and // equinox of date (AU, AU/s) // // Notes: // // 1 This routine is a full implementation of the algorithm // published by Meeus (see reference). // // 2 Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in // latitude and 0.2 arcsec in HP (equivalent to about 20 km in // distance). Comparison with JPL DE200 over the interval // 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in // longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km // and 81 mm/s in distance. The maximum errors over the same // interval are 18 arcsec and 0.50 arcsec/hour in longitude, // 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s // in distance. // // 3 The original algorithm is expressed in terms of the obsolete // timescale Ephemeris Time. Either TDB or TT can be used, but // not UT without incurring significant errors (30 arcsec at // the present time) due to the Moon's 0.5 arcsec/sec movement. // // 4 The algorithm is based on pre IAU 1976 standards. However, // the result has been moved onto the new (FK5) equinox, an // adjustment which is in any case much smaller than the // intrinsic accuracy of the procedure. // // 5 Velocity is obtained by a complete analytical differentiation // of the Meeus model. // // Reference: // Meeus, l'Astronomie, June 1984, p348. // // P.T.Wallace Starlink 22 January 1998 // // Copyright (C) 1998 Rutherford Appleton Laboratory // struct dmoon { double PV[6]; dmoon( const double DATE ) { sla_dmoon_( &DATE,PV ); } std::string __str__() const { std::ostringstream s; s << "PV : " << PV << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const dmoon& o ) { return s << o.__str__(); } // - - - // D T // - - - // // Estimate the offset between dynamical time and Universal Time // for a given historical epoch. // // Given: // EPOCH d (Julian) epoch (e.g. 1850D0) // // The result is a rough estimate of ET-UT (after 1984, TT-UT) at // the given epoch, in seconds. // // Notes: // // 1 Depending on the epoch, one of three parabolic approximations // is used: // // before 979 Stephenson & Morrison's 390 BC to AD 948 model // 979 to 1708 Stephenson & Morrison's 948 to 1600 model // after 1708 McCarthy & Babcock's post-1650 model // // The breakpoints are chosen to ensure continuity: they occur // at places where the adjacent models give the same answer as // each other. // // 2 The accuracy is modest, with errors of up to 20 sec during // the interval since 1650, rising to perhaps 30 min by 1000 BC. // Comparatively accurate values from AD 1600 are tabulated in // the Astronomical Almanac (see section K8 of the 1995 AA). // // 3 The use of double-precision for both argument and result is // purely for compatibility with other SLALIB time routines. // // 4 The models used are based on a lunar tidal acceleration value // of -26.00 arcsec per century. // // Reference: Explanatory Supplement to the Astronomical Almanac, // ed P.K.Seidelmann, University Science Books (1992), // section 2.553, p83. This contains references to // the Stephenson & Morrison and McCarthy & Babcock // papers. // // P.T.Wallace Starlink 1 March 1995 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline double dt( const double EPOCH ) { return sla_dt_( &EPOCH ); }; // - - - - - - // F K 5 2 4 // - - - - - - // // Convert J2000.0 FK5 star data to B1950.0 FK4 (double precision) // // This routine converts stars from the new, IAU 1976, FK5, Fricke // system, to the old, Bessel-Newcomb, FK4 system. The precepts // of Smith et al (Ref 1) are followed, using the implementation // by Yallop et al (Ref 2) of a matrix method due to Standish. // Kinoshita's development of Andoyer's post-Newcomb precession is // used. The numerical constants from Seidelmann et al (Ref 3) are // used canonically. // // Given: (all J2000.0,FK5) // R2000,D2000 dp J2000.0 RA,Dec (rad) // DR2000,DD2000 dp J2000.0 proper motions (rad/Jul.yr) // P2000 dp parallax (arcsec) // V2000 dp radial velocity (km/s, +ve = moving away) // // Returned: (all B1950.0,FK4) // R1950,D1950 dp B1950.0 RA,Dec (rad) // DR1950,DD1950 dp B1950.0 proper motions (rad/trop.yr) // P1950 dp parallax (arcsec) // V1950 dp radial velocity (km/s, +ve = moving away) // // Notes: // // 1) The proper motions in RA are dRA/dt rather than // cos(Dec)*dRA/dt, and are per year rather than per century. // // 2) Note that conversion from Julian epoch 2000.0 to Besselian // epoch 1950.0 only is provided for. Conversions involving // other epochs will require use of the appropriate precession, // proper motion, and E-terms routines before and/or after // FK524 is called. // // 3) In the FK4 catalogue the proper motions of stars within // 10 degrees of the poles do not embody the differential // E-term effect and should, strictly speaking, be handled // in a different manner from stars outside these regions. // However, given the general lack of homogeneity of the star // data available for routine astrometry, the difficulties of // handling positions that may have been determined from // astrometric fields spanning the polar and non-polar regions, // the likelihood that the differential E-terms effect was not // taken into account when allowing for proper motion in past // astrometry, and the undesirability of a discontinuity in // the algorithm, the decision has been made in this routine to // include the effect of differential E-terms on the proper // motions for all stars, whether polar or not. At epoch 2000, // and measuring on the sky rather than in terms of dRA, the // errors resulting from this simplification are less than // 1 milliarcsecond in position and 1 milliarcsecond per // century in proper motion. // // References: // // 1 Smith, C.A. et al, 1989. "The transformation of astrometric // catalog systems to the equinox J2000.0". Astron.J. 97, 265. // // 2 Yallop, B.D. et al, 1989. "Transformation of mean star places // from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". // Astron.J. 97, 274. // // 3 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to // the Astronomical Almanac", ISBN 0-935702-68-7. // // P.T.Wallace Starlink 19 December 1993 // // Copyright (C) 1995 Rutherford Appleton Laboratory // struct fk524 { double R1950;double D1950;double DR1950;double DD1950;double P1950;double V1950; fk524( const double R2000,const double D2000,const double DR2000,const double DD2000,const double P2000,const double V2000 ) { sla_fk524_( &R2000,&D2000,&DR2000,&DD2000,&P2000,&V2000,&R1950,&D1950,&DR1950,&DD1950,&P1950,&V1950 ); } std::string __str__() const { std::ostringstream s; s << "R1950,D1950,DR1950,DD1950,P1950,V1950 : " << R1950<<" "<(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__(); } // - - - - - - - // P L A N E L // - - - - - - - // // Heliocentric position and velocity of a planet, asteroid or comet, // starting from orbital elements. // // Given: // DATE d date, Modified Julian Date (JD - 2400000.5, Note 1) // JFORM i choice of element set (1-3; Note 3) // EPOCH d epoch of elements (TT MJD, Note 4) // 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: // PV d(6) heliocentric x,y,z,xdot,ydot,zdot of date, // J2000 equatorial triad (AU,AU/s) // JSTAT i status: 0 = OK // -1 = illegal JFORM // -2 = illegal E // -3 = illegal AORQ // -4 = illegal DM // -5 = numerical error // // Called: sla_EL2UE, sla_UE2PV // // Notes // // 1 DATE is the instant for which the prediction is required. It is // in the TT timescale (formerly Ephemeris Time, ET) and is a // Modified Julian Date (JD-2400000.5). // // 2 The elements are with respect to the J2000 ecliptic and equinox. // // 3 A choice of three different element-set options is 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 elements and 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) // // Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not // accessed. // // 4 Each of the three element sets defines an unperturbed heliocentric // orbit. For a given epoch of observation, the position of the body // in its orbit can be predicted from these elements, which are // called "osculating elements", using standard two-body analytical // solutions. However, due to planetary perturbations, a given set // of osculating elements remains usable for only as long as the // unperturbed orbit that it describes is an adequate approximation // to reality. Attached to such a set of elements is a date called // the "osculating epoch", at which the elements are, momentarily, // a perfect representation of the instantaneous position and // velocity of the body. // // Therefore, for any given problem there are up to three different // epochs in play, and it is vital to distinguish clearly between // them: // // . The epoch of observation: the moment in time for which the // position of the body is to be predicted. // // . The epoch defining the position of the body: the moment in time // at which, in the absence of purturbations, the specified // position (mean longitude, mean anomaly, or perihelion) is // reached. // // . The osculating epoch: the moment in time at which the given // elements are correct. // // For the major-planet and minor-planet cases it is usual to make // the epoch that defines the position of the body the same as the // epoch of osculation. Thus, only two different epochs are // involved: the epoch of the elements and the epoch of observation. // // For comets, the epoch of perihelion fixes the position in the // orbit and in general a different epoch of osculation will be // chosen. Thus, all three types of epoch are involved. // // For the present routine: // // . The epoch of observation is the argument DATE. // // . The epoch defining the position of the body is the argument // EPOCH. // // . The osculating epoch is not used and is assumed to be close // enough to the epoch of observation to deliver adequate accuracy. // If not, a preliminary call to sla_PERTEL may be used to update // the element-set (and its associated osculating epoch) by // applying planetary perturbations. // // 5 The reference frame for the result is with respect to the mean // equator and equinox of epoch J2000. // // 6 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, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983. // // P.T.Wallace Starlink 31 December 2002 // // Copyright (C) 2002 Rutherford Appleton Laboratory // struct planel { double PV[6];int JSTAT; planel( const double DATE,const int JFORM,const double EPOCH,const double ORBINC,const double ANODE,const double PERIH,const double AORQ,const double E,const double AORL,const double DM ) { sla_planel_( &DATE,&JFORM,&EPOCH,&ORBINC,&ANODE,&PERIH,&AORQ,&E,&AORL,&DM,PV,&JSTAT ); } std::string __str__() const { std::ostringstream s; s << "PV,JSTAT : " << PV<<" "<(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 aoppat& o ) { return s << o.__str__(); } // - - - - - - // C A F 2 R // - - - - - - // // Convert degrees, arcminutes, arcseconds to radians // (single precision) // // Given: // IDEG int degrees // IAMIN int arcminutes // ASEC real arcseconds // // Returned: // RAD real angle in radians // J int status: 0 = OK // 1 = IDEG outside range 0-359 // 2 = IAMIN outside range 0-59 // 3 = ASEC outside range 0-59.999... // // Notes: // // 1) The result is computed even if any of the range checks // fail. // // 2) The sign must be dealt with outside this routine. // // P.T.Wallace Starlink 23 August 1996 // // Copyright (C) 1996 Rutherford Appleton Laboratory // struct caf2r { float RAD;int J; caf2r( const int IDEG,const int IAMIN,const float ASEC ) { sla_caf2r_( &IDEG,&IAMIN,&ASEC,&RAD,&J ); } std::string __str__() const { std::ostringstream s; s << "RAD,J : " << RAD<<" "< 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<<" "< 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__(); } // - - - - - // S 2 T P // - - - - - // // Projection of spherical coordinates onto tangent plane: // "gnomonic" projection - "standard coordinates" // (single precision) // // Given: // RA,DEC real spherical coordinates of point to be projected // RAZ,DECZ real spherical coordinates of tangent point // // Returned: // XI,ETA real rectangular coordinates on tangent plane // J int status: 0 = OK, star on tangent plane // 1 = error, star too far from axis // 2 = error, antistar on tangent plane // 3 = error, antistar too far from axis // // P.T.Wallace Starlink 18 July 1996 // // Copyright (C) 1996 Rutherford Appleton Laboratory // struct s2tp { float XI;float ETA;int J; s2tp( const float RA,const float DEC,const float RAZ,const float DECZ ) { sla_s2tp_( &RA,&DEC,&RAZ,&DECZ,&XI,&ETA,&J ); } std::string __str__() const { std::ostringstream s; s << "XI,ETA,J : " << XI<<" "<(STRING.length()),&IPTR,&A,&J ); } std::string __str__() const { std::ostringstream s; s << "A,J : " << A<<" "<(STRING.length()),&IPTR,&A,&J ); } std::string __str__() const { std::ostringstream s; s << "A,J : " << A<<" "< 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<<" "< 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<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const pcd& o ) { return s << o.__str__(); } // - - - - // D T T // - - - - // // Increment to be applied to Coordinated Universal Time UTC to give // Terrestrial Time TT (formerly Ephemeris Time ET) // // (double precision) // // Given: // UTC d UTC date as a modified JD (JD-2400000.5) // // Result: TT-UTC in seconds // // Notes: // // 1 The UTC is specified to be a date rather than a time to indicate // that care needs to be taken not to specify an instant which lies // within a leap second. Though in most cases UTC can include the // fractional part, correct behaviour on the day of a leap second // can only be guaranteed up to the end of the second 23:59:59. // // 2 Pre 1972 January 1 a fixed value of 10 + ET-TAI is returned. // // 3 See also the routine sla_DT, which roughly estimates ET-UT for // historical epochs. // // Called: sla_DAT // // P.T.Wallace Starlink 6 December 1994 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline double dtt( const double UTC ) { return sla_dtt_( &UTC ); }; // - - - - // E P B // - - - - // // Conversion of Modified Julian Date to Besselian Epoch // (double precision) // // Given: // DATE dp Modified Julian Date (JD - 2400000.5) // // The result is the Besselian Epoch. // // Reference: // Lieske,J.H., 1979. Astron.Astrophys.,73,282. // // P.T.Wallace Starlink February 1984 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline double epb( const double DATE ) { return sla_epb_( &DATE ); }; // - - - - - - // D A V 2 M // - - - - - - // // Form the rotation matrix corresponding to a given axial vector. // (double 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 d(3) axial vector (radians) // // Returned: // RMAT d(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 dav2m { double RMAT[3][3]; dav2m( const double AXVEC[3] ) { sla_dav2m_( AXVEC,RMAT ); } std::string __str__() const { std::ostringstream s; s << "RMAT : " << RMAT << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const dav2m& o ) { return s << o.__str__(); } // - - - - - // G E 5 0 // - - - - - // // Transformation from IAU 1958 galactic coordinates to // B1950.0 'FK4' equatorial coordinates (double precision) // // Given: // DL,DB dp galactic longitude and latitude L2,B2 // // Returned: // DR,DD dp B1950.0 'FK4' RA,Dec // // (all arguments are radians) // // Called: // sla_DCS2C, sla_DIMXV, sla_DCC2S, sla_ADDET, sla_DRANRM, sla_DRANGE // // Note: // The equatorial coordinates are B1950.0 'FK4'. Use the // routine sla_GALEQ if conversion to J2000.0 coordinates // is required. // // Reference: // Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960) // // P.T.Wallace Starlink 5 September 1993 // // Copyright (C) 1995 Rutherford Appleton Laboratory // struct ge50 { double DR;double DD; ge50( const double DL,const double DB ) { sla_ge50_( &DL,&DB,&DR,&DD ); } std::string __str__() const { std::ostringstream s; s << "DR,DD : " << DR<<" "<
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 euler& o ) { return s << o.__str__(); } // - - - - - // R V L G // - - - - - // // Velocity component in a given direction due to the combination // of the rotation of the Galaxy and the motion of the Galaxy // relative to the mean motion of the local group (single precision) // // Given: // R2000,D2000 real J2000.0 mean RA,Dec (radians) // // Result: // Component of SOLAR motion in direction R2000,D2000 (km/s) // // Sign convention: // The result is +ve when the Sun is receding from the // given point on the sky. // // Reference: // IAU Trans 1976, 168, p201. // // Called: // sla_CS2C, sla_VDV // // P.T.Wallace Starlink June 1985 // // Copyright (C) 1995 Rutherford Appleton Laboratory // inline float rvlg( const float R2000,const float D2000 ) { return sla_rvlg_( &R2000,&D2000 ); }; // - - - - - - // I D C H F // - - - - - - // // Internal routine used by DFLTIN // // Identify next character in string // // Given: // STRING char string // NPTR int pointer to character to be identified // // Returned: // NPTR int incremented unless end of field // NVEC int vector for identified character // NDIGIT int 0-9 if character was a numeral // DIGIT double equivalent of NDIGIT // // NVEC takes the following values: // // 1 0-9 // 2 space or TAB !!! n.b. ASCII TAB assumed !!! // 3 D,d,E or e // 4 . // 5 + // 6 - // 7 , // 8 else // 9 outside field // // If the character is not 0-9, NDIGIT and DIGIT are either not // altered or are set to arbitrary values. // // P.T.Wallace Starlink 22 December 1992 // // Copyright (C) 1995 Rutherford Appleton Laboratory // struct _idchf { int NVEC;int NDIGIT;double DIGIT; _idchf( const char STRING,int NPTR ) { sla__idchf_( &STRING,&NPTR,&NVEC,&NDIGIT,&DIGIT ); } std::string __str__() const { std::ostringstream s; s << "NVEC,NDIGIT,DIGIT : " << NVEC<<" "< 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__(); } // - - - - - - - // P R E N U T // - - - - - - - // // Form the matrix of precession and nutation (SF2001) // (double precision) // // Given: // EPOCH dp Julian Epoch for mean coordinates // DATE dp Modified Julian Date (JD-2400000.5) // for true coordinates // // Returned: // RMATPN dp(3,3) combined precession/nutation matrix // // Called: sla_PREC, sla_EPJ, sla_NUT, sla_DMXM // // Notes: // // 1) The epoch and date are TDB (loosely ET). TT will do, or even // UTC. // // 2) The matrix is in the sense V(true) = RMATPN * V(mean) // // Last revision: 3 December 2005 // // Copyright P.T.Wallace. All rights reserved. // struct prenut { double RMATPN[3][3]; prenut( const double EPOCH,const double DATE ) { sla_prenut_( &EPOCH,&DATE,RMATPN ); } std::string __str__() const { std::ostringstream s; s << "RMATPN : " << RMATPN << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const prenut& o ) { return s << o.__str__(); } // - - - - - - - // P R E C E S // - - - - - - - // // Precession - either FK4 (Bessel-Newcomb, pre IAU 1976) or // FK5 (Fricke, post IAU 1976) as required. // // Given: // SYSTEM char precession to be applied: 'FK4' or 'FK5' // EP0,EP1 dp starting and ending epoch // RA,DC dp RA,Dec, mean equator & equinox of epoch EP0 // // Returned: // RA,DC dp RA,Dec, mean equator & equinox of epoch EP1 // // Called: sla_DRANRM, sla_PREBN, sla_PREC, sla_DCS2C, // sla_DMXV, sla_DCC2S // // Notes: // // 1) Lowercase characters in SYSTEM are acceptable. // // 2) The epochs are Besselian if SYSTEM='FK4' and Julian if 'FK5'. // For example, to precess coordinates in the old system from // equinox 1900.0 to 1950.0 the call would be: // CALL sla_PRECES ('FK4', 1900D0, 1950D0, RA, DC) // // 3) This routine will NOT correctly convert between the old and // the new systems - for example conversion from B1950 to J2000. // For these purposes see sla_FK425, sla_FK524, sla_FK45Z and // sla_FK54Z. // // 4) If an invalid SYSTEM is supplied, values of -99D0,-99D0 will // be returned for both RA and DC. // // P.T.Wallace Starlink 20 April 1990 // // Copyright (C) 1995 Rutherford Appleton Laboratory // struct preces { ; preces( const char SYSTEM,const double EP0,const double EP1,double RA,double DC ) { sla_preces_( &SYSTEM,&EP0,&EP1,&RA,&DC ); } std::string __str__() const { std::ostringstream s; s << " : " << "" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const preces& o ) { return s << o.__str__(); } // - - - - - // D J C L // - - - - - // // Modified Julian Date to Gregorian year, month, day, // and fraction of a day. // // Given: // DJM dp modified Julian Date (JD-2400000.5) // // Returned: // IY int year // IM int month // ID int day // FD dp fraction of day // J int status: // 0 = OK // -1 = unacceptable date (before 4701BC March 1) // // The algorithm is adapted from Hatcher 1984 (QJRAS 25, 53-55). // // Last revision: 22 July 2004 // // Copyright P.T.Wallace. All rights reserved. // struct djcl { int IY;int IM;int ID;double FD;int J; djcl( const double DJM ) { sla_djcl_( &DJM,&IY,&IM,&ID,&FD,&J ); } std::string __str__() const { std::ostringstream s; s << "IY,IM,ID,FD,J : " << IY<<" "<" << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const unpcd& o ) { return s << o.__str__(); } // - - - - - - // P R E B N // - - - - - - // // Generate the matrix of precession between two epochs, // using the old, pre-IAU1976, Bessel-Newcomb model, using // Kinoshita's formulation (double precision) // // Given: // BEP0 dp beginning Besselian epoch // BEP1 dp ending Besselian epoch // // Returned: // RMATP dp(3,3) precession matrix // // The matrix is in the sense V(BEP1) = RMATP * V(BEP0) // // Reference: // Kinoshita, H. (1975) 'Formulas for precession', SAO Special // Report No. 364, Smithsonian Institution Astrophysical // Observatory, Cambridge, Massachusetts. // // Called: sla_DEULER // // P.T.Wallace Starlink 23 August 1996 // // Copyright (C) 1996 Rutherford Appleton Laboratory // struct prebn { double RMATP[3][3]; prebn( const double BEP0,const double BEP1 ) { sla_prebn_( &BEP0,&BEP1,RMATP ); } std::string __str__() const { std::ostringstream s; s << "RMATP : " << RMATP << endl; return s.str(); } }; inline ostream& operator<<(ostream& s, const prebn& o ) { return s << o.__str__(); } // - - - - // D A T // - - - - // // Increment to be applied to Coordinated Universal Time UTC to give // International Atomic Time TAI (double precision) // // Given: // UTC d UTC date as a modified JD (JD-2400000.5) // // Result: TAI-UTC in seconds // // Notes: // // 1 The UTC is specified to be a date rather than a time to indicate // that care needs to be taken not to specify an instant which lies // within a leap second. Though in most cases UTC can include the // fractional part, correct behaviour on the day of a leap second // can only be guaranteed up to the end of the second 23:59:59. // // 2 For epochs from 1961 January 1 onwards, the expressions from the // file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used. // // 3 The 5ms time step at 1961 January 1 is taken from 2.58.1 (p87) of // the 1992 Explanatory Supplement. // // 4 UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper // to call the routine with an earlier epoch. However, if this // is attempted, the TAI-UTC expression for the year 1960 is used. // // // :-----------------------------------------: // : : // : IMPORTANT : // : : // : This routine must be updated on each : // : occasion that a leap second is : // : announced : // : : // : Latest leap second: 2016 December 31 : // : : // :-----------------------------------------: // // Last revision: 03 October 2016 // // Copyright P.T.Wallace. All rights reserved. // inline double dat( const double UTC ) { return sla_dat_( &UTC ); }; // - - - - - - // F K 4 5 Z // - - - - - - // // Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero // proper motion in the FK5 frame (double precision) // // This routine converts stars from the old, Bessel-Newcomb, FK4 // system to the new, IAU 1976, FK5, Fricke system, in such a // way that the FK5 proper motion is zero. Because such a star // has, in general, a non-zero proper motion in the FK4 system, // the routine requires the epoch at which the position in the // FK4 system was determined. // // The method is from Appendix 2 of Ref 1, but using the constants // of Ref 4. // // Given: // R1950,D1950 dp B1950.0 FK4 RA,Dec at epoch (rad) // BEPOCH dp Besselian epoch (e.g. 1979.3D0) // // Returned: // R2000,D2000 dp J2000.0 FK5 RA,Dec (rad) // // Notes: // // 1) The epoch BEPOCH is strictly speaking Besselian, but // if a Julian epoch is supplied the result will be // affected only to a negligible extent. // // 2) Conversion from Besselian epoch 1950.0 to Julian epoch // 2000.0 only is provided for. Conversions involving other // epochs will require use of the appropriate precession, // proper motion, and E-terms routines before and/or // after FK45Z is called. // // 3) In the FK4 catalogue the proper motions of stars within // 10 degrees of the poles do not embody the differential // E-term effect and should, strictly speaking, be handled // in a different manner from stars outside these regions. // However, given the general lack of homogeneity of the star // data available for routine astrometry, the difficulties of // handling positions that may have been determined from // astrometric fields spanning the polar and non-polar regions, // the likelihood that the differential E-terms effect was not // taken into account when allowing for proper motion in past // astrometry, and the undesirability of a discontinuity in // the algorithm, the decision has been made in this routine to // include the effect of differential E-terms on the proper // motions for all stars, whether polar or not. At epoch 2000, // and measuring on the sky rather than in terms of dRA, the // errors resulting from this simplification are less than // 1 milliarcsecond in position and 1 milliarcsecond per // century in proper motion. // // References: // // 1 Aoki,S., et al, 1983. Astron.Astrophys., 128, 263. // // 2 Smith, C.A. et al, 1989. "The transformation of astrometric // catalog systems to the equinox J2000.0". Astron.J. 97, 265. // // 3 Yallop, B.D. et al, 1989. "Transformation of mean star places // from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space". // Astron.J. 97, 274. // // 4 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to // the Astronomical Almanac", ISBN 0-935702-68-7. // // Called: sla_DCS2C, sla_EPJ, sla_EPB2D, sla_DCC2S, sla_DRANRM // // P.T.Wallace Starlink 21 September 1998 // // Copyright (C) 1998 Rutherford Appleton Laboratory // struct fk45z { double R2000;double D2000; fk45z( const double R1950,const double D1950,const double BEPOCH ) { sla_fk45z_( &R1950,&D1950,&BEPOCH,&R2000,&D2000 ); } std::string __str__() const { std::ostringstream s; s << "R2000,D2000 : " << R2000<<" "<(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<<" "<