/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Robert Wagner 12/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MVPPlotter // // // ///////////////////////////////////////////////////////////////////////////// #include "MVPPlotter.h" #include #include "TCanvas.h" #include "TH1.h" #include "TH2.h" #include "TPaveText.h" #include "TPaveLabel.h" #include "TGraph.h" #include "TString.h" #include "TStyle.h" #include "MParList.h" #include "MVPTime.h" ClassImp(MVPPlotter); // -------------------------------------------------------------------------- // // Default constructor. // MVPPlotter::MVPPlotter(const char *name, const char *title) : fUseSun(kTRUE), fUseMoon(kTRUE), fUsePlanets(kTRUE), fAstronomicalDarkness(-18.0) { // fName = name ? name : "MVPPlotter"; // fTitle = title ? title : "Generates visibility histograms and information"; fgSecPerDay = 86400; fgMJD010170 = 40586; // 01-01-70 is JD 40586 fgDegToRad = 2*TMath::Pi()/360; } MVPPlotter::~MVPPlotter() // -------------------------------------------------------------------------- // // Destructor. Deletes objects if allocated. // { if (fUseSun) delete fSun; if (fUseMoon) delete fMoon; if (fUsePlanets) { delete fVenus; delete fMars; delete fJupiter; delete fSaturn; } } void MVPPlotter::SetupObjects() // -------------------------------------------------------------------------- // // Create instances for Sun, Moon, Planets if requested. // { fTime = new MVPTime(); fPlist.AddToList(fTime); fPlist.AddToList(fObs); if (fObject==NULL) { cout << "You didn't specify an object!" << endl; } fObject->PreProcess(&fPlist); cout << "Processing object " << fObject->GetObjectName() << endl; if (fUseSun) { fSun = new MVPObject(); fSun->SetObject(0); // Sun fSun->PreProcess(&fPlist); fSun->SetCalcEc(kTRUE); } if (fUseMoon) { fMoon = new MVPObject(); fMoon->SetObject(3); // Moon fMoon->PreProcess(&fPlist); fMoon->SetCalcEc(kTRUE); } if (fUsePlanets) { fVenus = new MVPObject(); fVenus->SetObject(2); fVenus->PreProcess(&fPlist); fVenus->SetCalcEc(kTRUE); fMars = new MVPObject(); fMars->SetObject(4); fMars->PreProcess(&fPlist); fMars->SetCalcEc(kTRUE); fJupiter = new MVPObject(); fJupiter->SetObject(5); fJupiter->PreProcess(&fPlist); fJupiter->SetCalcEc(kTRUE); fSaturn = new MVPObject(); fSaturn->SetObject(6); fSaturn->PreProcess(&fPlist); fSaturn->SetCalcEc(kTRUE); } } // -------------------------------------------------------------------------- // // Plots for a single object and a whole year are generated. // // Currently we do the following: // - Create a plot MJD vs UTC for one year // - Create a plot maxObjHeight vs Week for one year // - Create visibility tables for one year (ZA, hours) // Bool_t MVPPlotter::CalcYear(Int_t year, UInt_t daySlices=96) { SetupObjects(); UInt_t startday = (UInt_t)fTime->MJDStartOfYear(year); UInt_t stopday = (UInt_t)fTime->MJDStartOfYear(year+1)-1; cout << "Processing period MJD "<< startday << " to MJD "<< stopday << endl; ULong_t startdayROOT = fgSecPerDay * (startday-fgMJD010170); ULong_t stopdayROOT = fgSecPerDay * (stopday-fgMJD010170); // 2D Plot ZA vs MJD and UTC fMjdUtcYear = new TH2D("fMjdUtcYear", "Visibility of object", stopday-startday+1,startdayROOT,stopdayROOT, daySlices+1,-450,fgSecPerDay+450); // Observability hours per day MJD fMjdObsHours = new TH1D("fMjdObsHours", "Observation hours per day", stopday-startday+1,startdayROOT,stopdayROOT); if (fUseMoon) { // 2D Plot ZA of moon vs MJD and UTC fMjdUtcYearMoon =new TH2D("fMjdUtcYearMoon", "Moon ZA", stopday-startday+1,startdayROOT,stopdayROOT, daySlices+1,-450,fgSecPerDay+450); // Moon phase vs MJD fMjdMoonPhase =new TH1D("fMjdMoonPhase", "Moon phase", stopday-startday+1,startdayROOT,stopdayROOT); // Moon distance of object vs MJD fMjdMoonDistance =new TH1D("fMjdMoonDistance", "Moon distance", stopday-startday+1,startdayROOT,stopdayROOT); // Moon intensity at object vs MJD fMjdMoonIntensity =new TH1D("fMjdMoonIntensity", "Moon intensity at locus of object", stopday-startday+1,startdayROOT,stopdayROOT); } if (fUsePlanets) { // Distance of closest planet vs MJD fMjdPlanetDistance =new TH1D("fMjdPlanetDistance", "PlanetDistance", stopday-startday+1,startdayROOT,stopdayROOT); } // Array which holds total visibility time for whole year [0] and // each month [1]..[12] Float_t visibility[13][18]; memset(visibility, 0, 13*18*sizeof(Float_t)); /* for (int m=0;m<13;m++) for (int z=0;z<18;z++) visibility[m][z]=0; */ int fday, ftime; Double_t phase=0; Double_t obsHours; for (UInt_t day=startday; day=(daySlices/2)) { fday=day; ftime=i-(daySlices/2); } else { fday=day-1; ftime=i+(daySlices/2); } // Objects access fTime via parameter list... fTime->SetMJD(day,(Double_t)i/daySlices); if (fUseSun) fSun->Process(); if (fUseMoon) fMoon->Process(); if (fUseSun && fUseMoon) { // Calculate moon phase... phase = fSun->GetEcLong() - fMoon->GetEcLong(); phase = TMath::Pi()-(acos(cos(phase)*cos(fMoon->GetEcLat()))); phase = phase*180/TMath::Pi(); todaysPhase+=phase; } // If sun is not up (or we should not use sun information...) if (fSun->GetAltitudeDeg() < fAstronomicalDarkness) { // Calculate Position of object: fObject->Process(); // Calculate moon brightness at locus of object // now this is gonna be fun... /* Evaluates predicted LUNAR part of sky brightness, in * V magnitudes per square arcsecond, */ moonIntensity = LunSkyBright(phase/180, fObject->GetDistance(fMoon), fMoon->GetAltitudeDeg(), fObject->GetAltitudeDeg()); fMjdMoonIntensity->Fill(fgSecPerDay*(fday-fgMJD010170),moonIntensity); // If moon is not up (or we should not use moon information...) if (!fUseMoon || fMoon->GetAltitudeDeg()<=0 || moonIntensity<60) { // Fill MJD-UTC histogram fMjdUtcYear->Fill(fgSecPerDay*(fday-fgMJD010170),fgSecPerDay*ftime/daySlices,fObject->GetAltitudeDeg()); } // Sum up visibility time (only if moon not visible or moon // info shouldn't be used at all) if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) { // Calculate for (int z=0;z<18;z++) { if (fObject->GetAltitudeDeg()>(z*5)) { visibility[0][z]+=(Double_t)(60*24/daySlices); visibility[(Int_t)fTime->GetMonth()][z]+=(Double_t)(60*24/daySlices); } } //for if ((fObject->GetAltitudeDeg())>40) { fMjdObsHours->Fill(fgSecPerDay*(fday-fgMJD010170),(Double_t)(24/(Double_t)daySlices)); } } } //fi sun // Optional: Fill histo with moon-up times... // if (fMoon->GetAltitudeDeg() >0) // fMjdUtcYearMoon->Fill(fgSecPerDay*(day-fgMJD010170),fgSecPerDay*i/daySlices,phase); // fMoon->Process(); // Double_t phase; // phase = fSun->GetEcLong() - moon->GetEcLong(); // phase = TMath::Pi()-(acos(cos(phase)*cos(moon->GetEcLat()))); // cout << "Phase: " << phase*180/TMath::Pi() << endl; } //for daySlices // Distance of Venus to object if (fUsePlanets) { fVenus->Process(); fMars->Process(); fJupiter->Process(); fSaturn->Process(); Double_t distance = fVenus->GetDistance(fObject); distance = TMath::Min(distance, fMars->GetDistance(fObject)); distance = TMath::Min(distance, fJupiter->GetDistance(fObject)); distance = TMath::Min(distance, fSaturn->GetDistance(fObject)); fMjdPlanetDistance->Fill(fgSecPerDay*(fday-fgMJD010170),distance); } fMjdMoonPhase->Fill(fgSecPerDay*(fday-fgMJD010170),todaysPhase/i); fMjdMoonDistance->Fill(fgSecPerDay*(fday-fgMJD010170),fMoon->GetDistance(fObject)); } //for days // Here we print the tables with visibilities... cout << "Visibility time [hours]: " << endl; for (int z=1;z<17;z++) { if (visibility[0][z]==0) break; printf("Alt>%2d|%6d|",z*5,(Int_t)(visibility[0][z]/60)); for (int m=1;m<13;m++) { printf("%5d ",(Int_t)(visibility[m][z]/60)); } printf("\n"); } int vistimestart=0; int vistimeend=0; for (int m=1; m<13; m++) { int n = (m==1 ? 12 : m-1); if (visibility[m][8]/60>20 && visibility[n][8]/60<=20) { vistimestart=m; // we're on the rising slope } } for (int m=1; m<13; m++) { int n = (m==1 ? 12 : m-1); if (visibility[m][8]/60<20 && visibility[n][8]/60>=20) { vistimeend=n; // we're on the falling slope } } cout << "Visibility (Alt>40) during months: " << vistimestart << "-" << vistimeend << " approx " << visibility[0][8]/60 << "hrs" <Reset(); TCanvas *cMjdUtcYear = new TCanvas("cMjdUtcYear", "Object Visibility MjdUtcYear", 1100,500); cMjdUtcYear->cd(0); gStyle->SetPalette(1); // gStyle->SetOptStat(0); gStyle->SetOptFit(0); gStyle->SetFillStyle(0); gStyle->SetFillColor(10); gStyle->SetCanvasColor(10); gStyle->SetDrawBorder(0); gStyle->SetPadColor(10); gStyle->SetPadBorderSize(0); gStyle->SetPadLeftMargin(0.12); gStyle->SetTitleYOffset(1.2); gStyle->SetTitleXOffset(1.2); gStyle->SetPadLeftMargin(0.01); gStyle->SetPadRightMargin(0.1); gStyle->SetPadTopMargin(0.03); gStyle->SetPadBottomMargin(0.12); fMjdUtcYear->SetTitle(fObject->GetObjectName()); fMjdUtcYear->GetXaxis()->SetTimeDisplay(1); fMjdUtcYear->GetXaxis()->SetTimeFormat("%b %y"); fMjdUtcYear->GetYaxis()->SetTimeDisplay(1); fMjdUtcYear->GetYaxis()->SetTimeFormat("%Hh"); gStyle->SetTimeOffset(43200); fMjdUtcYear->GetYaxis()->SetLabelOffset(0.01); fMjdUtcYear->SetMinimum(0); fMjdUtcYear->SetMaximum(90); fMjdUtcYear->GetYaxis()->SetTitle("Hour"); fMjdUtcYear->GetXaxis()->SetTitle("Day of year"); fMjdUtcYear->GetZaxis()->SetTitle("Altitude"); fMjdUtcYear->Draw("SURF2BB9Z"); gPad->SetPhi(0); gPad->SetTheta(90); gPad->Modified(); TPaveText *pt = new TPaveText(-0.74,-0.35,-0.58,0.52,""); char ptLine[80]; pt->AddText("Visibility time [hrs]:"); pt->SetTextAlign(13); pt->SetFillStyle(0); pt->SetBorderSize(0); pt->SetTextFont(42); for (int j=1;j<17;j++) { sprintf (ptLine, "Alt>%i: %i", j*5, (Int_t)visibility[0][j]/60); pt->AddText(ptLine); if (visibility[0][j]==0) break; } pt->Draw(); if (fUseMoon) { TCanvas *cMjdMoon = new TCanvas("cMjdMoon", "the Moon phase", 600,600); gStyle->SetOptTitle(1); cMjdMoon->Divide(1,3); cMjdMoon->cd(1); fMjdMoonPhase->Draw(); cMjdMoon->cd(2); fMjdMoonDistance->Draw(); cMjdMoon->cd(3); fMjdMoonIntensity->Draw(); } TCanvas *cObsHours = new TCanvas("cObsHours", "ObsHours per day", 600,400); fMjdObsHours->Draw(); if (fUsePlanets) { TCanvas *cMjdPlanets = new TCanvas("cMjdPlanets", "the distance to planets", 600,200); cMjdPlanets->cd(0); fMjdPlanetDistance->Draw(); } // next histogram Float_t objectHeight[55]; Float_t simpleCounter[55]; Int_t weekCounter=0; simpleCounter[weekCounter]=0; objectHeight[weekCounter]=0; weekCounter++; const Int_t startday2 = fTime->MJDStartOfYear(year); const Int_t stopday2 = fTime->MJDStartOfYear(year+1)-1; for (int day=startday2; day=48) { fday=day; ftime=i-48; } else { fday=day-1; ftime=i+48; } fTime->SetMJD(day,(Double_t)i/daySlices); fObject->Process(); fSun->Process(); if (fSun->GetAltitudeDeg() < -25.0) { if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg())) objectHeight[weekCounter]=fObject->GetAltitudeDeg(); } } //i weekCounter++; } //day simpleCounter[weekCounter]=weekCounter-2; objectHeight[weekCounter]=0; weekCounter++; TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"}; TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100); // gStyle->SetOptTitle(0); // gStyle->SetPadLeftMargin(0.000001); // gStyle->SetPadRightMargin(0.000001); // gStyle->SetPadTopMargin(0.001); // gStyle->SetPadBottomMargin(0.000001); gPad->SetGrid(); c2->SetGrid(); TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight); tg->SetMinimum(1); tg->SetMaximum(90); Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude()); if (maxza > 90) maxza=90; if (maxza < 0) maxza=00; cout << "MaxZA: "; cout << maxza << endl; if (maxza < 30) { //colors=green to yellow maxza *= 9/30; maxza += 80; } else { //colors=yellow to red maxza -= 30; maxza *= 11/60; maxza += 89; } tg->SetFillColor((Int_t)maxza); tg->SetLineColor((Int_t)maxza); tg->SetLineWidth((Int_t)maxza); tg->Draw("AF"); tg->GetXaxis()->SetLimits(0,52); tg->GetXaxis()->SetTickLength(0.1); tg->GetXaxis()->SetNdivisions(0); tg->GetYaxis()->SetNdivisions(202); TPaveLabel* l= new TPaveLabel(2,46,35,88, fObject->GetObjectName()); l->SetBorderSize(0); l->SetFillStyle(0); l->SetTextAlign(12); l->Draw(); if ((vistimestart<13)&&(vistimeend>0)) { TString label=months[vistimestart-1]+"-"+months[vistimeend-1]; TPaveLabel* l2= new TPaveLabel(35,46,50,88, label); l2->SetBorderSize(0); l2->SetFillStyle(0); l2->SetTextAlign(32); l2->Draw(); } c2->Modified(); c2->Update(); return kTRUE; } Double_t MVPPlotter::LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt) { /* Evaluates predicted LUNAR part of sky brightness, in * V magnitudes per square arcsecond, following K. Krisciunas * and B. E. Schaeffer (1991) PASP 103, 1033. * * moon_phase = Phase of the Moon, between 0. (no moon) and 1. (full moon), * rho (deg) = separation of moon and object, * altmoon (deg) = altitude of moon above horizon, * alt (deg) = altitude of object above horizon */ double kzen=1.; double rho_rad = rho*fgDegToRad; double alpha = 180.*(1. - moon_phase); double Zmoon = (90. - altmoon)*fgDegToRad; double Z = (90. - alt)*fgDegToRad; double istar = -0.4*(3.84 + 0.026*fabs(alpha) + 4.0e-9*pow(alpha,4.)); /*eqn 20*/ istar = pow(10.,istar); if(fabs(alpha) < 7.) /* crude accounting for opposition effect */ istar *= 1.35 - 0.05 * fabs(istar); /* 35 per cent brighter at full, effect tapering linearly to zero at 7 degrees away from full. mentioned peripherally in Krisciunas and Scheafer, p. 1035. */ double fofrho = 229087. * (1.06 + cos(rho_rad)*cos(rho_rad)); if(fabs(rho) > 10.) fofrho+=pow(10.,(6.15 - rho/40.)); /* eqn 21 */ else if (fabs(rho) > 0.25) fofrho+=6.2e7 / (rho*rho); /* eqn 19 */ else fofrho = fofrho+9.9e8; /*for 1/4 degree -- radius of moon! */ double Xzm = sqrt(1.0 - 0.96*sin(Zmoon)*sin(Zmoon)); if(Xzm != 0.) Xzm = 1./Xzm; else Xzm = 10000.; double Xo = sqrt(1.0 - 0.96*sin(Z)*sin(Z)); if(Xo != 0.) Xo = 1./Xo; else Xo = 10000.; double Bmoon = fofrho * istar * pow(10.,(-0.4*kzen*Xzm)) * (1. - pow(10.,(-0.4*kzen*Xo))); /* nanoLamberts */ // cout << " Bmoon=" << Bmoon; if(Bmoon > 0.001) return(Bmoon); else return(99999.); }