source: tags/Mars-V0.9.4.1/mtemp/MVPPlotter.cc

/* ======================================================================== *\
!
! *
! * 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 <magicdev@rwagner.de> 12/2002
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MVPPlotter                                                          //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MVPPlotter.h"

#include <stream.h>

#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<stopday+1; day++) {     
    Double_t todaysPhase=0;
    Double_t moonIntensity;
    obsHours=0;
    for (UInt_t i=0; i<daySlices; i++)
      {  
        // Rearrange filling of bins such that a "day" doesn't start at midnight,
        // but rather at noon. This results in plots in which a "day" covers
        // a whole night.
        if (i>=(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" <<endl;


  /*!!!!!!!!!!!!!!!!!!!!!!!*/gROOT->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<stopday2+1; day+=7)
    {
        simpleCounter[weekCounter]=weekCounter-1;
        objectHeight[weekCounter]=0;
        for (int i=0; i<daySlices; i++)
        {
            if (i>=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.);
}