Changeset 4563


Ignore:
Timestamp:
08/10/04 16:03:01 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r4562 r4563  
    2424   * mbadpixels/MBadPixelsTreat.[h,cc]
    2525     - replaced flag SetSloppyTreatment() with SetHardTreatment()
    26         Now the default behaviour consists on treating only the
    27         unsuitable pixels, and not also the ureliable, as it was
    28         before. If you want to keep on treating both unsuitable and
    29         unreliable pixels, please set the new flag on in your macro.
     26       Now the default behaviour consists on treating only the
     27       unsuitable pixels, and not also the ureliable, as it was
     28       before. If you want to keep on treating both unsuitable and
     29       unreliable pixels, please set the new flag on in your macro.
    3030
    3131   * mjobs/MJExtractCalibTest.cc
    3232     - removed line containing SetSloppyTreatment()
     33
    3334
    3435
     
    4142     mpedestal/MPedCalcFromLoGain.[h,cc]:
    4243     - removed dependancy on MBadPixelsCam
     44
     45   * mastro/MAstro.[h,cc]:
     46     - implemented GetMoonPeriod
     47     - implemented GetMoonPhase
     48
     49   * mbase/MTime.cc:
     50     - Fixed a floating point problem in Set(&tv)
     51     - added GetDateOfSunrise
     52
    4353
    4454
     
    5666     - set the possibility to have the datacheck display with the
    5767       function SetDataCheckDisplay()
     68
     69
    5870
    5971  2004/08/09: Wolfgang Wittek
  • trunk/MagicSoft/Mars/NEWS

    r4521 r4563  
    22 *** Version <cvs>
    33
    4    - new macros for database interaction: filldotrun.C
     4   - new macros for database interaction: filldotrun.C, filldotrbk.C
    55
    66
  • trunk/MagicSoft/Mars/macros/sql/filldotrun.C

    r4516 r4563  
    4545// info will be put into the DB, eg:
    4646//   "/data/MAGIC"                              would do it for all data
    47 //   "/data/MAGIC/Period019"                    would do it for one Period
     47//   "/data/MAGIC/Period019/ccdata"             would do it for one Period
    4848//   "/data/MAGIC/Period019/ccdata/2004_05_17"  would do it for a single day
    4949//
  • trunk/MagicSoft/Mars/macros/sql/readrbk.C

    r4513 r4563  
    173173// loop over all files in this path
    174174//
    175 void readrbk(TString path="/data/MAGIC/Period014")
     175void readrbk(TString path="/data/MAGIC/Period017/ccdata")
    176176{
    177177    MDirIter Next(path, "*.rbk", -1);
  • trunk/MagicSoft/Mars/mastro/MAstro.cc

    r3666 r4563  
    405405    return TMath::ATan2(cosl*sinp, snlt-cslt*cosp);
    406406}
     407
     408
     409// --------------------------------------------------------------------------
     410//
     411// Kepler - solve the equation of Kepler
     412//
     413Double_t MAstro::Kepler(Double_t m, Double_t ecc)
     414{
     415    m *= TMath::DegToRad();
     416
     417    Double_t delta = 0;
     418    Double_t e = m;
     419    do {
     420        delta = e - ecc * sin(e) - m;
     421        e -= delta / (1 - ecc * cos(e));
     422    } while (fabs(delta) > 1e-6);
     423
     424    return e;
     425}
     426
     427// --------------------------------------------------------------------------
     428//
     429// GetMoonPhase - calculate phase of moon as a fraction:
     430//  Returns -1 if calculation failed
     431//
     432Double_t MAstro::GetMoonPhase(Double_t mjd)
     433{
     434    /****** Calculation of the Sun's position. ******/
     435
     436    // date within epoch
     437    const Double_t epoch = 44238;       // 1980 January 0.0
     438    const Double_t day = mjd - epoch;
     439    if (day<0)
     440    {
     441        cout << "MAstro::GetMoonPhase - Day before Jan 1980" << endl;
     442        return -1;
     443    }
     444
     445    // mean anomaly of the Sun
     446    const Double_t n = fmod(day*360/365.2422, 360);
     447
     448    const Double_t elonge = 278.833540; // ecliptic longitude of the Sun at epoch 1980.0
     449    const Double_t elongp = 282.596403; // ecliptic longitude of the Sun at perigee
     450
     451    // convert from perigee co-ordinates to epoch 1980.0
     452    const Double_t m = fmod(n + elonge - elongp + 360, 360);
     453
     454    // solve equation of Kepler
     455    const Double_t eccent = 0.016718; // eccentricity of Earth's orbit
     456    const Double_t k   = Kepler(m, eccent);
     457    const Double_t ec0 = sqrt((1 + eccent) / (1 - eccent)) * tan(k / 2);
     458    // true anomaly
     459    const Double_t ec  = 2 * atan(ec0) * TMath::RadToDeg();
     460
     461    // Sun's geocentric ecliptic longitude
     462    const Double_t lambdasun = fmod(ec + elongp + 720, 360);
     463
     464
     465    /****** Calculation of the Moon's position. ******/
     466
     467    // Moon's mean longitude.
     468    const Double_t mmlong  = 64.975464;  // moon's mean lonigitude at the epoch
     469    const Double_t ml      = fmod(13.1763966*day + mmlong + 360, 360);
     470    // Moon's mean anomaly.
     471    const Double_t mmlongp = 349.383063; // mean longitude of the perigee at the epoch
     472    const Double_t mm      = fmod(ml - 0.1114041*day - mmlongp + 720, 360);
     473    // Evection.
     474    const Double_t ev   = 1.2739 * sin((2 * (ml - lambdasun) - mm)*TMath::DegToRad());
     475    // Annual equation.
     476    const Double_t sinm = TMath::Sin(m*TMath::DegToRad());
     477    const Double_t ae   = 0.1858 * sinm;
     478    // Correction term.
     479    const Double_t a3   = 0.37 * sinm;
     480    // Corrected anomaly.
     481    const Double_t mmp  = (mm + ev - ae - a3)*TMath::DegToRad();
     482    // Correction for the equation of the centre.
     483    const Double_t mec  = 6.2886 * sin(mmp);
     484    // Another correction term.
     485    const Double_t a4   = 0.214 * sin(2 * mmp);
     486    // Corrected longitude.
     487    const Double_t lp   = ml + ev + mec - ae + a4;
     488    // Variation.
     489    const Double_t v    = 0.6583 * sin(2 * (lp - lambdasun)*TMath::DegToRad());
     490    // True longitude.
     491    const Double_t lpp  = lp + v;
     492    // Age of the Moon in degrees.
     493    const Double_t age  = (lpp - lambdasun)*TMath::DegToRad();
     494
     495    // Calculation of the phase of the Moon.
     496    return (1 - TMath::Cos(age)) / 2;
     497}
     498
     499// --------------------------------------------------------------------------
     500//
     501// Calculate the Period to which the time belongs to. The Period is defined
     502// as the number of synodic months ellapsed since the first full moon
     503// after Jan 1st 1980 (which was @ MJD=44240.37917)
     504//
     505// In a addition to this calculation an offset of 284 periods is substracted.
     506// Taking this offset into account you get the Period as it is defined
     507// for the MAGIC observations. If you want to get the MAGIC observation
     508// period number make sure that your Mjd is a integer without a floating
     509// point reminder.
     510//
     511// For MAGIC observation period do eg:
     512//   GetMoonPeriod(53257)
     513// or
     514//   MTime t;
     515//   t.SetMjd(53257);
     516//   GetMoonPeriod(t.GetMjd());
     517// or
     518//   MTime t;
     519//   t.Set(2004, 1, 1, 0);
     520//   GetMoonPeriod(t.GetMjd());
     521// or
     522//   MTime t;
     523//   t.Now();
     524//   GetMoonPeriod(TMath::Floor(t.GetMjd()));
     525//
     526Double_t MAstro::GetMoonPeriod(Double_t mjd, Int_t offset)
     527{
     528    const Double_t synmonth = 29.53058868; // synodic month (new Moon to new Moon)
     529    const Double_t epoch0   = 44240.37917; // First full moon after 1980/1/1
     530
     531    const Double_t et = mjd-epoch0; // Ellapsed time
     532    return TMath::Floor(et/synmonth) + offset;
     533}
  • trunk/MagicSoft/Mars/mastro/MAstro.h

    r3568 r4563  
    1111    static Double_t Round(Double_t val);
    1212    static Double_t Trunc(Double_t val);
     13
     14    static Double_t Kepler(Double_t m, Double_t ecc);
    1315
    1416public:
     
    5456    static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi);
    5557
     58    static Double_t GetMoonPhase(Double_t mjd);
     59    static Double_t GetMoonPeriod(Double_t mjd, Int_t offset=-284);
     60
    5661    ClassDef(MAstro, 0)
    5762};
  • trunk/MagicSoft/Mars/mbase/MTime.cc

    r4452 r4563  
    113113// --------------------------------------------------------------------------
    114114//
     115// Return date as year(y), month(m), day(d). If the time is afternoon
     116// (>=13:00:00) the date of the next day is returned.
     117//
     118void MTime::GetDateOfSunrise(UShort_t &y, Byte_t &m, Byte_t &d) const
     119{
     120    MAstro::Mjd2Ymd(fMjd, y, m, d);
     121}
     122
     123// --------------------------------------------------------------------------
     124//
    115125// Return the time in the range [0h, 24h) = [0h0m0.000s - 23h59m59.999s]
    116126//
     
    227237void MTime::Set(const struct timeval &tv)
    228238{
    229     const UInt_t mjd = 1000*tv.tv_sec/kDay + 40587;
     239    const UInt_t mjd = (UInt_t)TMath::Floor(1000.*tv.tv_sec/kDay) + 40587;
     240
    230241    const Long_t tm  = tv.tv_sec%(24*3600)*1000 + tv.tv_usec/1000;
    231242    const UInt_t ms  = tv.tv_usec%1000;
  • trunk/MagicSoft/Mars/mbase/MTime.h

    r3967 r4563  
    9292    TString  GetFileName() const;
    9393    void     GetDate(UShort_t &y, Byte_t &m, Byte_t &d) const;
     94    void     GetDateOfSunrise(UShort_t &y, Byte_t &m, Byte_t &d) const;
    9495    TTime    GetRootTime() const;
    9596    Double_t GetAxisTime() const;
Note: See TracChangeset for help on using the changeset viewer.