#include <libnova/solar.h>
#include <libnova/lunar.h>
#include <libnova/rise_set.h>
#include <libnova/transform.h>

#include "Database.h"

#include "Time.h"
#include "Configuration.h"

#include <TROOT.h>
#include <TH1.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TLegend.h>

using namespace std;

// ------------------------------------------------------------------------

double Angle(double ra, double dec, double r, double d)
{
    const double theta0 = M_PI/2-d*M_PI/180;
    const double phi0   = r*M_PI/12;

    const double theta1 = M_PI/2-dec*M_PI/180;
    const double phi1   = ra*M_PI/12;

    const double x0 = sin(theta0) * cos(phi0);
    const double y0 = sin(theta0) * sin(phi0);
    const double z0 = cos(theta0);

    const double x1 = sin(theta1) * cos(phi1);
    const double y1 = sin(theta1) * sin(phi1);
    const double z1 = cos(theta1);

    double arg = x0*x1 + y0*y1 + z0*z1;
    if(arg >  1.0) arg =  1.0;
    if(arg < -1.0) arg = -1.0;

    return acos(arg) * 180/M_PI;
}

// ========================================================================
// ========================================================================
// ========================================================================

void SetupConfiguration(Configuration &conf)
{
    po::options_description control("Smart FACT");
    control.add_options()
        ("ra",        var<double>(), "Source right ascension")
        ("dec",       var<double>(), "Source declination")
        ("date-time", var<string>(), "SQL time (UTC)")
        ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
        ("max-current", var<double>(75), "Maximum current to display in other plots.")
        ("max-zd", var<double>(75), "Maximum zenith distance to display in other plots")
        ("no-limits", po_switch(), "Switch off limits in plots")
        ;

    po::positional_options_description p;
    p.add("date-time", 1); // The first positional options

    conf.AddOptions(control);
    conf.SetArgumentPositions(p);
}

void PrintUsage()
{
    cout <<
        "makeplots - The astronomy plotter\n"
        "\n"
        "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
    cout << endl;
}

