Index: /trunk/FACT++/src/makeplots.cc
===================================================================
--- /trunk/FACT++/src/makeplots.cc	(revision 14450)
+++ /trunk/FACT++/src/makeplots.cc	(revision 14450)
@@ -0,0 +1,254 @@
+#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;
+}