Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3370) +++ trunk/MagicSoft/Mars/Changelog (revision 3371) @@ -9,4 +9,6 @@ - now the arrival time is at the half max of the spline: half max = (maximum - pedestal)/2 + + 2004/03/01: Wolfgang Wittek @@ -28,7 +30,12 @@ - fixed a bug in determin the file-type + * mbase/MTime.[h,cc]: + - added GetGmst() member function + * mastro/MAstro.[h,cc]: - added Coordinate2Angle - added AngularDistance + - implemented UT2GMST + - replaced defined by TMath::Pi() * mastro/MObservatory.[h,cc]: Index: trunk/MagicSoft/Mars/mastro/MAstro.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 3370) +++ trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 3371) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 11/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -33,4 +33,6 @@ #include // TVector3 +#include "MTime.h" // MTime::GetGmst + ClassImp(MAstro); @@ -56,6 +58,6 @@ { /* pi/(180*3600): arcseconds to radians */ -#define DAS2R 4.8481368110953599358991410235794797595635330237270e-6 - return Hms2Sec(deg, min, sec, sgn)*DAS2R; + //#define DAS2R 4.8481368110953599358991410235794797595635330237270e-6 + return Hms2Sec(deg, min, sec, sgn)*TMath::Pi()/(180*3600)/**DAS2R*/; } @@ -63,6 +65,6 @@ { /* pi/(12*3600): seconds of time to radians */ -#define DS2R 7.2722052166430399038487115353692196393452995355905e-5 - return Hms2Sec(hor, min, sec, sgn)*DS2R; +//#define DS2R 7.2722052166430399038487115353692196393452995355905e-5 + return Hms2Sec(hor, min, sec, sgn)*TMath::Pi()/(12*3600)/**DS2R*/; } @@ -295,2 +297,11 @@ return v0.Angle(v1); } + +// -------------------------------------------------------------------------- +// +// Calls MTime::GetGmst() Better use MTime::GetGmst() directly +// +Double_t MAstro::UT2GMST(Double_t ut1) +{ + return MTime(ut1).GetGmst(); +} Index: trunk/MagicSoft/Mars/mastro/MAstro.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstro.h (revision 3370) +++ trunk/MagicSoft/Mars/mastro/MAstro.h (revision 3371) @@ -13,4 +13,5 @@ public: + // Angle/Time conversion functions static Double_t Hms2Sec(Int_t deg, UInt_t min, Double_t sec, char sgn='+'); static Double_t Dms2Rad(Int_t deg, UInt_t min, Double_t sec, Char_t sgn='+'); @@ -37,11 +38,15 @@ static void Hor2Hm(Double_t rad, Char_t &sgn, UShort_t °, Double_t &min); + // Angle treatment functions static Bool_t String2Angle(TString &str, Double_t &ret); static Bool_t Coordinate2Angle(const TString &str, Double_t &ret); static Double_t AngularDistance(Double_t theta0, Double_t phi0, Double_t theta1, Double_t phi1); - + + // Time treatment functions static void Mjd2Ymd(UInt_t mjd, UShort_t &y, Byte_t &m, Byte_t &d); static Int_t Ymd2Mjd(UShort_t y, Byte_t m, Byte_t d); + + static Double_t UT2GMST(Double_t ut1); ClassDef(MAstro, 0) Index: trunk/MagicSoft/Mars/mbase/MTime.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MTime.cc (revision 3370) +++ trunk/MagicSoft/Mars/mbase/MTime.cc (revision 3371) @@ -317,4 +317,50 @@ return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Conversion from Universal Time to Greenwich mean sidereal time, +// with rounding errors minimized. +// +// The result is the Greenwich Mean Sidereal Time (radians) +// +// There is no restriction on how the UT is apportioned between the +// date and ut1 arguments. Either of the two arguments could, for +// example, be zero and the entire date+time supplied in the other. +// However, the routine is designed to deliver maximum accuracy when +// the date argument is a whole number and the ut argument lies in +// the range 0 to 1, or vice versa. +// +// The algorithm is based on the IAU 1982 expression (see page S15 of +// the 1984 Astronomical Almanac). This is always described as giving +// the GMST at 0 hours UT1. In fact, it gives the difference between +// the GMST and the UT, the steady 4-minutes-per-day drawing-ahead of +// ST with respect to UT. When whole days are ignored, the expression +// happens to equal the GMST at 0 hours UT1 each day. +// +// In this routine, the entire UT1 (the sum of the two arguments date +// and ut) is used directly as the argument for the standard formula. +// The UT1 is then added, but omitting whole days to conserve accuracy. +// +// The extra numerical precision delivered by the present routine is +// unlikely to be important in an absolute sense, but may be useful +// when critically comparing algorithms and in applications where two +// sidereal times close together are differenced. +// +Double_t MTime::GetGmst() const +{ + const Double_t ut = (Double_t)(fNanoSec/1e6+(Long_t)fTime)/kDay; + + // Julian centuries since J2000. + const Double_t t = (ut -(51544.5-fMjd)) / 36525.0; + + // GMST at this UT1 + const Double_t r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t; + const Double_t r2 = 86400.0*ut; + + const Double_t sum = (r1+r2)/(24*3600); + + return fmod(sum, 1)*TMath::TwoPi()+TMath::TwoPi(); } Index: trunk/MagicSoft/Mars/mbase/MTime.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MTime.h (revision 3370) +++ trunk/MagicSoft/Mars/mbase/MTime.h (revision 3371) @@ -56,5 +56,5 @@ Set(tm); } - MTime(Double_t d); + MTime(Double_t mjd); MTime(const MTime& t) : fMjd(t.fMjd), fTime(t.fTime), fNanoSec(t.fNanoSec) { @@ -85,4 +85,5 @@ void SetMjd(Double_t m); Double_t GetMjd() const; + Double_t GetGmst() const; TString GetString() const; TString GetStringFmt(const char *fmt=0) const;