int main(int argc, const char* argv[])
{
    gROOT->SetBatch();

    Configuration conf(argv[0]);
    conf.SetPrintUsage(PrintUsage);
    SetupConfiguration(conf);

    if (!conf.DoParse(argc, argv))
        return 127;

    if (conf.Has("ra")^conf.Has("dec"))
    {
        cout << "ERROR - Either --ra or --dec missing." << endl;
        return 1;
    }

    // ------------------ Eval config ---------------------

    const double lon = -(17.+53./60+26.525/3600);
    const double lat =   28.+45./60+42.462/3600;

    ln_lnlat_posn observer;
    observer.lng = lon;
    observer.lat = lat;

    Time time;
    if (conf.Has("date-time"))
        time.SetFromStr(conf.Get<string>("date-time"));

    const double max_current = conf.Get<double>("max-current");
    const double max_zd      = conf.Get<double>("max-zd");
    const double no_limits   = conf.Get<bool>("no-limits");

    ln_rst_time sun_astronomical;
    ln_get_solar_rst_horizon(time.JD(), &observer, -12, &sun_astronomical);

    const double jd  = floor(time.JD());
    const double mjd = floor(time.Mjd())+49718-0.5;

    const double jd0 = fmod(sun_astronomical.set,  1);
    const double jd1 = fmod(sun_astronomical.rise, 1);

    cout << "JD0=" << jd+jd0 << endl;
    cout << "JD1=" << jd+jd1 << endl;

    const string fDatabase = conf.Get<string>("source-database");

    // ------------------ Precalc moon coord ---------------------

    vector<pair<ln_equ_posn, double>> fMoonCoords;

    for (double h=0; h<1; h+=1./(24*12))
    {
        ln_equ_posn  moon;
        ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);

        const double disk = ln_get_lunar_disk(jd+h);

        fMoonCoords.push_back(make_pair(moon, disk));
    }

    // ------------- Get Sources from databasse ---------------------

    const mysqlpp::StoreQueryResult res =
        Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM source").store();

    // ------------- Create canvases and frames ---------------------

    TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
    hframe.SetStats(kFALSE);
    hframe.GetXaxis()->SetTimeFormat("%Hh%M'");
    hframe.GetXaxis()->SetTitle("Time [UTC]");
    hframe.GetXaxis()->CenterTitle();
    hframe.GetYaxis()->CenterTitle();
    hframe.GetXaxis()->SetTimeDisplay(true);
    hframe.GetYaxis()->SetTitleSize(0.040);
    hframe.GetXaxis()->SetTitleSize(0.040);
    hframe.GetXaxis()->SetTitleOffset(1.1);
    hframe.GetYaxis()->SetLabelSize(0.040);
    hframe.GetXaxis()->SetLabelSize(0.040);

    TCanvas c1;
    gPad->SetLeftMargin(0.085);
    gPad->SetRightMargin(0.01);
    gPad->SetTopMargin(0.03);
    gPad->SetGrid();
    hframe.GetYaxis()->SetTitle("Altitude [deg]");
    hframe.SetMinimum(15);
    hframe.SetMaximum(90);
    hframe.DrawCopy();

    TCanvas c2;
    gPad->SetLeftMargin(0.085);
    gPad->SetRightMargin(0.01);
    gPad->SetTopMargin(0.03);
    gPad->SetGrid();
    hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
    hframe.SetMinimum(0);
    hframe.SetMaximum(100);
    hframe.DrawCopy();

    TCanvas c3;
    gPad->SetLeftMargin(0.085);
    gPad->SetRightMargin(0.01);
    gPad->SetTopMargin(0.03);
    gPad->SetGrid();
    gPad->SetLogy();
    hframe.GetYaxis()->SetTitle("Estimated relative threshold");
    hframe.SetMinimum(0.9);
    hframe.SetMaximum(180);
    hframe.DrawCopy();

    Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
    Int_t style[] = { kSolid, kDashed, kDotted };

    TLegend leg(0, 0, 1, 1);

    // ------------- Loop over sources ---------------------

    Int_t cnt=0;
    for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
    {
        // Eval row
        const string name = (*v)[0].c_str();

        ln_equ_posn pos;
        pos.ra  = double((*v)[1])*15;
        pos.dec = double((*v)[2]);

        // Create graphs
        TGraph g1, g2, g3;
        g1.SetName(name.data());
        g2.SetName(name.data());
        g3.SetName(name.data());
        g1.SetLineWidth(2);
        g2.SetLineWidth(2);
        g3.SetLineWidth(2);
        g1.SetLineStyle(style[cnt/6]);
        g2.SetLineStyle(style[cnt/6]);
        g3.SetLineStyle(style[cnt/6]);
        g1.SetLineColor(color[cnt%6]);
        g2.SetLineColor(color[cnt%6]);
        g3.SetLineColor(color[cnt%6]);

        // Loop over 24 hours
        int i=0;
        for (double h=0; h<1; h+=1./(24*12), i++)
        {
            // check if it is betwene sun-rise and sun-set
            if (h<jd0 || h>jd1)
                continue;

            // get local position of source
            ln_hrz_posn hrz;
            ln_get_hrz_from_equ(&pos, &observer, jd+h, &hrz);

            // Get moon properties and
            ln_equ_posn moon  = fMoonCoords[i].first;
            const double disk = fMoonCoords[i].second;

            ln_hrz_posn hrzm;
            ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrzm);

            // Distance between source and moon
            const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec);

            // Current prediction
            const double lc = angle*hrzm.alt*pow(disk, 6)/360/360;
            const double cur = lc>0 ? 7.7+4942*lc : 7.7;

            // Relative  energy threshold prediction
            const double cs = cos((90+hrz.alt)*M_PI/180);
            const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;

            // Add points to curve
            const double axis = (mjd+h)*24*3600;

            // If there is a gap of more than one bin, start a new curve
            if (g1.GetN()>0 && axis-g1.GetX()[g1.GetN()-1]>450)
            {
                c1.cd();
                ((TGraph*)g1.DrawClone("C"))->SetBit(kCanDelete);
                while (g1.GetN())
                    g1.RemovePoint(0);
            }

            if (g2.GetN()>0 && axis-g2.GetX()[g2.GetN()-1]>450)
            {
                c2.cd();
                ((TGraph*)g2.DrawClone("C"))->SetBit(kCanDelete);
                while (g2.GetN())
                    g2.RemovePoint(0);
            }

            if (g3.GetN()>0 && axis-g3.GetX()[g3.GetN()-1]>450)
            {
                c3.cd();
                ((TGraph*)g3.DrawClone("C"))->SetBit(kCanDelete);
                while (g3.GetN())
                    g3.RemovePoint(0);
            }

            // Add data
            if (no_limits || cur<max_current)
                g1.SetPoint(g1.GetN(), axis, hrz.alt);

            if (no_limits || 90-hrz.alt<max_zd)
                g2.SetPoint(g2.GetN(), axis, cur);

            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
                g3.SetPoint(g3.GetN(), axis, ratio*cur/7.7);
        }

        // Add graphs to canvases and add corresponding entry to legend
        TGraph *g = 0;

        c1.cd();
        g = (TGraph*)g1.DrawClone("C");
        g->SetBit(kCanDelete);

        c2.cd();
        g = (TGraph*)g2.DrawClone("C");
        g->SetBit(kCanDelete);

        c3.cd();
        g = (TGraph*)g3.DrawClone("C");
        g->SetBit(kCanDelete);

        leg.AddEntry(g, name.data(), "l");
    }

    // Save three plots
    TCanvas c4;
    leg.Draw();

    c1.SaveAs("test1.eps");
    c2.SaveAs("test2.eps");
    c3.SaveAs("test3.eps");
    c4.SaveAs("legend.eps");

    c1.SaveAs("test1.root");
    c2.SaveAs("test2.root");
    c3.SaveAs("test3.root");

    return 0;
}