Index: trunk/MagicSoft/Mars/mtemp/MVPObject.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPObject.cc (revision 1680) +++ trunk/MagicSoft/Mars/mtemp/MVPObject.cc (revision 1680) @@ -0,0 +1,312 @@ +/* ======================================================================== *\ +! +! * +! * 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, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MVPObject // +// // +// Class used by the visibility plotter to convert RA/Dec to Alt/Az // +// // +// This class represents an object and is used with the Visibility // +// macro. It must be provided with its RA/Dec coordinates and an // +// object name (cf. MVPObject::SetRA, MVPObject::SetDec, and // +// MVPObject::SetName). Alternatively, you can require the MVPObject // +// to be a solar system object like the Sun, Mars or the Moon // +// (cf. MVPObject::SetObject). // +// // +// MVPObject is ready to be used in a Mars Eventloop. You must provide // +// an Observatory Location as well as a time at which the position // +// of the MVPObject is to be calculated. MVPObject::PreProcess // +// checks the existence of the required containers and also makes sure // +// all necessary setters have been called. MVPObject::Process // +// then calculates the Alt/Az position of the object, as well as the // +// Zenith angle and the object diameter (Solar system objects). // +// // +// The astronomical algorithms used are taken from SLALIB 2.4-8. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MVPObject.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" + +#include "../../slalib/slalib.h" + +ClassImp(MVPObject); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MVPObject::MVPObject(const char *name, const char *title) : fDiameter(0), fCalcEc(kFALSE), fUT1(52000), fBody(10), fGotRA(kFALSE), fGotDec(kFALSE), fGotName(kFALSE) +{ + fName = name ? name : "MVPObject"; + fTitle = title ? title : "Task to calculate Alt, Az of a given object"; + + fgDegToRad=2*TMath::Pi()/360; + fgHrsToRad=2*TMath::Pi()/24; +} + +MVPObject::~MVPObject() +{ + //Destructor: nothing special yet. +} + +// -------------------------------------------------------------------------- +// +// Check if necessary containers exist in the parameter list already. +// We need an ObservatoryLocation and a MVPTime object. +// +Bool_t MVPObject::PreProcess(MParList *pList) +{ + fObservatory = (MObservatoryLocation*)pList->FindObject("MObservatoryLocation"); + if (!fObservatory) + { + *fLog << dbginf << "MObservatoryLocation not found... aborting." << endl; + return kFALSE; + } + + fTime = (MVPTime*)pList->FindObject("MVPTime"); + if (!fTime) + { + *fLog << dbginf << "MVPTime not found... aborting." << endl; + return kFALSE; + } + + if (!fGotRA || !fGotDec || !fGotName) + { + *fLog << dbginf << "Object information is not complete." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets coordinates from object name. Instead of providing RA, Dec and Name +// of an object, you may also just provide the object name in the from +// HHMMsDDT, where RA is given in hours and minutes and Declination is +// given by degrees DD and tenths of degrees T. "s" may be "+" or +// "-" +// +void MVPObject::SetObjectByName(char* object) +{ + fObjectName=object; + fGotName=kTRUE; + +// cout<<"OBJ:"< 9) + { + *fLog << dbginf << "No solar system object associated with value " << body <<"! Ignoring request." << endl; + return kFALSE; + } + else // We are working on a solar system body. + { + switch (body) + { + case 1: fObjectName="Mercury"; break; + case 2: fObjectName="Venus"; break; + case 3: fObjectName="Moon"; break; + case 4: fObjectName="Mars"; break; + case 5: fObjectName="Jupiter"; break; + case 6: fObjectName="Saturn"; break; + case 7: fObjectName="Uranus"; break; + case 8: fObjectName="Neptune"; break; + case 9: fObjectName="Pluto"; break; + default: fObjectName="Sun"; + } + } + + fBody = body; + fGotRA = fGotDec = fGotName = kTRUE; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Given RA, Dec or a solar system object as well as an observatory +// location and a MVPTime, MVPObject::Process() calculates +// Alt, Az, ZA and (in the case of solar system objects) the apparent +// object diameter +// +Bool_t MVPObject::Process() +{ + Double_t diameter = 0.0; + + if (fBody < 10) // We are working on a solar system body. + { + slaRdplan(fTime->GetMJD(), fBody, fObservatory->GetLongitudeRad(), fObservatory->GetLatitudeRad(), &fRA, &fDec, &diameter); + } + + if (fCalcEc) slaEqecl(fRA, fDec, fTime->GetMJD(), &fEcLong, &fEcLat); + + Float_t azimuth; + Float_t elevation; + + Float_t hourAngle = (Float_t)UT1ToGMST(fTime->GetMJD()) - fRA; + + // cout << "ha: " << hourAngle << " ra: " << fRA << " dec " << fDec <GetLatitudeRad(), &azimuth, &elevation); + + fZA = slaZd(hourAngle, fDec, fObservatory->GetLatitudeRad()); + fAlt = (Double_t)elevation; + fAz = (Double_t)azimuth; + fDiameter = diameter; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Returns distance of given object to this object in degrees +// +Double_t MVPObject::GetDistance(MVPObject* object) +{ + return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad())/fgDegToRad; +} + +// -------------------------------------------------------------------------- +// +// Returns distance of given object to this object in radians +// +Double_t MVPObject::GetDistanceRad(MVPObject* object) +{ + return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad()); +} + + +// -------------------------------------------------------------------------- +// +// Converts UT1 (given as MJD) to Greenwich mean star time in radians +// +Double_t MVPObject::UT1ToGMST(Double_t ut1) +{ + return slaGmst(ut1); +} + +void MVPObject::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Position of "<< fObjectName << + ": Dec " << fDec/fgDegToRad << " deg, " << + "RA " << fRA/fgHrsToRad << " hrs" << endl; + if (fCalcEc) *fLog << "Ecliptic Long: " << fEcLong/fgDegToRad << " deg, " << + "Ecliptic Lat: " << fEcLat/fgDegToRad << " deg, " << endl; +} + + Index: trunk/MagicSoft/Mars/mtemp/MVPObject.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/MVPObject.h (revision 1679) +++ (revision ) @@ -1,100 +1,0 @@ -#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 1679) +++ (revision ) @@ -1,579 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1679) +++ (revision ) @@ -1,85 +1,0 @@ -#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 1679) +++ (revision ) @@ -1,125 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 1679) +++ (revision ) @@ -1,37 +1,0 @@ -#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/TempIncl.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/TempIncl.h (revision 1679) +++ (revision ) @@ -1,5 +1,0 @@ -#ifndef __CINT__ - -#include - -#endif // __CINT__ Index: trunk/MagicSoft/Mars/mtemp/TempLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 1679) +++ (revision ) @@ -1,19 +1,0 @@ -#ifdef __CINT__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ class MApplyPadding+; - -#pragma link C++ class MHSigmaPixel+; -#pragma link C++ class MHSigmabarTheta+; - -#pragma link C++ class MSigmabar+; -#pragma link C++ class MSigmabarCalc+; -#pragma link C++ class MSigmabarParam+; - -#pragma link C++ class MObservatoryLocation+; -#pragma link C++ class MVPObject+; -#pragma link C++ class MVPPlotter+; -#pragma link C++ class MVPTime+;