Changeset 8066
- Timestamp:
- 10/14/06 19:20:08 (18 years ago)
- Location:
- trunk/MagicSoft
- Files:
-
- 2 added
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r8065 r8066 59 59 * mmc/MMcRunHeader.cxx: 60 60 - initalize fImpactMax to -1 61 62 * macros/starvisday.C, macros/starvisyear.C: 63 - added 61 64 62 65 -
trunk/MagicSoft/Mars/NEWS
r8056 r8066 2 2 3 3 *** <cvs> 4 5 - general: changed the building of the sequence files such, that 6 the source-/project name in case of wobble-sources is not the 7 source/project of the first data run anymore but the new generalized 8 wobble source/project 9 10 - general: Added some functions to MAstro to calculate sun and 11 moon properties 12 13 - general: Added two new macros (starvisday.C and starvisyear.C) 14 which allow to plot the visibility of several stars in one day 15 or of one star along one year 4 16 5 17 -
trunk/MagicSoft/Mars/mastro/AstroIncl.h
r3326 r8066 2 2 3 3 #include <TArrayC.h> 4 #include <TArrayD.h> 4 5 #include "MArrayB.h" 5 6 #include "MArrayS.h" -
trunk/MagicSoft/Mars/mastro/MAstro.cc
r7550 r8066 33 33 #include <iostream> 34 34 35 #include <TArrayD.h> // TArrayD 35 36 #include <TVector3.h> // TVector3 36 37 … … 585 586 // -------------------------------------------------------------------------- 586 587 // 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 // 595 Double_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 // 647 void 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 // 697 Double_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 // 719 Double_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 // 746 Double_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 // 773 TArrayD 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 // 587 791 // Returns the distance in x,y between two polar-vectors (eg. Alt/Az, Ra/Dec) 588 792 // projected on aplain in a distance dist. For Magic this this the distance -
trunk/MagicSoft/Mars/mastro/MAstro.h
r7780 r8066 9 9 #endif 10 10 11 class TArrayD; 11 12 class TVector3; 12 13 class MTime; … … 74 75 static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi); 75 76 77 static void GetMoonRaDec(Double_t mjd, Double_t &ra, Double_t &dec); 78 76 79 static Double_t GetMoonPhase(Double_t mjd); 77 80 static Double_t GetMoonPeriod(Double_t mjd); 78 81 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); 79 89 80 90 // Get distance between v1 and v0 in a plain perpendicular to v0 in distance dist -
trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc
r7887 r8066 141 141 #include <TPaveText.h> // TPaveText 142 142 143 #include <TGraph.h> 144 143 145 #include "MLog.h" 144 146 #include "MLogManip.h" … … 312 314 continue; 313 315 314 if (epoch !="2000")316 if (epoch.Atoi()!=2000) 315 317 { 316 318 gLog << warn << "Epoch != 2000... skipped." << endl; … … 646 648 fList.AddLast(star); 647 649 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 // 660 void 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); 648 709 } 649 710 -
trunk/MagicSoft/Mars/mastro/MAstroCatalog.h
r7789 r8066 22 22 class TArrayI; 23 23 class TGToolTip; 24 class TGraph; 24 25 25 26 class MAttLine : public TObject, public TAttLine … … 69 70 protected: 70 71 enum { 72 kMark = BIT(14), // A mark for the sources in the list 71 73 kHasChanged = BIT(15), // Display has changed 72 74 kGuiActive = BIT(16), // GUI is interactive … … 152 154 TList *GetList() { return &fList; } // Return list of stars 153 155 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; 154 161 155 162 // TObject -
trunk/MagicSoft/Mars/mastro/MObservatory.cc
r5932 r8066 34 34 #include "MObservatory.h" 35 35 36 #include <TArrayD.h> 36 37 #include <TVector3.h> 37 38 … … 92 93 fObservatoryName = "Wuerzburg City"; 93 94 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; 94 102 } 95 103 … … 131 139 } 132 140 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 // 151 TArrayD MObservatory::GetSunRiseSet(Double_t mjd, Double_t alt) const 152 { 153 return MAstro::GetSunRiseSet(mjd, GetLongitudeDeg(), GetLatitudeDeg(), alt); 154 } 133 155 134 156 // -------------------------------------------------------------------------- -
trunk/MagicSoft/Mars/mastro/MObservatory.h
r5932 r8066 7 7 8 8 class MTime; 9 class TArrayD; 9 10 10 11 class MObservatory : public MParContainer … … 14 15 { 15 16 kMagic1, 16 kWuerzburgCity 17 kWuerzburgCity, 18 kTuorla 17 19 }; 18 20 … … 67 69 Double_t GetHeight() const { return fHeight; } 68 70 71 TArrayD GetSunRiseSet(Double_t mjd, Double_t alt=0) const; 72 69 73 void RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const; 70 74 Double_t RotationAngle(Double_t theta, Double_t phi) const; -
trunk/MagicSoft/Mars/mbase/MTime.h
r7953 r8066 97 97 98 98 // Getter functions 99 Double_t GetJD() const { return GetMjd()+2400000.5; } 99 100 Double_t GetMjd() const; 100 101 Double_t GetGmst() const; -
trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc
r8049 r8066 19 19 ! Author(s): Harald Kornmayer 1/2001 20 20 ! 21 ! Copyright: MAGIC Software Development, 2000-200 221 ! Copyright: MAGIC Software Development, 2000-2006 22 22 ! 23 23 ! -
trunk/MagicSoft/Mars/mhist/MHCamera.h
r7867 r8066 105 105 106 106 // This is a trick to remove TH1 entries from the context menu 107 TH1 *Rebin(Int_t ngroup=2, const char*newname="") { return this; } 107 TH1 *Rebin(Int_t ngroup=2, const char *newname="", const Double_t *bin=0) { return this; } 108 TH1 *Rebin(Int_t ngroup=2, const char *newname="") { return this; } 108 109 void DrawPanel() {} 109 110 -
trunk/MagicSoft/include-Classes/MMcFormat/MMcRunHeader.cxx
r7880 r8066 126 126 fShowerPhiMax = 0.0; 127 127 fShowerPhiMin = 0.0; 128 129 fImpactMax = -1; 128 130 129 131 fCWaveLower = 0.0;
Note:
See TracChangeset
for help on using the changeset viewer.