#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.")
        ;

    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"));

    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");
    hframe.GetXaxis()->CenterTitle();
    hframe.GetYaxis()->CenterTitle();
    hframe.GetXaxis()->SetTimeDisplay(true);
    hframe.GetXaxis()->SetLabelSize(0.025);

    TCanvas c1;
    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->SetRightMargin(0.01);
    gPad->SetTopMargin(0.03);
    gPad->SetGrid();
    hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
    hframe.SetMinimum(0);
    hframe.SetMaximum(100);
    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;
        g1.SetName(name.data());
        g2.SetName(name.data());
        g1.SetLineWidth(2);
        g2.SetLineWidth(2);
        g1.SetLineStyle(style[cnt/6]);
        g2.SetLineStyle(style[cnt/6]);
        g1.SetLineColor(color[cnt%6]);
        g2.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);

            // Check if source is visible
            if (hrz.alt<15)
                continue;

            // Add point to curve
            const double axis = (mjd+h)*24*3600;
            g1.SetPoint(g1.GetN(), axis, hrz.alt);

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

            ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrz);

            const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec);

            const double lc = angle*hrz.alt*pow(disk, 6)/360/360;

            // Add point to curve
            g2.SetPoint(g2.GetN(), axis, lc>0 ? 7.7+4942*lc : 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);

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

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

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

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

    return 0;
}