/* ======================================================================== *\
!
! *
! * 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];
  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;
    UInt_t i;
    obsHours=0;
    for (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 = 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<<endl;
  }
  
  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++;
  
  Int_t startday2 = fTime->MJDStartOfYear(year);
  Int_t stopday2  = fTime->MJDStartOfYear(year+1)-1;
  
    for (int day=startday2; day<stopday2+1; day=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*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.);
}


