Context Navigation

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

Visit:

Last change on this file since 3092 was 1745, checked in by tbretz, 22 years ago
* empty log message *
File size: 16.9 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
46	ClassImp(MVPPlotter);
47
48	// --------------------------------------------------------------------------
49	//
50	// Default constructor.
51	//
52	MVPPlotter::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
63	MVPPlotter::~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
79	void 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	//
143	Bool_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	memset(visibility, 0, 1318sizeof(Float_t));
193	/*
194	for (int m=0;m<13;m++)
195	for (int z=0;z<18;z++)
196	visibility[m][z]=0;
197	*/
198	int fday, ftime;
199	Double_t phase=0;
200	Double_t obsHours;
201
202	for (UInt_t day=startday; day<stopday+1; day++) {
203	Double_t todaysPhase=0;
204	Double_t moonIntensity;
205	obsHours=0;
206	for (UInt_t i=0; i<daySlices; i++)
207	{
208	// Rearrange filling of bins such that a "day" doesn't start at midnight,
209	// but rather at noon. This results in plots in which a "day" covers
210	// a whole night.
211	if (i>=(daySlices/2)) {
212	fday=day;
213	ftime=i-(daySlices/2);
214	} else {
215	fday=day-1;
216	ftime=i+(daySlices/2);
217	}
218
219	// Objects access fTime via parameter list...
220	fTime->SetMJD(day,(Double_t)i/daySlices);
221
222	if (fUseSun) fSun->Process();
223	if (fUseMoon) fMoon->Process();
224
225	if (fUseSun && fUseMoon) {
226
227	// Calculate moon phase...
228	phase = fSun->GetEcLong() - fMoon->GetEcLong();
229	phase = TMath::Pi()-(acos(cos(phase)*cos(fMoon->GetEcLat())));
230	phase = phase*180/TMath::Pi();
231	todaysPhase+=phase;
232
233	}
234
235	// If sun is not up (or we should not use sun information...)
236	if (fSun->GetAltitudeDeg() < fAstronomicalDarkness) {
237	// Calculate Position of object:
238	fObject->Process();
239
240	// Calculate moon brightness at locus of object
241	// now this is gonna be fun...
242
243	/* Evaluates predicted LUNAR part of sky brightness, in
244	* V magnitudes per square arcsecond,
245	*/
246
247	moonIntensity = LunSkyBright(phase/180, fObject->GetDistance(fMoon),
248	fMoon->GetAltitudeDeg(), fObject->GetAltitudeDeg());
249	fMjdMoonIntensity->Fill(fgSecPerDay*(fday-fgMJD010170),moonIntensity);
250
251
252
253	// If moon is not up (or we should not use moon information...)
254	if (!fUseMoon \|\| fMoon->GetAltitudeDeg()<=0 \|\| moonIntensity<60) {
255	// Fill MJD-UTC histogram
256	fMjdUtcYear->Fill(fgSecPerDay(fday-fgMJD010170),fgSecPerDayftime/daySlices,fObject->GetAltitudeDeg());
257	}
258
259	// Sum up visibility time (only if moon not visible or moon
260	// info shouldn't be used at all)
261	if ((!fUseMoon)\|\|(fMoon->GetAltitudeDeg()<=0)\|\|(moonIntensity<60)) {
262	// Calculate
263	for (int z=0;z<18;z++) {
264	if (fObject->GetAltitudeDeg()>(z*5)) {
265	visibility[0][z]+=(Double_t)(60*24/daySlices);
266	visibility[(Int_t)fTime->GetMonth()][z]+=(Double_t)(60*24/daySlices);
267	}
268	} //for
269
270	if ((fObject->GetAltitudeDeg())>40) {
271	fMjdObsHours->Fill(fgSecPerDay*(fday-fgMJD010170),(Double_t)(24/(Double_t)daySlices));
272	}
273
274	}
275
276
277	} //fi sun
278
279	// Optional: Fill histo with moon-up times...
280	// if (fMoon->GetAltitudeDeg() >0)
281	// fMjdUtcYearMoon->Fill(fgSecPerDay(day-fgMJD010170),fgSecPerDayi/daySlices,phase);
282	// fMoon->Process();
283	// Double_t phase;
284	// phase = fSun->GetEcLong() - moon->GetEcLong();
285	// phase = TMath::Pi()-(acos(cos(phase)*cos(moon->GetEcLat())));
286	// cout << "Phase: " << phase*180/TMath::Pi() << endl;
287
288	} //for daySlices
289
290	// Distance of Venus to object
291
292	if (fUsePlanets) {
293	fVenus->Process();
294	fMars->Process();
295	fJupiter->Process();
296	fSaturn->Process();
297
298	Double_t distance = fVenus->GetDistance(fObject);
299	distance = TMath::Min(distance, fMars->GetDistance(fObject));
300	distance = TMath::Min(distance, fJupiter->GetDistance(fObject));
301	distance = TMath::Min(distance, fSaturn->GetDistance(fObject));
302
303	fMjdPlanetDistance->Fill(fgSecPerDay*(fday-fgMJD010170),distance);
304	}
305
306	fMjdMoonPhase->Fill(fgSecPerDay*(fday-fgMJD010170),todaysPhase/i);
307	fMjdMoonDistance->Fill(fgSecPerDay*(fday-fgMJD010170),fMoon->GetDistance(fObject));
308
309
310	} //for days
311
312
313	// Here we print the tables with visibilities...
314	cout << "Visibility time [hours]: " << endl;
315
316	for (int z=1;z<17;z++) {
317	if (visibility[0][z]==0) break;
318	printf("Alt>%2d\|%6d\|",z*5,(Int_t)(visibility[0][z]/60));
319	for (int m=1;m<13;m++) {
320	printf("%5d ",(Int_t)(visibility[m][z]/60));
321	}
322	printf("\n");
323	}
324
325	int vistimestart=0;
326	int vistimeend=0;
327	for (int m=1; m<13; m++) {
328	int n = (m==1 ? 12 : m-1);
329	if (visibility[m][8]/60>20 && visibility[n][8]/60<=20) {
330	vistimestart=m; // we're on the rising slope
331	}
332	}
333
334	for (int m=1; m<13; m++) {
335	int n = (m==1 ? 12 : m-1);
336	if (visibility[m][8]/60<20 && visibility[n][8]/60>=20) {
337	vistimeend=n; // we're on the falling slope
338	}
339	}
340
341	cout << "Visibility (Alt>40) during months: " << vistimestart << "-" << vistimeend << " approx " << visibility[0][8]/60 << "hrs" <<endl;
342
343
344	/!!!!!!!!!!!!!!!!!!!!!!!/gROOT->Reset();
345
346	TCanvas *cMjdUtcYear = new TCanvas("cMjdUtcYear", "Object Visibility MjdUtcYear", 1100,500);
347	cMjdUtcYear->cd(0);
348
349	gStyle->SetPalette(1);
350	// gStyle->SetOptStat(0);
351	gStyle->SetOptFit(0);
352	gStyle->SetFillStyle(0);
353	gStyle->SetFillColor(10);
354	gStyle->SetCanvasColor(10);
355	gStyle->SetDrawBorder(0);
356	gStyle->SetPadColor(10);
357	gStyle->SetPadBorderSize(0);
358	gStyle->SetPadLeftMargin(0.12);
359	gStyle->SetTitleYOffset(1.2);
360	gStyle->SetTitleXOffset(1.2);
361
362	gStyle->SetPadLeftMargin(0.01);
363	gStyle->SetPadRightMargin(0.1);
364	gStyle->SetPadTopMargin(0.03);
365	gStyle->SetPadBottomMargin(0.12);
366
367	fMjdUtcYear->SetTitle(fObject->GetObjectName());
368	fMjdUtcYear->GetXaxis()->SetTimeDisplay(1);
369	fMjdUtcYear->GetXaxis()->SetTimeFormat("%b %y");
370	fMjdUtcYear->GetYaxis()->SetTimeDisplay(1);
371	fMjdUtcYear->GetYaxis()->SetTimeFormat("%Hh");
372	gStyle->SetTimeOffset(43200);
373	fMjdUtcYear->GetYaxis()->SetLabelOffset(0.01);
374	fMjdUtcYear->SetMinimum(0);
375	fMjdUtcYear->SetMaximum(90);
376
377	fMjdUtcYear->GetYaxis()->SetTitle("Hour");
378	fMjdUtcYear->GetXaxis()->SetTitle("Day of year");
379	fMjdUtcYear->GetZaxis()->SetTitle("Altitude");
380	fMjdUtcYear->Draw("SURF2BB9Z");
381	gPad->SetPhi(0);
382	gPad->SetTheta(90);
383	gPad->Modified();
384
385	TPaveText *pt = new TPaveText(-0.74,-0.35,-0.58,0.52,"");
386	char ptLine[80];
387	pt->AddText("Visibility time [hrs]:");
388	pt->SetTextAlign(13);
389	pt->SetFillStyle(0);
390	pt->SetBorderSize(0);
391	pt->SetTextFont(42);
392	for (int j=1;j<17;j++) {
393	sprintf (ptLine, "Alt>%i: %i", j*5, (Int_t)visibility[0][j]/60);
394	pt->AddText(ptLine);
395	if (visibility[0][j]==0) break;
396	}
397	pt->Draw();
398
399
400	if (fUseMoon) {
401	TCanvas *cMjdMoon = new TCanvas("cMjdMoon", "the Moon phase", 600,600);
402	gStyle->SetOptTitle(1);
403	cMjdMoon->Divide(1,3);
404	cMjdMoon->cd(1);
405	fMjdMoonPhase->Draw();
406	cMjdMoon->cd(2);
407	fMjdMoonDistance->Draw();
408	cMjdMoon->cd(3);
409	fMjdMoonIntensity->Draw();
410	}
411
412	TCanvas *cObsHours = new TCanvas("cObsHours", "ObsHours per day", 600,400);
413	fMjdObsHours->Draw();
414
415
416
417	if (fUsePlanets) {
418	TCanvas *cMjdPlanets = new TCanvas("cMjdPlanets", "the distance to planets", 600,200);
419	cMjdPlanets->cd(0);
420	fMjdPlanetDistance->Draw();
421	}
422
423
424	// next histogram
425	Float_t objectHeight[55];
426	Float_t simpleCounter[55];
427	Int_t weekCounter=0;
428
429	simpleCounter[weekCounter]=0;
430	objectHeight[weekCounter]=0;
431	weekCounter++;
432
433	const Int_t startday2 = fTime->MJDStartOfYear(year);
434	const Int_t stopday2 = fTime->MJDStartOfYear(year+1)-1;
435
436	for (int day=startday2; day<stopday2+1; day+=7)
437	{
438	simpleCounter[weekCounter]=weekCounter-1;
439	objectHeight[weekCounter]=0;
440	for (int i=0; i<daySlices; i++)
441	{
442	if (i>=48)
443	{
444	fday=day;
445	ftime=i-48;
446	}
447	else
448	{
449	fday=day-1;
450	ftime=i+48;
451	}
452
453	fTime->SetMJD(day,(Double_t)i/daySlices);
454	fObject->Process();
455	fSun->Process();
456
457	if (fSun->GetAltitudeDeg() < -25.0)
458	{
459	if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg()))
460	objectHeight[weekCounter]=fObject->GetAltitudeDeg();
461	}
462
463	} //i
464	weekCounter++;
465	} //day
466	simpleCounter[weekCounter]=weekCounter-2;
467	objectHeight[weekCounter]=0;
468	weekCounter++;
469
470	TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
471
472	TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100);
473
474	// gStyle->SetOptTitle(0);
475	// gStyle->SetPadLeftMargin(0.000001);
476	// gStyle->SetPadRightMargin(0.000001);
477	// gStyle->SetPadTopMargin(0.001);
478	// gStyle->SetPadBottomMargin(0.000001);
479	gPad->SetGrid();
480
481	c2->SetGrid();
482	TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight);
483
484	tg->SetMinimum(1);
485	tg->SetMaximum(90);
486
487	Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude());
488
489	if (maxza > 90) maxza=90;
490	if (maxza < 0) maxza=00;
491
492	cout << "MaxZA: ";
493	cout << maxza << endl;
494
495	if (maxza < 30) { //colors=green to yellow
496	maxza *= 9/30;
497	maxza += 80;
498
499	} else { //colors=yellow to red
500	maxza -= 30;
501	maxza *= 11/60;
502	maxza += 89;
503	}
504
505	tg->SetFillColor((Int_t)maxza);
506	tg->SetLineColor((Int_t)maxza);
507	tg->SetLineWidth((Int_t)maxza);
508	tg->Draw("AF");
509
510	tg->GetXaxis()->SetLimits(0,52);
511	tg->GetXaxis()->SetTickLength(0.1);
512	tg->GetXaxis()->SetNdivisions(0);
513	tg->GetYaxis()->SetNdivisions(202);
514
515	TPaveLabel* l= new TPaveLabel(2,46,35,88, fObject->GetObjectName());
516	l->SetBorderSize(0);
517	l->SetFillStyle(0);
518	l->SetTextAlign(12);
519	l->Draw();
520
521	if ((vistimestart<13)&&(vistimeend>0)) {
522	TString label=months[vistimestart-1]+"-"+months[vistimeend-1];
523	TPaveLabel* l2= new TPaveLabel(35,46,50,88, label);
524	l2->SetBorderSize(0);
525	l2->SetFillStyle(0);
526	l2->SetTextAlign(32);
527	l2->Draw();
528
529	}
530
531	c2->Modified();
532	c2->Update();
533
534	return kTRUE;
535	}
536
537	Double_t MVPPlotter::LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt)
538	{
539	/* Evaluates predicted LUNAR part of sky brightness, in
540	* V magnitudes per square arcsecond, following K. Krisciunas
541	* and B. E. Schaeffer (1991) PASP 103, 1033.
542	*
543	* moon_phase = Phase of the Moon, between 0. (no moon) and 1. (full moon),
544	* rho (deg) = separation of moon and object,
545	* altmoon (deg) = altitude of moon above horizon,
546	* alt (deg) = altitude of object above horizon
547	*/
548
549	double kzen=1.;
550
551	double rho_rad = rho*fgDegToRad;
552	double alpha = 180.*(1. - moon_phase);
553	double Zmoon = (90. - altmoon)*fgDegToRad;
554	double Z = (90. - alt)*fgDegToRad;
555
556	double istar = -0.4(3.84 + 0.026fabs(alpha) + 4.0e-9pow(alpha,4.)); /eqn 20*/
557	istar = pow(10.,istar);
558	if(fabs(alpha) < 7.) /* crude accounting for opposition effect */
559	istar = 1.35 - 0.05 fabs(istar);
560	/* 35 per cent brighter at full, effect tapering linearly to
561	zero at 7 degrees away from full. mentioned peripherally in
562	Krisciunas and Scheafer, p. 1035. */
563	double fofrho = 229087. * (1.06 + cos(rho_rad)*cos(rho_rad));
564	if(fabs(rho) > 10.)
565	fofrho+=pow(10.,(6.15 - rho/40.)); /* eqn 21 */
566	else if (fabs(rho) > 0.25)
567	fofrho+=6.2e7 / (rhorho); / eqn 19 */
568	else fofrho = fofrho+9.9e8; /for 1/4 degree -- radius of moon! /
569
570	double Xzm = sqrt(1.0 - 0.96sin(Zmoon)sin(Zmoon));
571
572	if(Xzm != 0.) Xzm = 1./Xzm;
573	else Xzm = 10000.;
574
575	double Xo = sqrt(1.0 - 0.96sin(Z)sin(Z));
576	if(Xo != 0.) Xo = 1./Xo;
577	else Xo = 10000.;
578
579	double Bmoon = fofrho * istar * pow(10.,(-0.4kzenXzm))
580	* (1. - pow(10.,(-0.4kzenXo))); /* nanoLamberts */
581	// cout << " Bmoon=" << Bmoon;
582	if(Bmoon > 0.001)
583	return(Bmoon);
584	else return(99999.);
585	}
586
587

Note: See TracBrowser for help on using the repository browser.

Download in other formats: