Index: trunk/MagicSoft/Mars/mastro/MAstro.cc
===================================================================
--- trunk/MagicSoft/Mars/mastro/MAstro.cc	(revision 4521)
+++ 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 4521)
+++ 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)
 };