Index: /trunk/FACT++/Makefile.am
===================================================================
--- /trunk/FACT++/Makefile.am	(revision 18410)
+++ /trunk/FACT++/Makefile.am	(revision 18411)
@@ -69,5 +69,5 @@
 if HAS_SQL
 if HAS_NOVA
-bin_PROGRAMS += makedata
+bin_PROGRAMS += makedata makeschedule
 if HAS_ROOT
 bin_PROGRAMS += makeplots
@@ -332,4 +332,7 @@
 makedata_LDADD = libTime.la libConfiguration.la
 
+makeschedule_SOURCES = src/makeschedule.cc
+makeschedule_LDADD = libTime.la libConfiguration.la libTools.la
+
 
 chatserv_SOURCES = src/chatserv.cc src/LocalControl.h
Index: /trunk/FACT++/src/makeschedule.cc
===================================================================
--- /trunk/FACT++/src/makeschedule.cc	(revision 18411)
+++ /trunk/FACT++/src/makeschedule.cc	(revision 18411)
@@ -0,0 +1,772 @@
+#include "externals/Prediction.h"
+
+//#include <unordered_map>
+#include <boost/algorithm/string/join.hpp>
+
+#include "Database.h"
+
+#include "tools.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;
+using namespace Nova;
+
+// -----------------------------------------------------------------------
+
+/*
+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("Makeschedule");
+    control.add_options()
+        ("date", var<string>(), "SQL time (UTC), e.g. '2016-12-24'")
+        ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
+        ("max-current", var<double>(90), "Global maximum current limit in uA")
+        ("max-zd", var<double>(75), "Global zenith distance limit in degree")
+        ("source", vars<string>(), "List of all TeV sources to be included, names according to the database")
+        ("setup.*", var<double>(), "Setup for the sources to be observed")
+        ("preobs.*", vars<string>(), "Prescheduled observations")
+        ("startup.offset", var<double>(15), "Determines how many minutes the startup is scheduled before data-taking.start [0;120]")
+        ("data-taking.start", var<double>(-12), "Begin of data-taking in degree of sun below horizon")
+        ("data-taking.end", var<double>(-13.75), "End of data-taking in degree of sun below horizon")
+        ;
+
+    po::positional_options_description p;
+    p.add("date", 1); // The first positional options
+
+    conf.AddOptions(control);
+    conf.SetArgumentPositions(p);
+}
+
+void PrintUsage()
+{
+    cout <<
+        "makeschedule - Creates an automatic schedule\n"
+        "\n"
+        "Usage: makeschedule [yyyy-mm-dd]\n";
+    cout << endl;
+}
+
+void PrintHelp()
+{
+#ifdef HAVE_ROOT
+    cout <<
+        "\n"
+        "Examples:\n"
+        "\n"
+        "  makeschedule 2016-12-24\n"
+        "\n"
+        "Calculate the Christmas schedule for 2016 using all TeV sources from the\n"
+        "database. If the date is omitted the current date is used.\n"
+        "\n"
+        "  makeschedule --source='Mrk 421' --source='Mrk 501' --source='Crab'\n"
+        "\n"
+        "Use only the mentioned sources to calculate the schedule.\n"
+        "\n"
+        "  makeschedule --source=Crab --setup.Crab.max-zd=30\n"
+        "\n"
+        "Limit the zenith distance of Crab into the range [0;30]deg.\n"
+        "\n"
+        "  makeschedule --source=Crab --setup.Crab.max-current=50\n"
+        "\n"
+        "Limit the maximum estimated current of Crab at 50uA.\n"
+        "\n"
+        "  makeschedule --source='IC 310' '--setup.IC 310.penalty=1.2'\n"
+        "\n"
+        "Multiply IC310's estimated relative threshold by a factor 1.2\n"
+        "\n";
+    cout << endl;
+#endif
+}
+
+
+struct MyDouble
+{
+    double val;
+    bool valid;
+    MyDouble(Configuration &conf, const string &str) : val(0)
+    {
+        valid = conf.Has(str);
+        if (valid)
+            val = conf.Get<double>(str);
+    }
+    MyDouble() : val(0), valid(false) {}
+};
+
+/*
+struct MinMax
+{
+    MyDouble min;
+    MyDouble max;
+    MinMax(Configuration &conf, const string &str)
+    {
+        min = MyDouble(conf, str+".min");
+        max = MyDouble(conf, str+".max");
+    }
+    MinMax() {}
+};
+*/
+
+struct Source
+{
+    // Global limits
+    static double max_current;
+    static double max_zd;
+
+    // Source descrition
+    string name;
+    EquPosn equ;
+
+    // Source specific limits
+    MyDouble maxzd;
+    MyDouble maxcurrent;
+    double penalty;
+
+    // Possible observation time
+    double begin;
+    double end;
+
+    // Threshold (sorting reference)
+    double threshold;
+
+    double duration() const { return end-begin; };
+
+    // Pre-observations (e.g. ratescan)
+    vector<string> preobs;
+
+    Source(const string &n="") : name(n), begin(0), threshold(std::numeric_limits<double>::max()) { }
+
+    //bool IsSpecial() const { return threshold==std::numeric_limits<double>::max(); }
+
+    double zd(const double &jd) const
+    {
+        return 90-GetHrzFromEqu(equ, jd).alt;
+    }
+
+    bool valid(const SolarObjects &so) const
+    {
+        const HrzPosn hrz     = GetHrzFromEqu(equ, so.fJD);
+        const double  current = FACT::PredictI(so, equ);
+
+        if (current>max_current)
+            return false;
+
+        if (hrz.alt<=0 || 90-hrz.alt>max_zd)
+            return false;
+
+        if (maxzd.valid && 90-hrz.alt>maxzd.val)
+            return false;
+
+        if (maxcurrent.valid && current>maxcurrent.val)
+            return false;
+
+        return true;
+    }
+
+    bool IsRangeValid(const double &jd_begin, const double &jd_end) const
+    {
+        const double step = 1./24/60;
+        for (double jd=jd_begin; jd<jd_end+step/2; jd++)
+            if (!valid(SolarObjects(jd)))
+                return false;
+
+        return true;
+    }
+
+    double getThreshold(const SolarObjects &so) const
+    {
+        const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD);
+        const double  current = FACT::PredictI(so, equ);
+
+        const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
+
+        return penalty*ratio*pow(current/6.2, 0.394);
+    }
+};
+
+double Source::max_zd;
+double Source::max_current;
+
+bool SortByThreshold(const Source &i, const Source &j) { return i.threshold<j.threshold; }
+
+bool RescheduleFirstSources(vector<Source> &obs)
+{
+    if (obs.size()<2 || obs[0].duration()>=40./24/60 || obs[0].name=="SLEEP" || obs[1].name=="SLEEP")
+        return false;
+
+    cout << "First source [" << obs[0].name << "] detected < 40min" << endl;
+
+    const double obs1_duration = obs[1].end - obs[0].begin - 40./24/60;
+    const double obs0_end      =              obs[0].begin + 40./24/60;
+
+    // Check that:
+    //  - the duration for the shrunken source obs[1] is still above 40min
+    //  - obs[0] does not exceed 60deg at the end of its new window
+    //  - obs[0] does not exceed any limit within the new window
+
+    if (obs1_duration>=40./24/60 && obs[0].IsRangeValid(obs[0].end, obs0_end))
+    {
+        obs[0].end   = obs0_end;
+        obs[1].begin = obs0_end;
+
+        cout << "First source [" << obs[0].name << "] extended to 40min" << endl;
+
+        return false;
+    }
+
+    // Try to remove the first source, check if second source fullfills all limits
+    if (obs[1].IsRangeValid(obs[0].begin, obs[0].end))
+    {
+        cout << "First source [" << obs[0].name << "] removed" << endl;
+
+        obs[1].begin = obs[0].begin;
+        obs.erase(obs.begin());
+
+        return true;
+    }
+
+    // Try to remove the second source, check if the first source fullfills all limits
+    if (obs[0].IsRangeValid(obs[1].begin, obs[1].end))
+    {
+        cout << "Second source [" << obs[1].name << "] removed" << endl;
+
+        obs[0].end = obs[1].end;
+        obs.erase(obs.begin()+1);
+
+        if (obs.size()==0 || obs[0].name!=obs[1].name)
+            return true;
+
+        obs[0].end = obs[1].end;
+        obs.erase(obs.begin()+1);
+
+        cout << "Combined first two indentical sources [" << obs[0].name << "] into one observation" << endl;
+
+        return true;
+    }
+
+    cout << "No reschedule possible within limit." << endl;
+
+    return false;
+}
+
+bool RescheduleLastSources(vector<Source> &obs)
+{
+    // If observation time is smaller than 40min for the first source
+    // extend it to 40min if zenith angle will not go above 60deg.
+    const int last = obs.size()-1;
+    if (obs.size()<2 || obs[last].duration()>=40./24/60 || obs[last].name=="SLEEP" || obs[last-1].name=="SLEEP")
+        return false;
+
+    cout << "Last source [" << obs[last].name << "] detected < 40min" << endl;
+
+    const double obs1_duration = obs[last].end - 40./24/60 - obs[last-1].begin;
+    const double obs0_begin    = obs[last].end - 40./24/60;
+
+    // Check that:
+    //  - the duration for the shrunken source obs[1] is still above 40min
+    //  - obs[0] does not exceed 60deg at the end of its new window
+    //  - obs[0] does not exceed any limit within the new window
+
+    if (obs1_duration>=40./24/60 && obs[last].IsRangeValid(obs0_begin, obs[last].begin))
+    {
+        obs[last].begin = obs0_begin;
+        obs[last-1].end = obs0_begin;
+
+        cout << "Last source [" << obs[last].name << "] extended to 40min" << endl;
+
+        return false;
+    }
+
+    // Try to remove the last source, check if second source fullfills all limits
+    if (obs[last-1].IsRangeValid(obs[last].begin, obs[last].end))
+    {
+        cout << "Last source [" << obs[last].name << "] removed" << endl;
+
+        obs[last-1].end = obs[last].end;
+        obs.pop_back();
+
+        return true;
+    }
+
+    // Try to remove the second last source, check if the first source fullfills all limits
+    if (obs[last].IsRangeValid(obs[last-1].begin, obs[last-1].end))
+    {
+        cout << "Second last source [" << obs[last-1].name << "] removed" << endl;
+
+        obs[last].begin = obs[last-1].begin;
+        obs.erase(obs.begin()+obs.size()-2);
+
+        if (obs.size()==0 || obs[last-1].name!=obs[last-2].name)
+            return true;
+
+        obs[last-2].end = obs[last-1].end;
+        obs.pop_back();
+
+        cout << "Combined last two indentical sources [" << obs[last-1].name << "] into one observation" << endl;
+
+        return true;
+    }
+
+    cout << "No reschedule possible within limit." << endl;
+
+    return false;
+}
+
+bool RescheduleIntermediateSources(vector<Source> &obs)
+{
+    for (size_t i=1; i<obs.size()-1; i++)
+    {
+        if (obs[i].duration()>=40./24/60)
+            continue;
+
+        if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
+            continue;
+
+        cout << "Intermediate source [" << obs[i].name << "] detected < 40min" << endl;
+
+        double intersection = -1;
+
+        if (obs[i-1].name=="SLEEP")
+            intersection = obs[i].begin;
+
+        if (obs[i+1].name=="SLEEP")
+            intersection = obs[i].end;
+
+        if (obs[i-1].name==obs[i+1].name)
+            intersection = obs[i].begin;
+
+        if (intersection<0)
+        {
+            const double step = 1./24/60;
+            for (double jd=obs[i].begin; jd<obs[i].end-step/2; jd+=step)
+            {
+                const SolarObjects so(jd);
+                if (obs[i-1].getThreshold(so)>=obs[i+1].getThreshold(so))
+                {
+                    intersection = jd;
+                    break;
+                }
+            }
+        }
+
+        if ((obs[i-1].name!="SLEEP" && !obs[i-1].IsRangeValid(obs[i-1].end, intersection)) ||
+            (obs[i+1].name!="SLEEP" && !obs[i+1].IsRangeValid(intersection, obs[i+1].begin)))
+        {
+            cout << "No reschedule possible within limits." << endl;
+            continue;
+        }
+
+        cout << "Intermediate source [" << obs[i].name << "] removed" << endl;
+
+        const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;
+
+        obs[i-1].end   = intersection;
+        obs[i+1].begin = intersection;
+        obs.erase(obs.begin()+i);
+
+        i--;
+
+        if (obs.size()>1 && obs[i].name==obs[i+1].name)
+        {
+            obs[i].end = obs[i+1].end;
+            obs.erase(obs.begin()+i+1);
+
+            cout << "Combined two surrounding indentical sources [" << obs[i].name << "] into one observation" << endl;
+
+            i--;
+
+            continue;
+        }
+
+        if (underflow)
+            cout << "WARNING - Neighbor source < 40min as well." << endl;
+    }
+    return false;
+}
+
+void Print(const vector<Source> &obs, double startup_offset)
+{
+    cout << Time(obs[0].begin-startup_offset).GetAsStr() << "  STARTUP\n";
+    for (const auto& src: obs)
+    {
+        string tm = Time(src.begin).GetAsStr();
+
+        if (src.preobs.size()>0)
+        {
+            for (const auto& pre: src.preobs)
+            {
+                cout << tm << "  " << pre << "\n";
+                tm = "                   ";
+            }
+        }
+
+        cout << tm << "  " << src.name << " [";
+        cout << src.duration()*24*60 << "'";
+        if (src.name!="SLEEP")
+            cout << Tools::Form("; %.1f/%.1f", src.zd(src.begin), src.zd(src.end));
+        cout << "]";
+
+        if (src.duration()*24*60<40)
+            cout << " (!)";
+
+        cout << "\n";
+    }
+    cout << Time(obs.back().end).GetAsStr() << "  SHUTDOWN" << endl;
+}
+
+int main(int argc, const char* argv[])
+{
+//    gROOT->SetBatch();
+
+    Configuration conf(argv[0]);
+    conf.SetPrintUsage(PrintUsage);
+    SetupConfiguration(conf);
+
+    if (!conf.DoParse(argc, argv, PrintHelp))
+        return 127;
+
+    // ------------------ Eval config ---------------------
+
+    Time time;
+    if (conf.Has("date"))
+        time.SetFromStr(conf.Get<string>("date")+" 12:00:00");
+
+    Source::max_current = conf.Get<double>("max-current");
+    Source::max_zd      = conf.Get<double>("max-zd");
+
+    const double startup_offset = conf.Get<double>("startup.offset")/60/24;
+
+    const double angle_sun_set  = conf.Get<double>("data-taking.start");
+    const double angle_sun_rise = conf.Get<double>("data-taking.end");
+
+    if (startup_offset<0 || startup_offset>120)
+        throw runtime_error("Only values [0;120] are allowed for startup.offset");
+
+    if (angle_sun_set>-6)
+        throw runtime_error("datataking.start not allowed before sun at -6deg");
+
+    if (angle_sun_rise>-6)
+        throw runtime_error("datataking.end not allowed after sun at -6deg");
+
+    // -12: nautical
+    // Sun set with the same date than th provided date
+    // Sun rise on the following day
+    const RstTime sun_set  = GetSolarRst(time.JD()-0.5, angle_sun_set);
+    const RstTime sun_rise = GetSolarRst(time.JD()+0.5, angle_sun_rise);
+
+    const double sunset  = ceil(sun_set.set*24*60)/24/60;
+    const double sunrise = floor(sun_rise.rise*24*60)/24/60;
+
+    cout << "\n";
+    cout << "Date: " << Time(floor(sunset)) << "\n";
+    cout << "Set:  " << Time(sun_set.set)   << "\n";
+    cout << "Rise: " << Time(sun_rise.rise) << "\n";
+    cout << endl;
+
+    // ------------- Get Sources from databasse ---------------------
+
+    const vector<string> ourcenames = conf.Vec<string>("source");
+    const vector<string> sourcenames = conf.Vec<string>("source");
+    cout << "Nsources = " << sourcenames.size() << "\n";
+
+    string query = "SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1";
+    if (sourcenames.size()>0)
+        query += " AND fSourceName IN ('" + boost::algorithm::join(sourcenames, "', '")+"')";
+
+    const string fDatabase = conf.Get<string>("source-database");
+    const mysqlpp::StoreQueryResult res =
+        Database(fDatabase).query(query).store();
+
+    // ------------------ Eval config ---------------------
+
+    vector<Source> sources;
+    for (const auto &row: res)
+    {
+        const string name = string(row[0]);
+
+        Source src(name);
+
+        src.equ.ra  = double(row[1])*15;
+        src.equ.dec = double(row[2]);
+
+        src.maxzd = MyDouble(conf, "setup."+name+".max-zd");
+        src.maxcurrent = MyDouble(conf, "setup."+name+".max-current");
+        src.penalty = conf.Has("setup."+name+".penalty") ?
+            conf.Get<double>("setup."+name+".penalty") : 1;
+
+        src.preobs = conf.Vec<string>("preobs."+name);
+
+
+        cout << "[" << name << "]";
+
+        if (src.maxzd.valid)
+            cout << " Zd<" << src.maxzd.val;
+        if (src.penalty!=1)
+            cout << " Penalty=" << src.penalty;
+
+        cout << " " << boost::algorithm::join(src.preobs, "+") << endl;
+
+        /*
+         RstTime t1 = GetObjectRst(floor(sunset)-1, src.equ);
+         RstTime t2 = GetObjectRst(floor(sunset),   src.equ);
+
+         src.rst.transit = t1.transit<floor(sunset) ? t2.transit : t1.transit;
+         src.rst.rise = t1.rise>src.rst.transit ? t2.rise : t1.rise;
+         src.rst.set  = t1.set <src.rst.transit ? t2.set  : t1.set;
+         */
+
+        sources.emplace_back(src);
+    }
+    cout << endl;
+
+    // -------------------------------------------------------------------------
+
+    vector<Source> obs;
+
+    const double step = 1./24/60;
+    for (double jd=sunset; jd<sunrise-step/2; jd+=step)
+    {
+        const SolarObjects so(jd);
+
+        vector<Source> vis;
+        for (auto& src: sources)
+        {
+            if (src.valid(so))
+                vis.emplace_back(src);
+        }
+
+        // In case no source was found, add a sleep source
+        Source src("SLEEP");
+        vis.emplace_back(src);
+
+        // Source has higher priority if minimum observation time not yet fullfilled
+        sort(vis.begin(), vis.end(), SortByThreshold);
+
+        if (obs.size()>0 && obs.back().name==vis[0].name)
+            continue;
+
+        vis[0].begin = jd;
+        obs.emplace_back(vis[0]);
+    }
+
+    if (obs.size()==0)
+    {
+        cout << "No source found." << endl;
+        return 0;
+    }
+
+    // -------------------------------------------------------------------------
+
+    for (auto it=obs.begin(); it<obs.end()-1; it++)
+        it[0].end = it[1].begin;
+    obs.back().end = sunrise;
+
+    // -------------------------------------------------------------------------
+
+    Print(obs, startup_offset);
+    cout << endl;
+
+    // -------------------------------------------------------------------------
+
+    while (RescheduleFirstSources(obs));
+    while (RescheduleLastSources(obs));
+    while (RescheduleIntermediateSources(obs));
+
+    // ---------------------------------------------------------------------
+
+    Print(obs, startup_offset);
+    cout << endl;
+
+    // ---------------------------------------------------------------------
+
+    return 1;
+
+    // ------------- Create canvases and frames ---------------------
+/*
+    // It is important to use an offset which is larger than
+    // 1970-01-01 00:00:00. This one will not work if your
+    // local time zone is positive!
+    TH1S hframe("", "", 1, Time(sunset).Mjd()*24*3600, Time(sunrise).Mjd()*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.GetYaxis()->SetMoreLogLabels();
+    hframe.SetMinimum(0.9);
+    hframe.SetMaximum(11);
+    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 ---------------------
+
+    for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
+    {
+        // Eval row
+        const string name = (*v)[0].c_str();
+
+        EquPosn pos;
+        pos.ra  = double((*v)[1])*15;
+        pos.dec = double((*v)[2]);
+
+        // Loop over 24 hours
+        for (double h=0; h<1; h+=1./(24*12))
+        {
+            const SolarObjects so(jd+h);
+
+            // get local position of source
+            const HrzPosn hrz = GetHrzFromEqu(pos, so.fJD);
+
+            if (v==res.begin())
+                cout << Time(so.fJD) <<" " << 90-so.fMoonHrz.alt <<  endl;
+
+            const double cur = FACT::PredictI(so, pos);
+
+            // Relative  energy threshold prediction
+            const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
+
+            // Add points to curve
+            const double axis = Time(so.fJD).Mjd()*24*3600;
+
+            // If there is a gap of more than one bin, start a new curve
+
+            // 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*pow(cur/6.2, 0.394));
+
+            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
+            {
+                const double angle = GetAngularSeparation(so.fMoonEqu, pos);
+                g4.SetPoint(g4.GetN(), axis, angle);
+            }
+
+            if (cnt==0)
+                gm.SetPoint(gm.GetN(), axis, so.fMoonHrz.alt);
+        }
+    }
+
+    // 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());
+
+    c1.Print((t+".pdf(").c_str(), "pdf");
+    c2.Print((t+".pdf" ).c_str(), "pdf");
+    c3.Print((t+".pdf" ).c_str(), "pdf");
+    c4.Print((t+".pdf" ).c_str(), "pdf");
+    c5.Print((t+".pdf)").c_str(), "pdf");
+*/
+}