Changeset 4563 for trunk/MagicSoft/Mars

Timestamp:

08/10/04 16:03:01 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (3 diffs)
• NEWS (modified) (1 diff)
• macros/sql/filldotrun.C (modified) (1 diff)
• macros/sql/readrbk.C (modified) (1 diff)
• mastro/MAstro.cc (modified) (1 diff)
• mastro/MAstro.h (modified) (2 diffs)
• mbase/MTime.cc (modified) (2 diffs)
• mbase/MTime.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -24 +24 @@
    * 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 +42 @@
      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 +66 @@
      - set the possibility to have the datacheck display with the 
        function SetDataCheckDisplay()
+
+
 
   2004/08/09: Wolfgang Wittek

trunk/MagicSoft/Mars/NEWS

@@ -2 +2 @@
  *** Version <cvs>
 
-   - new macros for database interaction: filldotrun.C
+   - new macros for database interaction: filldotrun.C, filldotrbk.C
 
 

trunk/MagicSoft/Mars/macros/sql/filldotrun.C

@@ -45 +45 @@
 // 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
 //

trunk/MagicSoft/Mars/macros/sql/readrbk.C

@@ -173 +173 @@
 // 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);

trunk/MagicSoft/Mars/mastro/MAstro.cc

@@ -405 +405 @@
     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;
+}

trunk/MagicSoft/Mars/mastro/MAstro.h

@@ -11 +11 @@
     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 +56 @@
     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)
 };

trunk/MagicSoft/Mars/mbase/MTime.cc

@@ -113 +113 @@
 // --------------------------------------------------------------------------
 //
+// 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 +237 @@
 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;

trunk/MagicSoft/Mars/mbase/MTime.h

@@ -92 +92 @@
     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;