Index: trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc (revision 1681) @@ -0,0 +1,69 @@ +/* ======================================================================== *\ +! +! * +! * 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 10/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MObservatoryLocation // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MObservatoryLocation.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MObservatoryLocation); + +MObservatoryLocation::MObservatoryLocation(const char *name, const char *title) +{ + fName = name ? name : "MObservatoryLocation"; + fTitle = title ? title : "Storage container for coordinates of an observatory"; + // TH1F fHorizon=new TH1F; + fgDegToRad=2*TMath::Pi()/360; + fLatitude = 28.7594 * fgDegToRad; // rad; 28 45 34 + fLongitude = 17.8761 * fgDegToRad; // rad; 17 52 34; + // slalib uses + for WEST !!! + fElevation = 2300; // m + fObsName = "Observatorio del Roque de los Muchachos"; +} + +//Double_t GetHorizon(Double_t phi); +//TF1 SetHorizonLine() { return fHorizon; } +//TF1 GetHorizonLine() { return fHorizon; } + +MObservatoryLocation::~MObservatoryLocation() +{ + // do nothing special. +} + +void MObservatoryLocation::Print(Option_t *) const +{ + *fLog << all; + *fLog << fObsName << endl; + *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude/fgDegToRad) : -(fLatitude/fgDegToRad)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl; + *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude/fgDegToRad) : -(fLongitude/fgDegToRad)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl; + *fLog << "Elevation " << fElevation << "m" << endl; +} Index: trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h (revision 1681) @@ -0,0 +1,44 @@ +#ifndef MARS_MObservatoryLocation +#define MARS_MObservatoryLocation + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MObservatoryLocation : public MParContainer +{ +private: + char* fObsName; + Double_t fLatitude; + Double_t fLongitude; + Double_t fElevation; + Double_t fgDegToRad; + // TH1F fHorizon; + +public: + MObservatoryLocation(const char *name=NULL, const char *title=NULL); + ~MObservatoryLocation(); + + inline void SetLatitude(Double_t latitude) { fLatitude = latitude; } + inline void SetLongitude(Double_t longitude) { fLongitude = longitude; } + inline void SetElevation(Double_t elevation) { fElevation = elevation; } + inline void SetObservatoryName(char* name) { fObsName = name; } + + void MObservatoryLocation::Print(Option_t *) const; + + inline Double_t GetLatitude() { return fLatitude/fgDegToRad; } + inline Double_t GetLongitude() { return fLongitude/fgDegToRad; } + inline Double_t GetElevation() { return fElevation; } + Double_t GetLatitudeRad() { return fLatitude; } + Double_t GetLongitudeRad() { return fLongitude; } + char* GetObservatoryName() { return fObsName; } + // Double_t GetHorizon(Double_t phi); + // void SetHorizonLine(TF1 hor) { fHorizon = hor; } + // TH1F GetHorizonLine() { return fHorizon; } + + ClassDef(MObservatoryLocation, 1) + +}; + +#endif + Index: trunk/MagicSoft/Mars/mtemp/MVPObject.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPObject.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MVPObject.h (revision 1681) @@ -0,0 +1,100 @@ +#ifndef MARS_MVPObject +#define MARS_MVPObject + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MObservatoryLocation +#include "MObservatoryLocation.h" +#endif + +#ifndef MARS_MVPTime +#include "MVPTime.h" +#endif + +class MVPObject : public MTask +{ +private: + MObservatoryLocation *fObservatory; + + Double_t fRA; + Double_t fDec; + Double_t fAlt; + Double_t fAz; + Double_t fZA; + Double_t fDiameter; + Double_t fEcLat; + Double_t fEcLong; + Bool_t fCalcEc; + + Double_t fUT1; + MVPTime *fTime; + + char* fObjectName; + UInt_t fBody; + + Double_t UT1ToGMST(Double_t ut1); + + Double_t fgDegToRad; + Double_t fgHrsToRad; + + Bool_t fGotRA; + Bool_t fGotDec; + Bool_t fGotName; + +public: + + MVPObject(const char *name=NULL, const char *title=NULL); + ~MVPObject(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + + void SetObservatory(MObservatoryLocation *observatory); + + void SetRA(Double_t rad) { fRA = rad; fBody = 10; fGotRA = kTRUE; } + void SetDec(Double_t rad) { fDec = rad; fBody = 10; fGotDec = kTRUE; } + + void SetRA(Int_t rh, Int_t rm, Int_t rs, Int_t ru = 0); + void SetDec(Int_t dh, Int_t dm, Int_t ds, Int_t du = 0); + + void SetCalcEc(Bool_t calcEc) { fCalcEc = calcEc; } + + void SetObjectName(char* name) { fObjectName = name; fGotName = kTRUE; } + void SetObjectByName(char* object); + + Bool_t SetObject(UInt_t body); + + Double_t GetRA() { return fRA/fgHrsToRad; } + Double_t GetDec() { return fDec/fgDegToRad; } + + Double_t GetRARad() { return fRA; } + Double_t GetDecRad() { return fDec; } + + Double_t GetZA() { return fZA; } + Double_t GetAltitude() { return fAlt; } + Double_t GetAzimut() { return fAz; } + Double_t GetDiameter() { return fDiameter; } + + Double_t GetEcLat() { return fEcLat; } + Double_t GetEcLong() { return fEcLong; } + + Double_t GetZADeg() { return fZA/fgDegToRad; } + Double_t GetAltitudeDeg() { return fAlt/fgDegToRad; } + Double_t GetAzimutDeg() { return fAz/fgDegToRad; } + + Double_t GetDistance(MVPObject* object); + Double_t GetDistanceRad(MVPObject* object); + + char* GetObjectName() { return fObjectName; } + + void Print(Option_t *) const; + + Double_t MJDStartOfYear(UInt_t year); + + ClassDef(MVPObject, 1) +}; + +#endif + Index: trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc (revision 1681) @@ -0,0 +1,579 @@ +/* ======================================================================== *\ +! +! * +! * 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]; + 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 = distance < fMars->GetDistance(fObject) ? distance : fMars->GetDistance(fObject); + distance = distance < fJupiter->GetDistance(fObject) ? distance : fJupiter->GetDistance(fObject); + distance = distance < fSaturn->GetDistance(fObject) ? 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)); + } + cout<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++; + + Int_t startday2 = fTime->MJDStartOfYear(year); + 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*maxza/30; + maxza=80+maxza; + + } else { //colors=yellow to red + maxza-=30; + maxza=11*maxza/60; + maxza=89+maxza; + } + + 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; + 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 = 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=fofrho+pow(10.,(6.15 - rho/40.)); /* eqn 21 */ + else if (fabs(rho) > 0.25) + fofrho= 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.); +} + + Index: trunk/MagicSoft/Mars/mtemp/MVPPlotter.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPPlotter.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MVPPlotter.h (revision 1681) @@ -0,0 +1,85 @@ +#ifndef MARS_MVPPlotter +#define MARS_MVPPlotter + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MParList +#include "MParList.h" +#endif + +#ifndef MARS_MTaskList +#include "MTaskList.h" +#endif + +#ifndef MARS_MVPObject +#include "MVPObject.h" +#endif + +#ifndef MARS_MVPTime +#include "MVPTime.h" +#endif + +#ifndef MARS_MObservatoryLocation +#include "MObservatoryLocation.h" +#endif + +class TH1D; +class TH2D; +class TCanvas; + +class MVPPlotter +{ +private: + + MVPTime *fTime; + MObservatoryLocation *fObs; + + MParList fPlist; + MTaskList fTlist; + + Bool_t fUseSun; + Bool_t fUseMoon; + Bool_t fUsePlanets; + + Double_t fAstronomicalDarkness; + + Int_t fgSecPerDay; + Int_t fgMJD010170; // 01-01-70 is JD 40586 + Double_t fgDegToRad; + + TH2D *fMjdUtcYear; + TH2D *fMjdUtcYearMoon; + TH1D *fMjdMoonPhase; + TH1D *fMjdMoonDistance; + TH1D *fMjdPlanetDistance; + TH1D *fMjdMoonIntensity; + TH1D *fMjdObsHours; + + MVPObject *fObject; + MVPObject *fSun; + MVPObject *fMoon; + MVPObject *fVenus; + MVPObject *fMars; + MVPObject *fJupiter; + MVPObject *fSaturn; + + void SetupObjects(); + +public: + MVPPlotter(const char *name=NULL, const char *title=NULL); + ~MVPPlotter(); + + void SetObject(MVPObject *o) { fObject=o; } + void SetObservatory(MObservatoryLocation *l) { fObs=l; } + void SetAstronomicalDarkness(Double_t d) { fAstronomicalDarkness=d; } + + Double_t LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt); + Bool_t CalcYear(Int_t year, UInt_t daySlices); + + ClassDef(MVPPlotter, 2) // Visibility Plotter: The Plotter Routine itself +}; + +#endif + Index: trunk/MagicSoft/Mars/mtemp/MVPTime.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPTime.cc (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MVPTime.cc (revision 1681) @@ -0,0 +1,125 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MVPTime // +// // +// Time Container used by the visibility plotter storing times in MJD // +// // +// This is a simple time container for the Visibility Plotter. It takes // +// "normal" or MJD-formatted time, converts it to MJD and can return MJD. // +// // +// The astronomical algorithms used are taken from SLALIB 2.4-8. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MVPTime.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MVPTime); + +MVPTime::MVPTime(const char *name, const char *title) +{ + fName = name ? name : "MVPTime"; + fTitle = title ? title : "Storage container for time"; +} + +MVPTime::~MVPTime() +{ + // do nothing special. +} + +void MVPTime::SetTime(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se) +{ + fUT1 = CalendarToUT1(yr, mo, dy, hr, mi, se); + fYear = yr; + fMonth = mo; + fDay = dy; +} + +void MVPTime::SetMJD(Double_t mjd, Double_t fracMjd) +{ + fUT1 = mjd + fracMjd; + Double_t fd; + Int_t status; + slaDjcl (mjd, &fYear, &fMonth, &fDay, &fd, &status); +} + +// -------------------------------------------------------------------------- +// +// +Double_t MVPTime::CalendarToUT1(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se) +{ + int status; + Double_t mjd; + slaCldj(yr, mo, dy, &mjd, &status); + switch (status) + { + case 1: + *fLog << "CalendarToMJD Warn: bad year" << endl; + break; + case 2: + *fLog << "CalendarToMJD Warn: bad month" << endl; + break; + case 3: + *fLog << "CalendarToMJD Warn: bad day" << endl; + break; + } + + mjd += ((Double_t)hr/24) + ((Double_t)mi/(24*60)) + ((Double_t)se/(24*60*60)); + return mjd; +} + +// -------------------------------------------------------------------------- +// +// +Double_t MVPTime::MJDStartOfYear(UInt_t year) +{ + int status; + Double_t mjd; + slaCldj(year, 1, 1, &mjd, &status); + switch (status) + { + case 1: + *fLog << "CalendarToMJD Warn: bad year" << endl; + break; + case 2: + *fLog << "CalendarToMJD Warn: bad month" << endl; + break; + case 3: + *fLog << "CalendarToMJD Warn: bad day" << endl; + break; + } + return mjd; +} + +void MVPTime::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Time (MJD) is: " << fUT1 << endl; +} Index: trunk/MagicSoft/Mars/mtemp/MVPTime.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPTime.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/MVPTime.h (revision 1681) @@ -0,0 +1,37 @@ +#ifndef MARS_MVPTime +#define MARS_MVPTime + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MVPTime : public MParContainer +{ +private: + Double_t fUT1; + Int_t fYear; + Int_t fMonth; + Int_t fDay; + Double_t MVPTime::CalendarToUT1(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se); + +public: + MVPTime(const char *name=NULL, const char *title=NULL); + ~MVPTime(); + + void SetTime(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se); + void SetMJD(Double_t mjd, Double_t fracMjd = 0); + inline Double_t GetMJD() { return fUT1; } + + void Print(Option_t *) const; + + Double_t MJDStartOfYear(UInt_t year); + Double_t GetYear() { return fYear; } + Double_t GetMonth() { return fMonth; } + Double_t GetDay() { return fDay; } + + ClassDef(MVPTime, 1) + +}; + +#endif + Index: trunk/MagicSoft/Mars/mtemp/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/Makefile (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/Makefile (revision 1681) @@ -0,0 +1,51 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = Temp + +# +# Library name to creatre +# +LIB = mtemp.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mraw -I../manalysis -I../mmc \ + -I../mgui -I../mgeom -I../mdata -I../mhist -I../../slalib + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = \ + MObservatoryLocation.cc \ + MVPObject.cc \ + MVPPlotter.cc \ + MVPTime.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode + + Index: trunk/MagicSoft/Mars/mtemp/TempIncl.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/TempIncl.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/TempIncl.h (revision 1681) @@ -0,0 +1,5 @@ +#ifndef __CINT__ + +#include + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mtemp/TempLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 1681) +++ trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 1681) @@ -0,0 +1,10 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MObservatoryLocation+; +#pragma link C++ class MVPObject+; +#pragma link C++ class MVPPlotter+; +#pragma link C++ class MVPTime+;