#include "Prediction.h" #include #include "Database.h" #include "tools.h" #include "Time.h" #include "Configuration.h" using namespace std; using namespace Nova; // ----------------------------------------------------------------------- void SetupConfiguration(Configuration &conf) { po::options_description control("Makeschedule"); control.add_options() ("date", var(), "SQL time (UTC), e.g. '2016-12-24' (equiv. '2016-12-24 12:00:00' to '2016-12-25 11:59:59')") ("source-database", var()->required(), "Database link as in\n\tuser:password@server[:port]/database.") ("schedule-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") ("enter-schedule-into-database", var(), "Enter schedule into database (required schedule-database, false: dry-run)") ; 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 << "First for each minute of the night a list is calculated of all " "selected sources fulfiling all global and all source specific " "constraints, e.g. on the zenith distance or the current.\n" "\n" "The remaining source list is sorted by the relative threshold, " "while the threshold is weighted with a user defined source " "specific penalty. The first source of the list is taken to " "be observed.\n" "\n" "In a next step the first and last source of the resulting schedule " "are evaluated. If their observation time is below 40', it is tried " "to extend it to 40min. If this violates one of the criteria mentioned " "above or gives an observation time for the neighbouring source of " "less than 40min, try to replace it by the neighbouring source. " "If this also does not fulfil the requirements, the original " "schedule remains unchanged.\n" "\n" "Now a similar check is run for all intermediate sources. They are " "checked (from the beginning to the end, one by one), if they have " "an observation time of less than 40min. In this case, it is tried " "to remove them. The observation of the two neighbouring sources is " "extended to their penalized point of equal relative threshold. " "If this solution would not fulfil all criteria, no change is made.\n" "\n" "In a last step, all remaining sources with less than 5min " "observation time are replaced with sleep and sleep after startup " "or before shutdown are removed.\n" "\n" "\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; uint16_t key; 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="", uint16_t k=-1) : name(n), key(k), 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 uint32_t n = nearbyint((jd_end-jd_begin)*24*60); for (uint32_t i=0; imax_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 uint32_t n = nearbyint((obs[i].end-obs[i].begin)*24*60); for (uint32_t ii=0; ii=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 RemoveMiniSources(vector &obs) { for (size_t i=1; i=5./24/60) continue; if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP") continue; cout << "Mini source [" << obs[i].name << "] detected < 5min" << endl; if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP") { obs[i-1].end = obs[i+2].begin; obs.erase(obs.begin()+i+1); obs.erase(obs.begin()+i); i -= 2; cout << "Combined two surrounding sleep into one" << endl; continue; } if (obs[i-1].name=="SLEEP") { obs[i-1].end = obs[i+1].begin; obs.erase(obs.begin()+i); i--; cout << "Extended previous sleep" << endl; continue; } if (obs[i+1].name=="SLEEP") { obs[i+1].begin = obs[i-1].end; obs.erase(obs.begin()+i); cout << "Extended following sleep" << endl; i--; continue; } } } void CheckStartupAndShutdown(vector &obs) { if (obs.front().name=="SLEEP") { obs.erase(obs.begin()); cout << "Detected sleep after startup... removed." << endl; } if (obs.back().name=="SLEEP") { obs.pop_back(); cout << "Detected sleep before shutdown... removed." << endl; } } 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 FillSql(Database &db, int enter, const vector &obs, double startup_offset) { const string query0 = "SELECT COUNT(*) FROM Schedule WHERE DATE(ADDTIME(fStart, '-12:00')) = '"+Time(obs[0].begin).GetAsStr("%Y-%m-%d")+"'"; const mysqlpp::StoreQueryResult res0 = db.query(query0).store(); if (res0.num_rows()!=1) { cout << "Check for schedule size failed." << endl; return 10; } if (uint32_t(res0[0][0])!=0) { cout << "Schedule not empty." << endl; return 11; } const mysqlpp::StoreQueryResult res1 = db.query("SELECT fMeasurementTypeName, fMeasurementTypeKEY FROM MeasurementType").store(); map types; for (const auto &row: res1) types.insert(make_pair(string(row[0]), uint32_t(row[1]))); ostringstream str; str << "INSERT INTO Schedule (fStart, fUser, fMeasurementID, fMeasurementTypeKEY, fSourceKEY) VALUES "; str << "('" << Time(obs[0].begin-startup_offset).GetAsStr() << "', 'auto', 0, " << types["Startup"] << ", NULL),\n"; // [Startup]\n"; for (const auto& src: obs) { string tm = Time(src.begin).GetAsStr(); /* if (src.preobs.size()>0) { for (const auto& pre: src.preobs) { str << tm << " " << pre << "\n"; tm = " "; } }*/ if (src.name!="SLEEP") str << "('" << tm << "', 'auto', 0, " << types["Data"] << ", " << src.key << "),\n"; // [Data: " << src.name << "]\n"; else str << "('" << tm << "', 'auto', 0, " << types["Sleep"] << ", NULL),\n"; // [Sleep]\n"; } str << "('" << Time(obs.back().end).GetAsStr() << "', 'auto', 0, " << types["Shutdown"] << ", NULL)";// [Shutdown]"; if (enter<0) { cout << str.str() << endl; return 0; } db.query(str.str()).exec(); cout << "Schedule entered successfully into database." << endl; return 0; } 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 --------------------- const int enter = conf.Has("enter-schedule-into-database") ? (conf.Get("enter-schedule-into-database") ? 1 : -1) : 0; if (enter && !conf.Has("schedule-database")) throw runtime_error("enter-schedule-into-database required schedule-database."); Time time; if (conf.Has("date")) time.SetFromStr(conf.Get("date")+" 12:00:00"); if (enter && floor(time.JD())("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(floor(time.JD())-0.5, angle_sun_set); const RstTime sun_rise = GetSolarRst(floor(time.JD())+0.5, angle_sun_rise); const double sunset = ceil(sun_set.set*24*60) /24/60 + 1e-9; const double sunrise = floor(sun_rise.rise*24*60)/24/60 + 1e-9; cout << "\n"; cout << "Date: " << Time(floor(sunset)).GetAsStr() << "\n"; cout << "Set: " << Time(sunset).GetAsStr() << " [" << Time(sun_set.set) << "]\n"; cout << "Rise: " << Time(sunrise).GetAsStr() << " [" << Time(sun_rise.rise) << "]\n"; cout << "\n"; cout << "Global maximum current: " << Source::max_current << " uA/pix\n"; cout << "Global zenith distance: " << Source::max_zd << " deg\n"; cout << "\n"; // ------------- 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, fSourceKEY, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1"; if (sourcenames.size()>0) query += " AND fSourceName IN ('" + boost::algorithm::join(sourcenames, "', '")+"')"; const string sourcedb = conf.Get("source-database"); const mysqlpp::StoreQueryResult res = Database(sourcedb).query(query).store(); // ------------------ Eval config --------------------- vector sources; for (const auto &row: res) { const string name = string(row[0]); Source src(name, row[1]); src.equ.ra = double(row[2])*15; src.equ.dec = double(row[3]); 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 uint32_t n = nearbyint((sunrise-sunset)*24*60); for (uint32_t i=0; i 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 1; } // ------------------------------------------------------------------------- for (auto it=obs.begin(); it("schedule-database"); Database db(scheduledb); if (enter>0) db.query("LOCK TABLES Schedule WRITE"); const int rc = FillSql(db, enter, obs, startup_offset); if (enter>0) db.query("UNLOCK TABLES"); // --------------------------------------------------------------------- return rc; }