#include "externals/Prediction.h" //#include #include #include "Database.h" #include "tools.h" #include "Time.h" #include "Configuration.h" /* #include #include #include #include #include */ 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(), "SQL time (UTC), e.g. '2016-12-24'") ("source-database", var(""), "Database link as in\n\tuser:password@server[:port]/database.") ("max-current", var(90), "Global maximum current limit in uA") ("max-zd", var(75), "Global zenith distance limit in degree") ("source", vars(), "List of all TeV sources to be included, names according to the database") ("setup.*", var(), "Setup for the sources to be observed") ("preobs.*", vars(), "Prescheduled observations") ("startup.offset", var(15), "Determines how many minutes the startup is scheduled before data-taking.start [0;120]") ("data-taking.start", var(-12), "Begin of data-taking in degree of sun below horizon") ("data-taking.end", var(-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(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 preobs; Source(const string &n="") : name(n), begin(0), threshold(std::numeric_limits::max()) { } //bool IsSpecial() const { return threshold==std::numeric_limits::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; jdmax_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; const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664); threshold = penalty*ratio*pow(current/6.2, 0.394); return true; } }; double Source::max_zd; double Source::max_current; bool SortByThreshold(const Source &i, const Source &j) { return i.threshold &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 &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 &obs) { for (size_t i=1; i=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+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 &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("date")+" 12:00:00"); Source::max_current = conf.Get("max-current"); Source::max_zd = conf.Get("max-zd"); const double startup_offset = conf.Get("startup.offset")/60/24; const double angle_sun_set = conf.Get("data-taking.start"); const double angle_sun_rise = conf.Get("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 ourcenames = conf.Vec("source"); const vector sourcenames = conf.Vec("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("source-database"); const mysqlpp::StoreQueryResult res = Database(fDatabase).query(query).store(); // ------------------ Eval config --------------------- vector 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("setup."+name+".penalty") : 1; src.preobs = conf.Vec("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.transitsrc.rst.transit ? t2.rise : t1.rise; src.rst.set = t1.set obs; const double step = 1./24/60; for (double jd=sunset; jd vis; for (auto& src: sources) { if (src.calcThreshold(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(); itSetTimeFormat("%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::const_iterator v=res.begin(); v