Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4562) +++ /trunk/MagicSoft/Mars/Changelog (revision 4563) @@ -24,11 +24,12 @@ * mbadpixels/MBadPixelsTreat.[h,cc] - replaced flag SetSloppyTreatment() with SetHardTreatment() - Now the default behaviour consists on treating only the - unsuitable pixels, and not also the ureliable, as it was - before. If you want to keep on treating both unsuitable and - unreliable pixels, please set the new flag on in your macro. + Now the default behaviour consists on treating only the + unsuitable pixels, and not also the ureliable, as it was + before. If you want to keep on treating both unsuitable and + unreliable pixels, please set the new flag on in your macro. * mjobs/MJExtractCalibTest.cc - removed line containing SetSloppyTreatment() + @@ -41,4 +42,13 @@ mpedestal/MPedCalcFromLoGain.[h,cc]: - removed dependancy on MBadPixelsCam + + * mastro/MAstro.[h,cc]: + - implemented GetMoonPeriod + - implemented GetMoonPhase + + * mbase/MTime.cc: + - Fixed a floating point problem in Set(&tv) + - added GetDateOfSunrise + @@ -56,4 +66,6 @@ - set the possibility to have the datacheck display with the function SetDataCheckDisplay() + + 2004/08/09: Wolfgang Wittek Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 4562) +++ /trunk/MagicSoft/Mars/NEWS (revision 4563) @@ -2,5 +2,5 @@ *** Version - - new macros for database interaction: filldotrun.C + - new macros for database interaction: filldotrun.C, filldotrbk.C Index: /trunk/MagicSoft/Mars/macros/sql/filldotrun.C =================================================================== --- /trunk/MagicSoft/Mars/macros/sql/filldotrun.C (revision 4562) +++ /trunk/MagicSoft/Mars/macros/sql/filldotrun.C (revision 4563) @@ -45,5 +45,5 @@ // info will be put into the DB, eg: // "/data/MAGIC" would do it for all data -// "/data/MAGIC/Period019" would do it for one Period +// "/data/MAGIC/Period019/ccdata" would do it for one Period // "/data/MAGIC/Period019/ccdata/2004_05_17" would do it for a single day // Index: /trunk/MagicSoft/Mars/macros/sql/readrbk.C =================================================================== --- /trunk/MagicSoft/Mars/macros/sql/readrbk.C (revision 4562) +++ /trunk/MagicSoft/Mars/macros/sql/readrbk.C (revision 4563) @@ -173,5 +173,5 @@ // loop over all files in this path // -void readrbk(TString path="/data/MAGIC/Period014") +void readrbk(TString path="/data/MAGIC/Period017/ccdata") { MDirIter Next(path, "*.rbk", -1); Index: /trunk/MagicSoft/Mars/mastro/MAstro.cc =================================================================== --- /trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 4562) +++ /trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 4563) @@ -405,2 +405,129 @@ return TMath::ATan2(cosl*sinp, snlt-cslt*cosp); } + + +// -------------------------------------------------------------------------- +// +// Kepler - solve the equation of Kepler +// +Double_t MAstro::Kepler(Double_t m, Double_t ecc) +{ + m *= TMath::DegToRad(); + + Double_t delta = 0; + Double_t e = m; + do { + delta = e - ecc * sin(e) - m; + e -= delta / (1 - ecc * cos(e)); + } while (fabs(delta) > 1e-6); + + return e; +} + +// -------------------------------------------------------------------------- +// +// GetMoonPhase - calculate phase of moon as a fraction: +// Returns -1 if calculation failed +// +Double_t MAstro::GetMoonPhase(Double_t mjd) +{ + /****** Calculation of the Sun's position. ******/ + + // date within epoch + const Double_t epoch = 44238; // 1980 January 0.0 + const Double_t day = mjd - epoch; + if (day<0) + { + cout << "MAstro::GetMoonPhase - Day before Jan 1980" << endl; + return -1; + } + + // mean anomaly of the Sun + const Double_t n = fmod(day*360/365.2422, 360); + + const Double_t elonge = 278.833540; // ecliptic longitude of the Sun at epoch 1980.0 + const Double_t elongp = 282.596403; // ecliptic longitude of the Sun at perigee + + // convert from perigee co-ordinates to epoch 1980.0 + const Double_t m = fmod(n + elonge - elongp + 360, 360); + + // solve equation of Kepler + const Double_t eccent = 0.016718; // eccentricity of Earth's orbit + const Double_t k = Kepler(m, eccent); + const Double_t ec0 = sqrt((1 + eccent) / (1 - eccent)) * tan(k / 2); + // true anomaly + const Double_t ec = 2 * atan(ec0) * TMath::RadToDeg(); + + // Sun's geocentric ecliptic longitude + const Double_t lambdasun = fmod(ec + elongp + 720, 360); + + + /****** Calculation of the Moon's position. ******/ + + // Moon's mean longitude. + const Double_t mmlong = 64.975464; // moon's mean lonigitude at the epoch + const Double_t ml = fmod(13.1763966*day + mmlong + 360, 360); + // Moon's mean anomaly. + const Double_t mmlongp = 349.383063; // mean longitude of the perigee at the epoch + const Double_t mm = fmod(ml - 0.1114041*day - mmlongp + 720, 360); + // Evection. + const Double_t ev = 1.2739 * sin((2 * (ml - lambdasun) - mm)*TMath::DegToRad()); + // Annual equation. + const Double_t sinm = TMath::Sin(m*TMath::DegToRad()); + const Double_t ae = 0.1858 * sinm; + // Correction term. + const Double_t a3 = 0.37 * sinm; + // Corrected anomaly. + const Double_t mmp = (mm + ev - ae - a3)*TMath::DegToRad(); + // Correction for the equation of the centre. + const Double_t mec = 6.2886 * sin(mmp); + // Another correction term. + const Double_t a4 = 0.214 * sin(2 * mmp); + // Corrected longitude. + const Double_t lp = ml + ev + mec - ae + a4; + // Variation. + const Double_t v = 0.6583 * sin(2 * (lp - lambdasun)*TMath::DegToRad()); + // True longitude. + const Double_t lpp = lp + v; + // Age of the Moon in degrees. + const Double_t age = (lpp - lambdasun)*TMath::DegToRad(); + + // Calculation of the phase of the Moon. + return (1 - TMath::Cos(age)) / 2; +} + +// -------------------------------------------------------------------------- +// +// Calculate the Period to which the time belongs to. The Period is defined +// as the number of synodic months ellapsed since the first full moon +// after Jan 1st 1980 (which was @ MJD=44240.37917) +// +// In a addition to this calculation an offset of 284 periods is substracted. +// Taking this offset into account you get the Period as it is defined +// for the MAGIC observations. If you want to get the MAGIC observation +// period number make sure that your Mjd is a integer without a floating +// point reminder. +// +// For MAGIC observation period do eg: +// GetMoonPeriod(53257) +// or +// MTime t; +// t.SetMjd(53257); +// GetMoonPeriod(t.GetMjd()); +// or +// MTime t; +// t.Set(2004, 1, 1, 0); +// GetMoonPeriod(t.GetMjd()); +// or +// MTime t; +// t.Now(); +// GetMoonPeriod(TMath::Floor(t.GetMjd())); +// +Double_t MAstro::GetMoonPeriod(Double_t mjd, Int_t offset) +{ + const Double_t synmonth = 29.53058868; // synodic month (new Moon to new Moon) + const Double_t epoch0 = 44240.37917; // First full moon after 1980/1/1 + + const Double_t et = mjd-epoch0; // Ellapsed time + return TMath::Floor(et/synmonth) + offset; +} Index: /trunk/MagicSoft/Mars/mastro/MAstro.h =================================================================== --- /trunk/MagicSoft/Mars/mastro/MAstro.h (revision 4562) +++ /trunk/MagicSoft/Mars/mastro/MAstro.h (revision 4563) @@ -11,4 +11,6 @@ static Double_t Round(Double_t val); static Double_t Trunc(Double_t val); + + static Double_t Kepler(Double_t m, Double_t ecc); public: @@ -54,4 +56,7 @@ static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi); + static Double_t GetMoonPhase(Double_t mjd); + static Double_t GetMoonPeriod(Double_t mjd, Int_t offset=-284); + ClassDef(MAstro, 0) }; Index: /trunk/MagicSoft/Mars/mbase/MTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4562) +++ /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4563) @@ -113,4 +113,14 @@ // -------------------------------------------------------------------------- // +// Return date as year(y), month(m), day(d). If the time is afternoon +// (>=13:00:00) the date of the next day is returned. +// +void MTime::GetDateOfSunrise(UShort_t &y, Byte_t &m, Byte_t &d) const +{ + MAstro::Mjd2Ymd(fMjd, y, m, d); +} + +// -------------------------------------------------------------------------- +// // Return the time in the range [0h, 24h) = [0h0m0.000s - 23h59m59.999s] // @@ -227,5 +237,6 @@ void MTime::Set(const struct timeval &tv) { - const UInt_t mjd = 1000*tv.tv_sec/kDay + 40587; + const UInt_t mjd = (UInt_t)TMath::Floor(1000.*tv.tv_sec/kDay) + 40587; + const Long_t tm = tv.tv_sec%(24*3600)*1000 + tv.tv_usec/1000; const UInt_t ms = tv.tv_usec%1000; Index: /trunk/MagicSoft/Mars/mbase/MTime.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.h (revision 4562) +++ /trunk/MagicSoft/Mars/mbase/MTime.h (revision 4563) @@ -92,4 +92,5 @@ TString GetFileName() const; void GetDate(UShort_t &y, Byte_t &m, Byte_t &d) const; + void GetDateOfSunrise(UShort_t &y, Byte_t &m, Byte_t &d) const; TTime GetRootTime() const; Double_t GetAxisTime() const;