#include <iostream>
#include <iomanip.h>
#include <fstream>
#include <string>
#include <cmath>
#include <TObject.h>
#include <TROOT.h>
#include <TChain.h> 
#include <TFile.h>
#include <TH3D.h>
#include <TH1D.h>
#include <TLine.h>
#include <TLatex.h>
#include <TCut.h>
#include <TLegend.h>
#include <TCanvas.h>
#include "TRandom.h"

// TM July 2005
// AJH July 2005

void moon (long JD, long secs, double *dec, double *ra, double *alt, double *az, double *hAngle);
Int_t moon(){

	Local_to_Equa1 localfunction1;
	localfunction1.set_with_equatorial_coordinates(0.,0.);

/********************************************************
 *Reading in JPL Data
 *Begins Jan 1, 2005 00:00 increments 1 day for 1 year
 *In file moon4.dat
 ********************************************************/

	double moon_eph[365];
        ifstream myMoon;
        myMoon.open("moon4.dat");
        if (!myMoon)
                cout << "Error opening file: " << "moon4.dat";
        else {
                double temp = 0;
                for (int i = 0; i < 365; i++){
                        myMoon >> temp;
                        moon_eph[i] = (double)temp;
                        if ((i % 5) == 0){cout << endl;}
                }
                myMoon.close();
        }

/********************************************************
 *Get Date and time input				*
 *Calculate MJD						*
 *Set conversion factors				*
 ********************************************************/

  double year, month, day, hours, minutes, seconds, nsecs;
  cout << "Year:    "; cin >> year;
  cout << "Month:   "; cin >> month;
  cout << "Day:     "; cin >> day;
  cout << "Hours:   "; cin >> hours;
  cout << "Minutes: "; cin >> minutes;
  cout << "Seconds: "; cin >> seconds;
	
  Double_t time_jd = Date(year, month, day, hours, minutes, seconds).get_modified_julian();    
  Double_t RADDEG=TMath::RadToDeg(); //180.0 / Pi()
  Double_t degrad = TMath::DegToRad(); //Pi() / 180.0
  Double_t PI = TMath::Pi();
  Double_t deghour = 24./360.;

/*************************************************
 *Defining canvas and histograms
 *and setting some of their properties
 *************************************************/

	c1 = new TCanvas("c1", "Declination");
	c2 = new TCanvas("c2", "Declination Difference");
	TH1D *h_moon_jpl = new TH1D("JPL Dec", "Declination with JPL", 365, 0, 365);
	TH1D *h_moon_rdp = new TH1D("RDP Dec", "Declination with RDP", 365, 0, 365);
	TH1D *h_moon_alb = new TH1D("ALB Dec", "Declination with ALB", 365, 0, 365);
	TH1D *h_diff_jr =  new TH1D("JPL vs RDP", "Declination Difference", 365, 0, 365);
	TH1D *h_diff_ja =  new TH1D("JPL vs ALB", "Declination Difference", 365, 0, 365);
	TH1D *h_diff_ra =  new TH1D("RDP vs ALB", "Declination Difference", 365, 0, 365); 
	h_moon_rdp->SetLineColor(2);
	h_moon_alb->SetLineColor(3);
	h_diff_jr->SetLineColor(2);
	h_diff_ja->SetLineColor(3);
//	gStyle->SetOptStat(0);

/************************************************
 *Start calculating declinations and 
 *filling the histograms
 ************************************************/ 

	Double_t ra;
	Double_t dec;
	Double_t diam;
	Double_t hAngle;
	Double_t alt;
	Double_t az;
	double ddec;

	for (int a = 0; a < 365; a++) {
		h_moon_jpl->Fill(a, moon_eph[a]);

		localfunction1.set_Moon(time_jd + a);
                localfunction1.get_Moon(ra,dec,diam);
		h_moon_rdp->Fill(a, dec * RADDEG);
		ddec = dec;
	
		moon(time_jd + a, 0, &dec, &ra, &alt, &az, &hAngle);
		h_moon_alb->Fill(a, dec);

		h_diff_jr->Fill(a, moon_eph[a]-ddec*RADDEG);
		h_diff_ja->Fill(a, moon_eph[a]-dec);
		h_diff_ra->Fill(a, ddec*RADDEG-dec);
	}

/**********************************************
 *Set the Legend and final display settings
 **********************************************/
	leg = new TLegend(0.4,0.6,0.89,0.89);
	leg->AddEntry(h_diff_jr,"JPL vs RDP","l");
	leg->AddEntry(h_diff_ja,"JPL vs ALB","l");
	leg->AddEntry(h_diff_ra,"RDP vs ALB","l");
	leg->SetHeader("Data Comparaisons");


/**********************************************
 *Draw and Display the histograms
 **********************************************/

	c1->cd();
	h_moon_jpl->Draw();
	h_moon_rdp->Draw("same");
	h_moon_alb->Draw("same");
	c1->Update();

	c2->cd();
	h_diff_ja->Draw();
	h_diff_jr->Draw("same");
	h_diff_ra->Draw("same");
	leg->Draw();
	c2->Update();

  return 1;
}