#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 CheckForGap(TCanvas &c, TGraph &g, double axis)
{
    if (g.GetN()==0 || axis-g.GetX()[g.GetN()-1]<450)
        return;

    c.cd();
    ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
    while (g.GetN())
        g.RemovePoint(0);
}

void DrawClone(TCanvas &c, TGraph &g)
{
    if (g.GetN()==0)
        return;

    c.cd();
    ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
}

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

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"
        "Calculates several plots for the sources in the database\n"
        "helpful or needed for scheduling. The Plot is always calculated\n"
        "for the night which starts at the same so. So no matter if\n"
        "you specify '1974-09-09 00:00:00' or '1974-09-09 21:56:00'\n"
        "the plots will refer to the night 1974-09-09/1974-09-10.\n"
        "The advantage is that specification of the date as in\n"
        "1974-09-09 is enough. Time axis starts and ends at nautical\n"
        "twilight which is 12deg below horizon.\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");

    // -12: nautical
    ln_rst_time sun_set;   // Sun set with the same date than th provided date
    ln_rst_time sun_rise;  // Sun rise on the following day
    ln_get_solar_rst_horizon(time.JD()-0.5, &observer, -12, &sun_set);
    ln_get_solar_rst_horizon(time.JD()+0.5, &observer, -12, &sun_rise);

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

    cout << "Time: " << time     << endl;
    cout << "Base: " << Time(jd-0.5) << endl;
    cout << "Set:  " << Time(sun_set.set)   << endl;
    cout << "Rise: " << Time(sun_rise.rise) << endl;

    const double jd0 = fmod(sun_set.set,   1);   // ~0.3
    const double jd1 = fmod(sun_rise.rise, 1);   // ~0.8

    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%F1995-01-01 00:00:00 GMT");
    hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+"  -  "+Time(jd+1).GetAsStr("%d/%m/%Y")+"  [UTC]").c_str());
    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;
    c1.SetFillColor(kWhite);
    c1.SetBorderMode(0);
    c1.SetFrameBorderMode(0);
    c1.SetLeftMargin(0.085);
    c1.SetRightMargin(0.01);
    c1.SetTopMargin(0.03);
    c1.SetGrid();
    hframe.GetYaxis()->SetTitle("Altitude [deg]");
    hframe.SetMinimum(15);
    hframe.SetMaximum(90);
    hframe.DrawCopy();

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

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

    TCanvas c4;
    c4.SetFillColor(kWhite);
    c4.SetBorderMode(0);
    c4.SetFrameBorderMode(0);
    c4.SetLeftMargin(0.085);
    c4.SetRightMargin(0.01);
    c4.SetTopMargin(0.03);
    c4.SetGrid();
    hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
    hframe.SetMinimum(0);
    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, g4, gr, gm;
        g1.SetName(name.data());
        g2.SetName(name.data());
        g3.SetName(name.data());
        g4.SetName(name.data());
        g1.SetLineWidth(2);
        g2.SetLineWidth(2);
        g3.SetLineWidth(2);
        g4.SetLineWidth(2);
        gm.SetLineWidth(1);
        g1.SetLineStyle(style[cnt/6]);
        g2.SetLineStyle(style[cnt/6]);
        g3.SetLineStyle(style[cnt/6]);
        g4.SetLineStyle(style[cnt/6]);
        g1.SetLineColor(color[cnt%6]);
        g2.SetLineColor(color[cnt%6]);
        g3.SetLineColor(color[cnt%6]);
        g4.SetLineColor(color[cnt%6]);
        gm.SetLineColor(kYellow);

        if (cnt==0)
            leg.AddEntry(gm.Clone(), "Moon", "l");
        leg.AddEntry(g1.Clone(), name.data(), "l");

        // 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 cang = sin(angle   *M_PI/180);
            const double calt = sin(hrzm.alt*M_PI/180);

            const double lc = cang==0 ? -1 : calt*pow(disk, 3)/sqrt(cang);
            const double cur = lc>0 ? 8.1+94.6*lc : 8.1;

            // 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
            CheckForGap(c1, g1, axis);
            CheckForGap(c1, gm, axis);
            CheckForGap(c2, g2, axis);
            CheckForGap(c3, g3, axis);
            CheckForGap(c4, g4, axis);

            // 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/8.1);

            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
                g4.SetPoint(g4.GetN(), axis, angle);

            if (cnt==0)
                gm.SetPoint(gm.GetN(), axis, hrzm.alt);
        }

        if (cnt==0)
            DrawClone(c1, gm);

        DrawClone(c1, g1);
        DrawClone(c2, g2);
        DrawClone(c3, g3);
        DrawClone(c4, g4);
    }


    // Save three plots
    TCanvas c5;
    c5.SetFillColor(kWhite);
    c5.SetBorderMode(0);
    c5.SetFrameBorderMode(0);
    leg.Draw();

    const string t = Time(jd).GetAsStr("%Y%m%d");

    c1.SaveAs((t+"-ZenithDistance.eps").c_str());
    c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
    c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
    c4.SaveAs((t+"-MoonDist.eps").c_str());
    c5.SaveAs((t+"-Legend.eps").c_str());

    c1.SaveAs((t+"-ZenithDistance.root").c_str());
    c2.SaveAs((t+"-PredictedCurrent.root").c_str());
    c3.SaveAs((t+"-RelativeThreshold.root").c_str());
    c4.SaveAs((t+"-MoonDist.root").c_str());

    return 0;
}