source: branches/Mars_use_drstimefiles/macros/starvisyear.C@ 20051

Last change on this file since 20051 was 14152, checked in by tbretz, 12 years ago
Added output of the observatory; fixed a typo in a comment.
File size: 6.5 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 10/2006 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2004-2006
21!
22!
23\* ======================================================================== */
24
25///////////////////////////////////////////////////////////////////////////
26//
27// starvisyear.C
28//
29// displays the star visibility along one year
30//
31// See the code for ducumentation and setup
32//
33///////////////////////////////////////////////////////////////////////////
34
35#include <TString.h>
36#include <TCanvas.h>
37#include <TH2D.h>
38#include <TGraph.h>
39
40#include "MObservatory.h"
41#include "MTime.h"
42#include "MAstro.h"
43#include "MAstroSky2Local.h"
44#include "MVector3.h"
45#include "MH.h"
46
47
48void starvisyear()
49{
50 // Setup the observatory location (see MObservatory for details)
51 MObservatory obs(MObservatory::kMagic1);;
52
53 // Setup the start time of the year (365 days) to be displayed
54 MTime time(-1); // -1 mean NOW
55
56 // Setup a different time if you like
57 time.Set(2006, 1, 1, 0);
58
59 // Setup right ascension and declination of the source to display
60 //Crab
61 const Double_t ra = MAstro::Hms2Rad(5, 34, 31.9);
62 const Double_t dec = MAstro::Dms2Rad(22, 0, 52.0);
63
64 cout << obs.Print() << endl;
65
66 return;
67
68 // Setup the name of the source
69 TString source("CrabNebula");
70
71 // Setup palette you would like to see (99 is fixed)
72 // for more details see MH::SetPalette
73 MH::SetPalette("glow2", 99);
74
75 // --------------------------------------------------------------------------
76 //
77 // Start producing the data for the plot
78 //
79
80 // Setup the source in the appropriate format
81 MVector3 v;
82 v.SetRaDec(ra, dec);
83
84 // For simplicity get mjd of time
85 Double_t mjd = time.GetMjd();
86
87 // Setup axis ranges
88 Double_t first = time.GetAxisTime();
89 Double_t last = MTime(mjd+365).GetAxisTime();
90
91 Int_t h0 = 13-obs.GetLongitudeDeg()/15;
92
93 // Define histograms
94 TH2D hist( "1", "", 365, first, last, 24*12, first+h0*60*60, first+(h0+24)*60*60);//h0/24.+1, (h0+24)/24.+1);
95 TH2D hmoon("2", "", 365, first, last, 24*12, first+h0*60*60, first+(h0+24)*60*60);//h0/24.+1, (h0+24)/24.+1);
96
97 // Fill histograms
98 TGraph sunset[4], sunrise[4];
99 for (int x=0; x<365; x++)
100 {
101 // Get moon phase and axis time for mjd+x
102 Double_t moon = MAstro::GetMoonPhase(mjd+x);
103 Double_t tmx = MTime(mjd+x).GetAxisTime();
104
105 // get sunrise and sunset for the current day
106 Double_t sunst[4], sunri[4];
107 for (int i=0; i<4; i++)
108 {
109 // get sunset and sunrise of next day
110 sunri[i] = obs.GetSunRiseSet(mjd+x+1,-18+6*i)[0];
111 sunst[i] = obs.GetSunRiseSet(mjd+x, -18+6*i)[1];
112
113 // fill both into the TGraphs
114 sunset[i].SetPoint( sunset[i].GetN(), tmx, MTime(sunst[i]).GetAxisTime()-x*24*60*60);
115 sunrise[i].SetPoint(sunrise[i].GetN(), tmx, MTime(sunri[i]).GetAxisTime()-x*24*60*60);
116 }
117
118 // If nautical twilight doesn't take place skip the rest
119 if (sunri[1]<0 || sunst[1]<0)
120 continue;
121
122 // Now calculate in steps of five minutes
123 for (int hor=0; hor<24*12; hor++)
124 {
125 Double_t hr = hor/12. + h0;
126 Double_t offset = x + hr/24.;
127
128 // mjd of the current time
129 MTime tm(mjd+offset);
130
131 // Check if current time is within darkness
132 if (tm.GetMjd()>sunri[1] || tm.GetMjd()<sunst[1])
133 continue;
134
135 // Initialize conversion from sky coordinates to local coordinates
136 MAstroSky2Local conv(tm.GetGmst(), obs);
137
138 // get moon position on the sky
139 Double_t moonra, moondec;
140 MAstro::GetMoonRaDec(tm.GetMjd(), moonra, moondec);
141
142 // calculate moon position as seen by the telescope
143 v.SetRaDec(moonra, moondec);
144 v *= conv;
145
146 // check if moon is above or below horizont
147 if (v.Theta()*TMath::RadToDeg()<90)
148 moon = 0;
149
150 // fill moon visibility into histpogram
151 hmoon.SetBinContent(x+1, hor+1, moon);
152
153 // calculate source position as seen by the telescope
154 v.SetRaDec(ra, dec);
155 v *= conv;
156
157 // fill altitude of source into histogram
158 hist.SetBinContent(x+1, hor+1, v.Theta()*TMath::RadToDeg());
159 }
160 }
161
162 // --------------------------------------------------------------------------
163 //
164 // Start producing the nice plot
165 //
166
167 // Setup canvas
168 TCanvas *c = new TCanvas;
169 c->SetBorderMode(0);
170 c->SetGridx();
171 c->SetGridy();
172
173 // Setup histogram
174 TString name = obs.GetObservatoryName();
175 hist.SetTitle(Form("Zenith distance of %s at %s", source.Data(), name.Data()));
176
177 hist.SetXTitle("Day");
178 hist.SetYTitle("UTC");
179
180 hist.GetXaxis()->CenterTitle();
181 hist.GetYaxis()->CenterTitle();
182
183 hist.GetXaxis()->SetTimeFormat("%d/%m/%y %F1995-01-01 00:00:00 GMT");
184 hist.GetXaxis()->SetTimeDisplay(1);
185 hist.GetXaxis()->SetLabelSize(0.033);
186
187 hist.GetYaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
188 hist.GetYaxis()->SetTimeDisplay(1);
189 hist.GetYaxis()->SetLabelSize(0.033);
190 hist.GetYaxis()->SetTitleOffset(1.3);
191
192 hist.SetContour(99);
193 hist.SetMinimum(0);
194 hist.SetMaximum(90);
195
196 hist.SetStats(kFALSE);
197
198 // draw histograms
199 hist.DrawCopy("colz");
200 hmoon.DrawCopy("scat=7 same");
201
202 // Draw rise and set time of sun
203 Int_t style[] = { kSolid, 5, 4, kDotted };
204 for (int i=0; i<4; i++)
205 {
206 sunset[i].SetLineColor(kBlue);
207 sunrise[i].SetLineColor(kBlue);
208 sunset[i].SetLineWidth(2);
209 sunrise[i].SetLineWidth(2);
210 sunset[i].SetLineStyle(style[i]);
211 sunrise[i].SetLineStyle(style[i]);
212 sunset[i].DrawClone("L");
213 sunrise[i].DrawClone("L");
214 }
215}
Note: See TracBrowser for help on using the repository browser.