source: trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc@ 1681

Last change on this file since 1681 was 1681, checked in by rwagner, 22 years ago
Preliminary version of classes for the Visibility Plotter. You need slalib installed in Mars/.. Makefile therefore is currently not included in central MARS Makefile.
File size: 16.7 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): Robert Wagner <magicdev@rwagner.de> 12/2002
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26// //
27// MVPPlotter //
28// //
29/////////////////////////////////////////////////////////////////////////////
30#include "MVPPlotter.h"
31
32#include <stream.h>
33
34#include "TCanvas.h"
35#include "TH1.h"
36#include "TH2.h"
37#include "TPaveText.h"
38#include "TPaveLabel.h"
39#include "TGraph.h"
40#include "TString.h"
41#include "TStyle.h"
42
43#include "MParList.h"
44#include "MVPTime.h"
45
46ClassImp(MVPPlotter);
47
48// --------------------------------------------------------------------------
49//
50// Default constructor.
51//
52MVPPlotter::MVPPlotter(const char *name, const char *title) : fUseSun(kTRUE), fUseMoon(kTRUE), fUsePlanets(kTRUE), fAstronomicalDarkness(-18.0)
53{
54 // fName = name ? name : "MVPPlotter";
55 // fTitle = title ? title : "Generates visibility histograms and information";
56
57 fgSecPerDay = 86400;
58 fgMJD010170 = 40586; // 01-01-70 is JD 40586
59 fgDegToRad = 2*TMath::Pi()/360;
60
61}
62
63MVPPlotter::~MVPPlotter()
64// --------------------------------------------------------------------------
65//
66// Destructor. Deletes objects if allocated.
67//
68{
69 if (fUseSun) delete fSun;
70 if (fUseMoon) delete fMoon;
71 if (fUsePlanets) {
72 delete fVenus;
73 delete fMars;
74 delete fJupiter;
75 delete fSaturn;
76 }
77}
78
79void MVPPlotter::SetupObjects()
80// --------------------------------------------------------------------------
81//
82// Create instances for Sun, Moon, Planets if requested.
83//
84{
85
86 fTime = new MVPTime();
87 fPlist.AddToList(fTime);
88 fPlist.AddToList(fObs);
89
90 if (fObject==NULL) {
91 cout << "You didn't specify an object!" << endl;
92 }
93
94 fObject->PreProcess(&fPlist);
95 cout << "Processing object " << fObject->GetObjectName() << endl;
96
97 if (fUseSun) {
98 fSun = new MVPObject();
99 fSun->SetObject(0); // Sun
100 fSun->PreProcess(&fPlist);
101 fSun->SetCalcEc(kTRUE);
102 }
103
104 if (fUseMoon) {
105 fMoon = new MVPObject();
106 fMoon->SetObject(3); // Moon
107 fMoon->PreProcess(&fPlist);
108 fMoon->SetCalcEc(kTRUE);
109 }
110
111 if (fUsePlanets) {
112 fVenus = new MVPObject();
113 fVenus->SetObject(2);
114 fVenus->PreProcess(&fPlist);
115 fVenus->SetCalcEc(kTRUE);
116
117 fMars = new MVPObject();
118 fMars->SetObject(4);
119 fMars->PreProcess(&fPlist);
120 fMars->SetCalcEc(kTRUE);
121
122 fJupiter = new MVPObject();
123 fJupiter->SetObject(5);
124 fJupiter->PreProcess(&fPlist);
125 fJupiter->SetCalcEc(kTRUE);
126
127 fSaturn = new MVPObject();
128 fSaturn->SetObject(6);
129 fSaturn->PreProcess(&fPlist);
130 fSaturn->SetCalcEc(kTRUE);
131 }
132}
133
134// --------------------------------------------------------------------------
135//
136// Plots for a single object and a whole year are generated.
137//
138// Currently we do the following:
139// - Create a plot MJD vs UTC for one year
140// - Create a plot maxObjHeight vs Week for one year
141// - Create visibility tables for one year (ZA, hours)
142//
143Bool_t MVPPlotter::CalcYear(Int_t year, UInt_t daySlices=96)
144{
145 SetupObjects();
146
147 UInt_t startday = (UInt_t)fTime->MJDStartOfYear(year);
148 UInt_t stopday = (UInt_t)fTime->MJDStartOfYear(year+1)-1;
149
150 cout << "Processing period MJD "<< startday << " to MJD "<< stopday << endl;
151
152 ULong_t startdayROOT = fgSecPerDay * (startday-fgMJD010170);
153 ULong_t stopdayROOT = fgSecPerDay * (stopday-fgMJD010170);
154
155 // 2D Plot ZA vs MJD and UTC
156 fMjdUtcYear = new TH2D("fMjdUtcYear", "Visibility of object",
157 stopday-startday+1,startdayROOT,stopdayROOT,
158 daySlices+1,-450,fgSecPerDay+450);
159
160 // Observability hours per day MJD
161 fMjdObsHours = new TH1D("fMjdObsHours", "Observation hours per day",
162 stopday-startday+1,startdayROOT,stopdayROOT);
163
164 if (fUseMoon) {
165 // 2D Plot ZA of moon vs MJD and UTC
166 fMjdUtcYearMoon =new TH2D("fMjdUtcYearMoon", "Moon ZA",
167 stopday-startday+1,startdayROOT,stopdayROOT,
168 daySlices+1,-450,fgSecPerDay+450);
169
170 // Moon phase vs MJD
171 fMjdMoonPhase =new TH1D("fMjdMoonPhase", "Moon phase",
172 stopday-startday+1,startdayROOT,stopdayROOT);
173 // Moon distance of object vs MJD
174 fMjdMoonDistance =new TH1D("fMjdMoonDistance", "Moon distance",
175 stopday-startday+1,startdayROOT,stopdayROOT);
176 // Moon intensity at object vs MJD
177 fMjdMoonIntensity =new TH1D("fMjdMoonIntensity", "Moon intensity at locus of object",
178 stopday-startday+1,startdayROOT,stopdayROOT);
179
180 }
181
182 if (fUsePlanets) {
183 // Distance of closest planet vs MJD
184 fMjdPlanetDistance =new TH1D("fMjdPlanetDistance", "PlanetDistance",
185 stopday-startday+1,startdayROOT,stopdayROOT);
186 }
187
188
189 // Array which holds total visibility time for whole year [0] and
190 // each month [1]..[12]
191 Float_t visibility[13][18];
192 for (int m=0;m<13;m++)
193 for (int z=0;z<18;z++)
194 visibility[m][z]=0;
195
196 int fday, ftime;
197 Double_t phase=0;
198 Double_t obsHours;
199
200 for (UInt_t day=startday; day<stopday+1; day++) {
201 Double_t todaysPhase=0;
202 Double_t moonIntensity;
203 UInt_t i;
204 obsHours=0;
205 for (i=0; i<daySlices; i++)
206 {
207 // Rearrange filling of bins such that a "day" doesn't start at midnight,
208 // but rather at noon. This results in plots in which a "day" covers
209 // a whole night.
210 if (i>=(daySlices/2)) {
211 fday=day;
212 ftime=i-(daySlices/2);
213 } else {
214 fday=day-1;
215 ftime=i+(daySlices/2);
216 }
217
218 // Objects access fTime via parameter list...
219 fTime->SetMJD(day,(Double_t)i/daySlices);
220
221 if (fUseSun) fSun->Process();
222 if (fUseMoon) fMoon->Process();
223
224 if (fUseSun && fUseMoon) {
225
226 // Calculate moon phase...
227 phase = fSun->GetEcLong() - fMoon->GetEcLong();
228 phase = TMath::Pi()-(acos(cos(phase)*cos(fMoon->GetEcLat())));
229 phase = phase*180/TMath::Pi();
230 todaysPhase+=phase;
231
232 }
233
234 // If sun is not up (or we should not use sun information...)
235 if (fSun->GetAltitudeDeg() < fAstronomicalDarkness) {
236 // Calculate Position of object:
237 fObject->Process();
238
239 // Calculate moon brightness at locus of object
240 // now this is gonna be fun...
241
242 /* Evaluates predicted LUNAR part of sky brightness, in
243 * V magnitudes per square arcsecond,
244 */
245
246 moonIntensity = LunSkyBright(phase/180, fObject->GetDistance(fMoon),
247 fMoon->GetAltitudeDeg(), fObject->GetAltitudeDeg());
248 fMjdMoonIntensity->Fill(fgSecPerDay*(fday-fgMJD010170),moonIntensity);
249
250
251
252 // If moon is not up (or we should not use moon information...)
253 if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) {
254 // Fill MJD-UTC histogram
255 fMjdUtcYear->Fill(fgSecPerDay*(fday-fgMJD010170),fgSecPerDay*ftime/daySlices,fObject->GetAltitudeDeg());
256 }
257
258 // Sum up visibility time (only if moon not visible or moon
259 // info shouldn't be used at all)
260 if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) {
261 // Calculate
262 for (int z=0;z<18;z++) {
263 if (fObject->GetAltitudeDeg()>(z*5)) {
264 visibility[0][z]+=(Double_t)(60*24/daySlices);
265 visibility[(Int_t)fTime->GetMonth()][z]+=(Double_t)(60*24/daySlices);
266 }
267 } //for
268
269 if ((fObject->GetAltitudeDeg())>40) {
270 fMjdObsHours->Fill(fgSecPerDay*(fday-fgMJD010170),(Double_t)(24/(Double_t)daySlices));
271 }
272
273 }
274
275
276 } //fi sun
277
278 // Optional: Fill histo with moon-up times...
279 // if (fMoon->GetAltitudeDeg() >0)
280 // fMjdUtcYearMoon->Fill(fgSecPerDay*(day-fgMJD010170),fgSecPerDay*i/daySlices,phase);
281 // fMoon->Process();
282 // Double_t phase;
283 // phase = fSun->GetEcLong() - moon->GetEcLong();
284 // phase = TMath::Pi()-(acos(cos(phase)*cos(moon->GetEcLat())));
285 // cout << "Phase: " << phase*180/TMath::Pi() << endl;
286
287 } //for daySlices
288
289 // Distance of Venus to object
290
291 if (fUsePlanets) {
292 fVenus->Process();
293 fMars->Process();
294 fJupiter->Process();
295 fSaturn->Process();
296
297 Double_t distance = fVenus->GetDistance(fObject);
298 distance = distance < fMars->GetDistance(fObject) ? distance : fMars->GetDistance(fObject);
299 distance = distance < fJupiter->GetDistance(fObject) ? distance : fJupiter->GetDistance(fObject);
300 distance = distance < fSaturn->GetDistance(fObject) ? distance : fSaturn->GetDistance(fObject);
301
302 fMjdPlanetDistance->Fill(fgSecPerDay*(fday-fgMJD010170),distance);
303 }
304
305 fMjdMoonPhase->Fill(fgSecPerDay*(fday-fgMJD010170),todaysPhase/i);
306 fMjdMoonDistance->Fill(fgSecPerDay*(fday-fgMJD010170),fMoon->GetDistance(fObject));
307
308
309 } //for days
310
311
312 // Here we print the tables with visibilities...
313 cout << "Visibility time [hours]: " << endl;
314
315 for (int z=1;z<17;z++) {
316 if (visibility[0][z]==0) break;
317 printf("Alt>%2d|%6d|",z*5,(Int_t)(visibility[0][z]/60));
318 for (int m=1;m<13;m++) {
319 printf("%5d ",(Int_t)(visibility[m][z]/60));
320 }
321 cout<<endl;
322 }
323
324 int vistimestart=0;
325 int vistimeend=0;
326 for (int m=1;m<13;m++) {
327 int n = (m==1 ? 12 : m-1);
328 if (((visibility[m][8]/60)>20)&&((visibility[n][8])/60 <= 20)) {
329 vistimestart=m; // we're on the rising slope
330 }
331 }
332
333 for (int m=1;m<13;m++) {
334 int n = (m==1 ? 12 : m-1);
335 if (((visibility[m][8]/60)<20)&&((visibility[n][8]/60)>=20)) {
336 vistimeend=n; // we're on the falling slope
337 }
338 }
339
340 cout << "Visibility (Alt>40) during months: " << vistimestart << "-" << vistimeend << " approx " << visibility[0][8]/60 << "hrs" <<endl;
341
342
343 gROOT->Reset();
344
345 TCanvas *cMjdUtcYear = new TCanvas("cMjdUtcYear", "Object Visibility MjdUtcYear", 1100,500);
346 cMjdUtcYear->cd(0);
347
348 gStyle->SetPalette(1);
349 // gStyle->SetOptStat(0);
350 gStyle->SetOptFit(0);
351 gStyle->SetFillStyle(0);
352 gStyle->SetFillColor(10);
353 gStyle->SetCanvasColor(10);
354 gStyle->SetDrawBorder(0);
355 gStyle->SetPadColor(10);
356 gStyle->SetPadBorderSize(0);
357 gStyle->SetPadLeftMargin(0.12);
358 gStyle->SetTitleYOffset(1.2);
359 gStyle->SetTitleXOffset(1.2);
360
361 gStyle->SetPadLeftMargin(0.01);
362 gStyle->SetPadRightMargin(0.1);
363 gStyle->SetPadTopMargin(0.03);
364 gStyle->SetPadBottomMargin(0.12);
365
366 fMjdUtcYear->SetTitle(fObject->GetObjectName());
367 fMjdUtcYear->GetXaxis()->SetTimeDisplay(1);
368 fMjdUtcYear->GetXaxis()->SetTimeFormat("%b %y");
369 fMjdUtcYear->GetYaxis()->SetTimeDisplay(1);
370 fMjdUtcYear->GetYaxis()->SetTimeFormat("%Hh");
371 gStyle->SetTimeOffset(43200);
372 fMjdUtcYear->GetYaxis()->SetLabelOffset(0.01);
373 fMjdUtcYear->SetMinimum(0);
374 fMjdUtcYear->SetMaximum(90);
375
376 fMjdUtcYear->GetYaxis()->SetTitle("Hour");
377 fMjdUtcYear->GetXaxis()->SetTitle("Day of year");
378 fMjdUtcYear->GetZaxis()->SetTitle("Altitude");
379 fMjdUtcYear->Draw("SURF2BB9Z");
380 gPad->SetPhi(0);
381 gPad->SetTheta(90);
382 gPad->Modified();
383
384 TPaveText *pt = new TPaveText(-0.74,-0.35,-0.58,0.52,"");
385 char ptLine[80];
386 pt->AddText("Visibility time [hrs]:");
387 pt->SetTextAlign(13);
388 pt->SetFillStyle(0);
389 pt->SetBorderSize(0);
390 pt->SetTextFont(42);
391 for (int j=1;j<17;j++) {
392 sprintf (ptLine, "Alt>%i: %i", j*5, (Int_t)visibility[0][j]/60);
393 pt->AddText(ptLine);
394 if (visibility[0][j]==0) break;
395 }
396 pt->Draw();
397
398
399 if (fUseMoon) {
400 TCanvas *cMjdMoon = new TCanvas("cMjdMoon", "the Moon phase", 600,600);
401 gStyle->SetOptTitle(1);
402 cMjdMoon->Divide(1,3);
403 cMjdMoon->cd(1);
404 fMjdMoonPhase->Draw();
405 cMjdMoon->cd(2);
406 fMjdMoonDistance->Draw();
407 cMjdMoon->cd(3);
408 fMjdMoonIntensity->Draw();
409 }
410
411 TCanvas *cObsHours = new TCanvas("cObsHours", "ObsHours per day", 600,400);
412 fMjdObsHours->Draw();
413
414
415
416 if (fUsePlanets) {
417 TCanvas *cMjdPlanets = new TCanvas("cMjdPlanets", "the distance to planets", 600,200);
418 cMjdPlanets->cd(0);
419 fMjdPlanetDistance->Draw();
420 }
421
422
423 // next histogram
424 Float_t objectHeight[55];
425 Float_t simpleCounter[55];
426 Int_t weekCounter=0;
427
428 simpleCounter[weekCounter]=0;
429 objectHeight[weekCounter]=0;
430 weekCounter++;
431
432 Int_t startday2 = fTime->MJDStartOfYear(year);
433 Int_t stopday2 = fTime->MJDStartOfYear(year+1)-1;
434
435 for (int day=startday2; day<stopday2+1; day=day+7)
436 {
437 simpleCounter[weekCounter]=weekCounter-1;
438 objectHeight[weekCounter]=0;
439 for (int i=0; i<daySlices; i++)
440 {
441 if (i>=48) {
442 fday=day;
443 ftime=i-48;
444 } else {
445 fday=day-1;
446 ftime=i+48;
447 }
448
449 fTime->SetMJD(day,(Double_t)i/daySlices);
450 fObject->Process();
451 fSun->Process();
452
453 if (fSun->GetAltitudeDeg() < -25.0) {
454 if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg()))
455 objectHeight[weekCounter]=fObject->GetAltitudeDeg();
456 }
457
458 } //i
459 weekCounter++;
460 } //day
461 simpleCounter[weekCounter]=weekCounter-2;
462 objectHeight[weekCounter]=0;
463 weekCounter++;
464
465 TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
466
467 TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100);
468
469 // gStyle->SetOptTitle(0);
470// gStyle->SetPadLeftMargin(0.000001);
471// gStyle->SetPadRightMargin(0.000001);
472// gStyle->SetPadTopMargin(0.001);
473// gStyle->SetPadBottomMargin(0.000001);
474 gPad->SetGrid();
475
476 c2->SetGrid();
477 TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight);
478
479 tg->SetMinimum(1);
480 tg->SetMaximum(90);
481
482 Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude());
483
484 if (maxza > 90) maxza=90;
485 if (maxza < 0) maxza=00;
486
487 cout << "MaxZA: ";
488 cout << maxza << endl;
489
490 if (maxza < 30) { //colors=green to yellow
491 maxza=9*maxza/30;
492 maxza=80+maxza;
493
494 } else { //colors=yellow to red
495 maxza-=30;
496 maxza=11*maxza/60;
497 maxza=89+maxza;
498 }
499
500 tg->SetFillColor((Int_t)maxza);
501 tg->SetLineColor((Int_t)maxza);
502 tg->SetLineWidth((Int_t)maxza);
503 tg->Draw("AF");
504
505 tg->GetXaxis()->SetLimits(0,52);
506 tg->GetXaxis()->SetTickLength(0.1);
507 tg->GetXaxis()->SetNdivisions(0);
508 tg->GetYaxis()->SetNdivisions(202);
509
510 TPaveLabel* l= new TPaveLabel(2,46,35,88, fObject->GetObjectName());
511 l->SetBorderSize(0);
512 l->SetFillStyle(0);
513 l->SetTextAlign(12);
514 l->Draw();
515
516 if ((vistimestart<13)&&(vistimeend>0)) {
517 TString label;
518 label=months[vistimestart-1]+"-"+months[vistimeend-1];
519 TPaveLabel* l2= new TPaveLabel(35,46,50,88, label);
520 l2->SetBorderSize(0);
521 l2->SetFillStyle(0);
522 l2->SetTextAlign(32);
523 l2->Draw();
524
525 }
526
527 c2->Modified();
528 c2->Update();
529
530return kTRUE;
531}
532
533Double_t MVPPlotter::LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt)
534{
535/* Evaluates predicted LUNAR part of sky brightness, in
536 * V magnitudes per square arcsecond, following K. Krisciunas
537 * and B. E. Schaeffer (1991) PASP 103, 1033.
538 *
539 * moon_phase = Phase of the Moon, between 0. (no moon) and 1. (full moon),
540 * rho (deg) = separation of moon and object,
541 * altmoon (deg) = altitude of moon above horizon,
542 * alt (deg) = altitude of object above horizon
543 */
544
545 double kzen=1.;
546
547 double rho_rad = rho*fgDegToRad;
548 double alpha = 180.*(1. - moon_phase);
549 double Zmoon = (90. - altmoon)*fgDegToRad;
550 double Z = (90. - alt)*fgDegToRad;
551
552 double istar = -0.4*(3.84 + 0.026*fabs(alpha) + 4.0e-9*pow(alpha,4.)); /*eqn 20*/
553 istar = pow(10.,istar);
554 if(fabs(alpha) < 7.) /* crude accounting for opposition effect */
555 istar = istar * (1.35 - 0.05 * fabs(istar));
556 /* 35 per cent brighter at full, effect tapering linearly to
557 zero at 7 degrees away from full. mentioned peripherally in
558 Krisciunas and Scheafer, p. 1035. */
559 double fofrho = 229087. * (1.06 + cos(rho_rad)*cos(rho_rad));
560 if(fabs(rho) > 10.)
561 fofrho=fofrho+pow(10.,(6.15 - rho/40.)); /* eqn 21 */
562 else if (fabs(rho) > 0.25)
563 fofrho= fofrho+ 6.2e7 / (rho*rho); /* eqn 19 */
564 else fofrho = fofrho+9.9e8; /*for 1/4 degree -- radius of moon! */
565 double Xzm = sqrt(1.0 - 0.96*sin(Zmoon)*sin(Zmoon));
566 if(Xzm != 0.) Xzm = 1./Xzm;
567 else Xzm = 10000.;
568 double Xo = sqrt(1.0 - 0.96*sin(Z)*sin(Z));
569 if(Xo != 0.) Xo = 1./Xo;
570 else Xo = 10000.;
571 double Bmoon = fofrho * istar * pow(10.,(-0.4*kzen*Xzm))
572 * (1. - pow(10.,(-0.4*kzen*Xo))); /* nanoLamberts */
573 // cout << " Bmoon=" << Bmoon;
574 if(Bmoon > 0.001)
575 return(Bmoon);
576 else return(99999.);
577}
578
579
Note: See TracBrowser for help on using the repository browser.