Changeset 8066 for trunk/MagicSoft/Mars/mastro

Timestamp:

10/14/06 19:20:08 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mastro

Files:

• AstroIncl.h (modified) (1 diff)
• MAstro.cc (modified) (2 diffs)
• MAstro.h (modified) (2 diffs)
• MAstroCatalog.cc (modified) (3 diffs)
• MAstroCatalog.h (modified) (3 diffs)
• MObservatory.cc (modified) (3 diffs)
• MObservatory.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mastro/AstroIncl.h

@@ -2 +2 @@
 
 #include <TArrayC.h>
+#include <TArrayD.h>
 #include "MArrayB.h"
 #include "MArrayS.h"

trunk/MagicSoft/Mars/mastro/MAstro.cc

@@ -33 +33 @@
 #include <iostream>
 
+#include <TArrayD.h>  // TArrayD
 #include <TVector3.h> // TVector3
 
@@ -585 +586 @@
 // --------------------------------------------------------------------------
 //
+// Returns right ascension and declination [rad] of the sun at the
+// given mjd (ra, dec).
+//
+// returns the mean longitude [rad].
+//
+// from http://xoomer.alice.it/vtomezzo/sunriset/formulas/index.html
+//
+Double_t MAstro::GetSunRaDec(Double_t mjd, Double_t &ra, Double_t &dec)
+{
+    const Double_t T = (mjd-51544.5)/36525;// +  (h-12)/24.0;
+
+    const Double_t T2 = T<0 ? -T*T : T*T;
+    const Double_t T3 = T*T*T;
+
+    // Find the ecliptic longitude of the Sun
+
+    // Geometric mean longitude of the Sun
+    const Double_t L = 280.46646 + 36000.76983*T + 0.0003032*T2;
+
+    // mean anomaly of the Sun
+    Double_t g = 357.52911 + 35999.05029*T - 0.0001537*T2;
+    g *= TMath::DegToRad();
+
+    // Longitude of the moon's ascending node
+    Double_t omega = 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000;
+    omega *= TMath::DegToRad();
+
+    const Double_t coso = cos(omega);
+    const Double_t sino = sin(omega);
+
+    // Equation of the center
+    const Double_t C = (1.914602 - 0.004817*T - 0.000014*T2)*sin(g) +
+        (0.019993 - 0.000101*T)*sin(2*g) + 0.000289*sin(3*g);
+
+    // True longitude of the sun
+    const Double_t tlong = L + C;
+
+    // Apperent longitude of the Sun (ecliptic)
+    Double_t lambda = tlong - 0.00569 - 0.00478*sino;
+    lambda *= TMath::DegToRad();
+
+    // Obliquity of the ecliptic
+    Double_t obliq = 23.4392911 - 0.01300416667*T - 0.00000016389*T2 + 0.00000050361*T3 + 0.00255625*coso;
+    obliq *= TMath::DegToRad();
+
+    // Find the RA and DEC of the Sun
+    const Double_t sinl =  sin(lambda);
+
+    ra  = atan2(cos(obliq) * sinl, cos(lambda));
+    dec = asin(sin(obliq) * sinl);
+
+    return L*TMath::DegToRad();
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns right ascension and declination [rad] of the moon at the
+// given mjd (ra, dec).
+//
+void MAstro::GetMoonRaDec(Double_t mjd, Double_t &ra, Double_t &dec)
+{
+    // Mean Moon orbit elements as of 1990.0
+    const Double_t l0 = 318.351648 * TMath::DegToRad();
+    const Double_t P0 =  36.340410 * TMath::DegToRad();
+    const Double_t N0 = 318.510107 * TMath::DegToRad();
+    const Double_t i  =   5.145396 * TMath::DegToRad();
+
+    Double_t sunra, sundec, g;
+    {
+        const Double_t T = (mjd-51544.5)/36525;// +  (h-12)/24.0;
+        const Double_t T2 = T<0 ? -T*T : T*T;
+
+        GetSunRaDec(mjd, sunra, sundec);
+
+        // mean anomaly of the Sun
+        g = 357.52911 + 35999.05029*T - 0.0001537*T2;
+        g *= TMath::DegToRad();
+    }
+
+    const Double_t sing   = sin(g)*TMath::DegToRad();
+
+    const Double_t D      = (mjd-47891) * TMath::DegToRad();
+    const Double_t l      =    13.1763966*D + l0;
+    const Double_t MMoon  = l  -0.1114041*D - P0; // Moon's mean anomaly M
+    const Double_t N      = N0 -0.0529539*D;      // Moon's mean ascending node longitude
+
+    const Double_t C      = l-sunra;
+    const Double_t Ev     = 1.2739 * sin(2*C-MMoon) * TMath::DegToRad();
+    const Double_t Ae     = 0.1858 * sing;
+    const Double_t A3     = 0.37   * sing;
+    const Double_t MMoon2 = MMoon+Ev-Ae-A3;  // corrected Moon anomaly
+
+    const Double_t Ec     = 6.2886 * sin(MMoon2)  * TMath::DegToRad();  // equation of centre
+    const Double_t A4     = 0.214  * sin(2*MMoon2)* TMath::DegToRad();
+    const Double_t l2     = l+Ev+Ec-Ae+A4; // corrected Moon's longitude
+
+    const Double_t V      = 0.6583 * sin(2*(l2-sunra)) * TMath::DegToRad();
+    const Double_t l3     = l2+V; // true orbital longitude;
+
+    const Double_t N2     = N -0.16*sing;
+
+    ra  = fmod( N2 + atan2( sin(l3-N2)*cos(i), cos(l3-N2) ), TMath::TwoPi() );
+    dec = asin(sin(l3-N2)*sin(i) );
+}
+
+// --------------------------------------------------------------------------
+//
+// Return Euqation of time in hours for given mjd
+//
+Double_t MAstro::GetEquationOfTime(Double_t mjd)
+{
+    Double_t ra, dec;
+    const Double_t L = fmod(GetSunRaDec(mjd, ra, dec), TMath::TwoPi());
+
+    if (L-ra>TMath::Pi())
+        ra += TMath::TwoPi();
+
+    return 24*(L - ra)/TMath::TwoPi();
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns noon time (the time of the highest altitude of the sun)
+// at the given mjd and at the given observers longitude [deg]
+//
+// The maximum altitude reached at noon time is
+//   altmax = 90.0 + dec - latit;
+//   if (dec > latit)
+//      altmax = 90.0 + latit - dec;
+// dec=Declination of the sun
+//
+Double_t MAstro::GetNoonTime(Double_t mjd, Double_t longit)
+{
+    const Double_t equation = GetEquationOfTime(TMath::Floor(mjd));
+    return 12. + equation - longit/15;
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the time (in hours) between noon (the sun culmination)
+// and the sun being at height alt[deg] (90=zenith, 0=horizont)
+//
+//       civil twilight:      0deg to  -6deg
+//       nautical twilight:  -6deg to -12deg
+//       astronom twilight: -12deg to -18deg
+//
+// latit is the observers latitude in rad
+//
+// returns -1 in case the sun doesn't reach this altitude.
+// (eg. alt=0: Polarnight or -day)
+//
+// To get the sun rise/set:
+//    double timediff = MAstro::GetTimeFromNoonToAlt(mjd, latit*TMath::DegToRad(), par[0]);
+//    double noon     = MAstro::GetNoonTime(mjd, longit);
+//    double N        = TMath::Floor(mjd)+noon/24.;
+//    double risetime = N-timediff/24.;
+//    double settime  = N+timediff/24.;
+//
+Double_t MAstro::GetTimeFromNoonToAlt(Double_t mjd, Double_t latit, Double_t alt)
+{
+    Double_t ra, dec;
+    GetSunRaDec(mjd, ra, dec);
+
+    const Double_t h = alt*TMath::DegToRad();
+
+    const Double_t arg = (sin(h) - sin(latit)*sin(dec))/(cos(latit)*cos(dec));
+
+    return TMath::Abs(arg)>1 ? -1 : 12*acos(arg)/TMath::Pi();
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the time of the sunrise/set calculated before and after
+// the noon of floor(mjd) (TO BE IMPROVED)
+//
+// Being longit and latit the longitude and latitude of the observer
+// in deg and alt the hight above or below the horizont in deg.
+//
+//       civil twilight:      0deg to  -6deg
+//       nautical twilight:  -6deg to -12deg
+//       astronom twilight: -12deg to -18deg
+//
+// A TArrayD(2) is returned with the the mjd of the sunrise in
+// TArray[0] and the mjd of the sunset in TArrayD[1].
+//
+TArrayD MAstro::GetSunRiseSet(Double_t mjd, Double_t longit, Double_t latit, Double_t alt)
+{
+    const Double_t timediff = MAstro::GetTimeFromNoonToAlt(mjd, latit*TMath::DegToRad(), alt);
+    const Double_t noon     = MAstro::GetNoonTime(mjd, longit);
+
+    const Double_t N = TMath::Floor(mjd)+noon/24.;
+
+    const Double_t rise = timediff<0 ? N-0.5 : N-timediff/24.;
+    const Double_t set  = timediff<0 ? N+0.5 : N+timediff/24.;
+
+    TArrayD rc(2);
+    rc[0] = rise;
+    rc[1] = set;
+    return rc;
+}
+
+// --------------------------------------------------------------------------
+//
 // Returns the distance in x,y between two polar-vectors (eg. Alt/Az, Ra/Dec)
 // projected on aplain in a distance dist. For Magic this this the distance

trunk/MagicSoft/Mars/mastro/MAstro.h

@@ -9 +9 @@
 #endif
 
+class TArrayD;
 class TVector3;
 class MTime;
@@ -74 +75 @@
     static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi);
 
+    static void     GetMoonRaDec(Double_t mjd, Double_t &ra, Double_t &dec);
+
     static Double_t GetMoonPhase(Double_t mjd);
     static Double_t GetMoonPeriod(Double_t mjd);
     static Int_t    GetMagicPeriod(Double_t mjd);
+
+    // Estimate some parameters around the sun
+    static Double_t GetSunRaDec(Double_t mjd, Double_t &ra, Double_t &dec);
+    static Double_t GetEquationOfTime(Double_t mjd);
+    static Double_t GetNoonTime(Double_t mjd, Double_t longit);
+    static Double_t GetTimeFromNoonToAlt(Double_t mjd, Double_t latit, Double_t alt);
+    static TArrayD  GetSunRiseSet(Double_t mjd, Double_t longit, Double_t latit, Double_t alt=0);
 
     // Get distance between v1 and v0 in a plain perpendicular to v0 in distance dist

trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc

@@ -141 +141 @@
 #include <TPaveText.h>  // TPaveText
 
+#include <TGraph.h>
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -312 +314 @@
             continue;
 
-        if (epoch!="2000")
+        if (epoch.Atoi()!=2000)
         {
             gLog << warn << "Epoch != 2000... skipped." << endl;
@@ -646 +648 @@
     fList.AddLast(star);
     return 1;
+}
+
+// --------------------------------------------------------------------------
+//
+// Get the visibility curve (altitude vs time) for the current time
+// and observatory for the catalog entry with name name.
+// If name==0 the name of the TGraph is taken instead.
+// The day is divided into as many points as the graph has
+// points. If the graph has no points the default is 96.
+//
+void MAstroCatalog::GetVisibilityCurve(TGraph &g, const char *name) const
+{
+    if (!fTime || !fObservatory)
+    {
+        g.Set(0);
+        return;
+    }
+
+    MVector3 *star = static_cast<MVector3*>(FindObject(name ? name : g.GetName()));
+    if (!star)
+        return;
+
+    const Double_t mjd = TMath::Floor(fTime->GetMjd());
+    const Double_t lng = fObservatory->GetLongitudeDeg()/360;
+
+    if (g.GetN()==0)
+        g.Set(96);
+
+    for (int i=0; i<g.GetN(); i++)
+    {
+        const Double_t offset = (Double_t)i/g.GetN() - 0.5;
+
+        const MTime tm(mjd-lng+offset);
+
+        MVector3 v(*star);
+        v *= MAstroSky2Local(tm.GetGmst(), *fObservatory);
+
+        g.SetPoint(i, tm.GetAxisTime(), 90-v.Theta()*TMath::RadToDeg());
+    }
+
+    TH1   &h = *g.GetHistogram();
+    TAxis &x = *h.GetXaxis();
+    TAxis &y = *h.GetYaxis();
+
+    y.SetTitle("Altitude [\\circ]");
+    y.CenterTitle();
+
+    x.SetTitle("UTC");
+    x.CenterTitle();
+    x.SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
+    x.SetTimeDisplay(1);
+    x.SetLabelSize(0.033);
+
+    const Double_t atm = MTime(mjd).GetAxisTime();
+
+    x.SetRangeUser(atm-(0.5+lng)*24*60*60+15*60, atm+(0.5-lng)*24*60*60-15*60);
+
+    g.SetMinimum(5);
+    g.SetMaximum(90);
 }
 

trunk/MagicSoft/Mars/mastro/MAstroCatalog.h

@@ -22 +22 @@
 class TArrayI;
 class TGToolTip;
+class TGraph;
 
 class MAttLine : public TObject, public TAttLine
@@ -69 +70 @@
 protected:
     enum {
+        kMark         = BIT(14), // A mark for the sources in the list
         kHasChanged   = BIT(15), // Display has changed
         kGuiActive    = BIT(16), // GUI is interactive
@@ -152 +154 @@
     TList *GetList() { return &fList; } // Return list of stars
     UInt_t GetNumStars() const { return fList.GetEntries(); }
+    TObject *FindObject(const char *name) const { return fList.FindObject(name); }
+    TObject *FindObject(const TObject *obj) const { return fList.FindObject(obj); }
+    void MarkObject(const char *name) const { TObject *o=FindObject(name); if (o) o->SetBit(kMark); }
+
+    void GetVisibilityCurve(TGraph &g, const char *name=0) const;
 
     // TObject

trunk/MagicSoft/Mars/mastro/MObservatory.cc

@@ -34 +34 @@
 #include "MObservatory.h"
 
+#include <TArrayD.h>
 #include <TVector3.h>
 
@@ -92 +93 @@
         fObservatoryName = "Wuerzburg City";
         break;
+
+    case kTuorla:
+        fLatitude  = MAstro::Dms2Rad(60, 24, 57.0);
+        fLongitude = MAstro::Dms2Rad(22, 26, 42.0);
+        fHeight    = 53;
+        fObservatoryName = "Tuorla";
+        break;
     }
 
@@ -131 +139 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Get the time (as mjd) of sunrise/sunset at the day floor(mjd)
+// above/below alt[deg]
+//
+// For more information see MAstro::GetSunRiseSet
+//
+// A TArrayD(2) is returned with the sunrise in TArray[0] and the
+// sunset in TArrayD[1].
+//
+TArrayD MObservatory::GetSunRiseSet(Double_t mjd, Double_t alt) const
+{
+    return MAstro::GetSunRiseSet(mjd, GetLongitudeDeg(), GetLatitudeDeg(), alt);
+}
 
 // --------------------------------------------------------------------------

trunk/MagicSoft/Mars/mastro/MObservatory.h

@@ -7 +7 @@
 
 class MTime;
+class TArrayD;
 
 class MObservatory : public MParContainer
@@ -14 +15 @@
     {
         kMagic1,
-        kWuerzburgCity
+        kWuerzburgCity,
+        kTuorla
     };
 
@@ -67 +69 @@
     Double_t GetHeight() const          { return fHeight; }
 
+    TArrayD GetSunRiseSet(Double_t mjd, Double_t alt=0) const;
+
     void RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const;
     Double_t RotationAngle(Double_t theta, Double_t phi) const;