Ignore:
Timestamp:
10/14/06 19:20:08 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mastro
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mastro/AstroIncl.h

    r3326 r8066  
    22
    33#include <TArrayC.h>
     4#include <TArrayD.h>
    45#include "MArrayB.h"
    56#include "MArrayS.h"
  • trunk/MagicSoft/Mars/mastro/MAstro.cc

    r7550 r8066  
    3333#include <iostream>
    3434
     35#include <TArrayD.h>  // TArrayD
    3536#include <TVector3.h> // TVector3
    3637
     
    585586// --------------------------------------------------------------------------
    586587//
     588// Returns right ascension and declination [rad] of the sun at the
     589// given mjd (ra, dec).
     590//
     591// returns the mean longitude [rad].
     592//
     593// from http://xoomer.alice.it/vtomezzo/sunriset/formulas/index.html
     594//
     595Double_t MAstro::GetSunRaDec(Double_t mjd, Double_t &ra, Double_t &dec)
     596{
     597    const Double_t T = (mjd-51544.5)/36525;// +  (h-12)/24.0;
     598
     599    const Double_t T2 = T<0 ? -T*T : T*T;
     600    const Double_t T3 = T*T*T;
     601
     602    // Find the ecliptic longitude of the Sun
     603
     604    // Geometric mean longitude of the Sun
     605    const Double_t L = 280.46646 + 36000.76983*T + 0.0003032*T2;
     606
     607    // mean anomaly of the Sun
     608    Double_t g = 357.52911 + 35999.05029*T - 0.0001537*T2;
     609    g *= TMath::DegToRad();
     610
     611    // Longitude of the moon's ascending node
     612    Double_t omega = 125.04452 - 1934.136261*T + 0.0020708*T2 + T3/450000;
     613    omega *= TMath::DegToRad();
     614
     615    const Double_t coso = cos(omega);
     616    const Double_t sino = sin(omega);
     617
     618    // Equation of the center
     619    const Double_t C = (1.914602 - 0.004817*T - 0.000014*T2)*sin(g) +
     620        (0.019993 - 0.000101*T)*sin(2*g) + 0.000289*sin(3*g);
     621
     622    // True longitude of the sun
     623    const Double_t tlong = L + C;
     624
     625    // Apperent longitude of the Sun (ecliptic)
     626    Double_t lambda = tlong - 0.00569 - 0.00478*sino;
     627    lambda *= TMath::DegToRad();
     628
     629    // Obliquity of the ecliptic
     630    Double_t obliq = 23.4392911 - 0.01300416667*T - 0.00000016389*T2 + 0.00000050361*T3 + 0.00255625*coso;
     631    obliq *= TMath::DegToRad();
     632
     633    // Find the RA and DEC of the Sun
     634    const Double_t sinl =  sin(lambda);
     635
     636    ra  = atan2(cos(obliq) * sinl, cos(lambda));
     637    dec = asin(sin(obliq) * sinl);
     638
     639    return L*TMath::DegToRad();
     640}
     641
     642// --------------------------------------------------------------------------
     643//
     644// Returns right ascension and declination [rad] of the moon at the
     645// given mjd (ra, dec).
     646//
     647void MAstro::GetMoonRaDec(Double_t mjd, Double_t &ra, Double_t &dec)
     648{
     649    // Mean Moon orbit elements as of 1990.0
     650    const Double_t l0 = 318.351648 * TMath::DegToRad();
     651    const Double_t P0 =  36.340410 * TMath::DegToRad();
     652    const Double_t N0 = 318.510107 * TMath::DegToRad();
     653    const Double_t i  =   5.145396 * TMath::DegToRad();
     654
     655    Double_t sunra, sundec, g;
     656    {
     657        const Double_t T = (mjd-51544.5)/36525;// +  (h-12)/24.0;
     658        const Double_t T2 = T<0 ? -T*T : T*T;
     659
     660        GetSunRaDec(mjd, sunra, sundec);
     661
     662        // mean anomaly of the Sun
     663        g = 357.52911 + 35999.05029*T - 0.0001537*T2;
     664        g *= TMath::DegToRad();
     665    }
     666
     667    const Double_t sing   = sin(g)*TMath::DegToRad();
     668
     669    const Double_t D      = (mjd-47891) * TMath::DegToRad();
     670    const Double_t l      =    13.1763966*D + l0;
     671    const Double_t MMoon  = l  -0.1114041*D - P0; // Moon's mean anomaly M
     672    const Double_t N      = N0 -0.0529539*D;      // Moon's mean ascending node longitude
     673
     674    const Double_t C      = l-sunra;
     675    const Double_t Ev     = 1.2739 * sin(2*C-MMoon) * TMath::DegToRad();
     676    const Double_t Ae     = 0.1858 * sing;
     677    const Double_t A3     = 0.37   * sing;
     678    const Double_t MMoon2 = MMoon+Ev-Ae-A3;  // corrected Moon anomaly
     679
     680    const Double_t Ec     = 6.2886 * sin(MMoon2)  * TMath::DegToRad();  // equation of centre
     681    const Double_t A4     = 0.214  * sin(2*MMoon2)* TMath::DegToRad();
     682    const Double_t l2     = l+Ev+Ec-Ae+A4; // corrected Moon's longitude
     683
     684    const Double_t V      = 0.6583 * sin(2*(l2-sunra)) * TMath::DegToRad();
     685    const Double_t l3     = l2+V; // true orbital longitude;
     686
     687    const Double_t N2     = N -0.16*sing;
     688
     689    ra  = fmod( N2 + atan2( sin(l3-N2)*cos(i), cos(l3-N2) ), TMath::TwoPi() );
     690    dec = asin(sin(l3-N2)*sin(i) );
     691}
     692
     693// --------------------------------------------------------------------------
     694//
     695// Return Euqation of time in hours for given mjd
     696//
     697Double_t MAstro::GetEquationOfTime(Double_t mjd)
     698{
     699    Double_t ra, dec;
     700    const Double_t L = fmod(GetSunRaDec(mjd, ra, dec), TMath::TwoPi());
     701
     702    if (L-ra>TMath::Pi())
     703        ra += TMath::TwoPi();
     704
     705    return 24*(L - ra)/TMath::TwoPi();
     706}
     707
     708// --------------------------------------------------------------------------
     709//
     710// Returns noon time (the time of the highest altitude of the sun)
     711// at the given mjd and at the given observers longitude [deg]
     712//
     713// The maximum altitude reached at noon time is
     714//   altmax = 90.0 + dec - latit;
     715//   if (dec > latit)
     716//      altmax = 90.0 + latit - dec;
     717// dec=Declination of the sun
     718//
     719Double_t MAstro::GetNoonTime(Double_t mjd, Double_t longit)
     720{
     721    const Double_t equation = GetEquationOfTime(TMath::Floor(mjd));
     722    return 12. + equation - longit/15;
     723}
     724
     725// --------------------------------------------------------------------------
     726//
     727// Returns the time (in hours) between noon (the sun culmination)
     728// and the sun being at height alt[deg] (90=zenith, 0=horizont)
     729//
     730//       civil twilight:      0deg to  -6deg
     731//       nautical twilight:  -6deg to -12deg
     732//       astronom twilight: -12deg to -18deg
     733//
     734// latit is the observers latitude in rad
     735//
     736// returns -1 in case the sun doesn't reach this altitude.
     737// (eg. alt=0: Polarnight or -day)
     738//
     739// To get the sun rise/set:
     740//    double timediff = MAstro::GetTimeFromNoonToAlt(mjd, latit*TMath::DegToRad(), par[0]);
     741//    double noon     = MAstro::GetNoonTime(mjd, longit);
     742//    double N        = TMath::Floor(mjd)+noon/24.;
     743//    double risetime = N-timediff/24.;
     744//    double settime  = N+timediff/24.;
     745//
     746Double_t MAstro::GetTimeFromNoonToAlt(Double_t mjd, Double_t latit, Double_t alt)
     747{
     748    Double_t ra, dec;
     749    GetSunRaDec(mjd, ra, dec);
     750
     751    const Double_t h = alt*TMath::DegToRad();
     752
     753    const Double_t arg = (sin(h) - sin(latit)*sin(dec))/(cos(latit)*cos(dec));
     754
     755    return TMath::Abs(arg)>1 ? -1 : 12*acos(arg)/TMath::Pi();
     756}
     757
     758// --------------------------------------------------------------------------
     759//
     760// Returns the time of the sunrise/set calculated before and after
     761// the noon of floor(mjd) (TO BE IMPROVED)
     762//
     763// Being longit and latit the longitude and latitude of the observer
     764// in deg and alt the hight above or below the horizont in deg.
     765//
     766//       civil twilight:      0deg to  -6deg
     767//       nautical twilight:  -6deg to -12deg
     768//       astronom twilight: -12deg to -18deg
     769//
     770// A TArrayD(2) is returned with the the mjd of the sunrise in
     771// TArray[0] and the mjd of the sunset in TArrayD[1].
     772//
     773TArrayD MAstro::GetSunRiseSet(Double_t mjd, Double_t longit, Double_t latit, Double_t alt)
     774{
     775    const Double_t timediff = MAstro::GetTimeFromNoonToAlt(mjd, latit*TMath::DegToRad(), alt);
     776    const Double_t noon     = MAstro::GetNoonTime(mjd, longit);
     777
     778    const Double_t N = TMath::Floor(mjd)+noon/24.;
     779
     780    const Double_t rise = timediff<0 ? N-0.5 : N-timediff/24.;
     781    const Double_t set  = timediff<0 ? N+0.5 : N+timediff/24.;
     782
     783    TArrayD rc(2);
     784    rc[0] = rise;
     785    rc[1] = set;
     786    return rc;
     787}
     788
     789// --------------------------------------------------------------------------
     790//
    587791// Returns the distance in x,y between two polar-vectors (eg. Alt/Az, Ra/Dec)
    588792// projected on aplain in a distance dist. For Magic this this the distance
  • trunk/MagicSoft/Mars/mastro/MAstro.h

    r7780 r8066  
    99#endif
    1010
     11class TArrayD;
    1112class TVector3;
    1213class MTime;
     
    7475    static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi);
    7576
     77    static void     GetMoonRaDec(Double_t mjd, Double_t &ra, Double_t &dec);
     78
    7679    static Double_t GetMoonPhase(Double_t mjd);
    7780    static Double_t GetMoonPeriod(Double_t mjd);
    7881    static Int_t    GetMagicPeriod(Double_t mjd);
     82
     83    // Estimate some parameters around the sun
     84    static Double_t GetSunRaDec(Double_t mjd, Double_t &ra, Double_t &dec);
     85    static Double_t GetEquationOfTime(Double_t mjd);
     86    static Double_t GetNoonTime(Double_t mjd, Double_t longit);
     87    static Double_t GetTimeFromNoonToAlt(Double_t mjd, Double_t latit, Double_t alt);
     88    static TArrayD  GetSunRiseSet(Double_t mjd, Double_t longit, Double_t latit, Double_t alt=0);
    7989
    8090    // Get distance between v1 and v0 in a plain perpendicular to v0 in distance dist
  • trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc

    r7887 r8066  
    141141#include <TPaveText.h>  // TPaveText
    142142
     143#include <TGraph.h>
     144
    143145#include "MLog.h"
    144146#include "MLogManip.h"
     
    312314            continue;
    313315
    314         if (epoch!="2000")
     316        if (epoch.Atoi()!=2000)
    315317        {
    316318            gLog << warn << "Epoch != 2000... skipped." << endl;
     
    646648    fList.AddLast(star);
    647649    return 1;
     650}
     651
     652// --------------------------------------------------------------------------
     653//
     654// Get the visibility curve (altitude vs time) for the current time
     655// and observatory for the catalog entry with name name.
     656// If name==0 the name of the TGraph is taken instead.
     657// The day is divided into as many points as the graph has
     658// points. If the graph has no points the default is 96.
     659//
     660void MAstroCatalog::GetVisibilityCurve(TGraph &g, const char *name) const
     661{
     662    if (!fTime || !fObservatory)
     663    {
     664        g.Set(0);
     665        return;
     666    }
     667
     668    MVector3 *star = static_cast<MVector3*>(FindObject(name ? name : g.GetName()));
     669    if (!star)
     670        return;
     671
     672    const Double_t mjd = TMath::Floor(fTime->GetMjd());
     673    const Double_t lng = fObservatory->GetLongitudeDeg()/360;
     674
     675    if (g.GetN()==0)
     676        g.Set(96);
     677
     678    for (int i=0; i<g.GetN(); i++)
     679    {
     680        const Double_t offset = (Double_t)i/g.GetN() - 0.5;
     681
     682        const MTime tm(mjd-lng+offset);
     683
     684        MVector3 v(*star);
     685        v *= MAstroSky2Local(tm.GetGmst(), *fObservatory);
     686
     687        g.SetPoint(i, tm.GetAxisTime(), 90-v.Theta()*TMath::RadToDeg());
     688    }
     689
     690    TH1   &h = *g.GetHistogram();
     691    TAxis &x = *h.GetXaxis();
     692    TAxis &y = *h.GetYaxis();
     693
     694    y.SetTitle("Altitude [\\circ]");
     695    y.CenterTitle();
     696
     697    x.SetTitle("UTC");
     698    x.CenterTitle();
     699    x.SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
     700    x.SetTimeDisplay(1);
     701    x.SetLabelSize(0.033);
     702
     703    const Double_t atm = MTime(mjd).GetAxisTime();
     704
     705    x.SetRangeUser(atm-(0.5+lng)*24*60*60+15*60, atm+(0.5-lng)*24*60*60-15*60);
     706
     707    g.SetMinimum(5);
     708    g.SetMaximum(90);
    648709}
    649710
  • trunk/MagicSoft/Mars/mastro/MAstroCatalog.h

    r7789 r8066  
    2222class TArrayI;
    2323class TGToolTip;
     24class TGraph;
    2425
    2526class MAttLine : public TObject, public TAttLine
     
    6970protected:
    7071    enum {
     72        kMark         = BIT(14), // A mark for the sources in the list
    7173        kHasChanged   = BIT(15), // Display has changed
    7274        kGuiActive    = BIT(16), // GUI is interactive
     
    152154    TList *GetList() { return &fList; } // Return list of stars
    153155    UInt_t GetNumStars() const { return fList.GetEntries(); }
     156    TObject *FindObject(const char *name) const { return fList.FindObject(name); }
     157    TObject *FindObject(const TObject *obj) const { return fList.FindObject(obj); }
     158    void MarkObject(const char *name) const { TObject *o=FindObject(name); if (o) o->SetBit(kMark); }
     159
     160    void GetVisibilityCurve(TGraph &g, const char *name=0) const;
    154161
    155162    // TObject
  • trunk/MagicSoft/Mars/mastro/MObservatory.cc

    r5932 r8066  
    3434#include "MObservatory.h"
    3535
     36#include <TArrayD.h>
    3637#include <TVector3.h>
    3738
     
    9293        fObservatoryName = "Wuerzburg City";
    9394        break;
     95
     96    case kTuorla:
     97        fLatitude  = MAstro::Dms2Rad(60, 24, 57.0);
     98        fLongitude = MAstro::Dms2Rad(22, 26, 42.0);
     99        fHeight    = 53;
     100        fObservatoryName = "Tuorla";
     101        break;
    94102    }
    95103
     
    131139}
    132140
     141// --------------------------------------------------------------------------
     142//
     143// Get the time (as mjd) of sunrise/sunset at the day floor(mjd)
     144// above/below alt[deg]
     145//
     146// For more information see MAstro::GetSunRiseSet
     147//
     148// A TArrayD(2) is returned with the sunrise in TArray[0] and the
     149// sunset in TArrayD[1].
     150//
     151TArrayD MObservatory::GetSunRiseSet(Double_t mjd, Double_t alt) const
     152{
     153    return MAstro::GetSunRiseSet(mjd, GetLongitudeDeg(), GetLatitudeDeg(), alt);
     154}
    133155
    134156// --------------------------------------------------------------------------
  • trunk/MagicSoft/Mars/mastro/MObservatory.h

    r5932 r8066  
    77
    88class MTime;
     9class TArrayD;
    910
    1011class MObservatory : public MParContainer
     
    1415    {
    1516        kMagic1,
    16         kWuerzburgCity
     17        kWuerzburgCity,
     18        kTuorla
    1719    };
    1820
     
    6769    Double_t GetHeight() const          { return fHeight; }
    6870
     71    TArrayD GetSunRiseSet(Double_t mjd, Double_t alt=0) const;
     72
    6973    void RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const;
    7074    Double_t RotationAngle(Double_t theta, Double_t phi) const;
Note: See TracChangeset for help on using the changeset viewer.