Changeset 19684 for trunk/FACT++
- Timestamp:
- 09/27/19 16:27:20 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/FACT++/src/makeschedule.cc
r19139 r19684 19 19 control.add_options() 20 20 ("date", var<string>(), "SQL time (UTC), e.g. '2016-12-24' (equiv. '2016-12-24 12:00:00' to '2016-12-25 11:59:59')") 21 ("svg", var<string>(), "File name to write a SVG file. Use together with --loop") 22 ("hist", var<string>(), "File name to write a histogram file. Use together with --loop; default: false") 23 ("loop", var<uint16_t>(uint16_t(0)), "Number of days to loop over (0:default, >=1:turns off database entry)") 21 24 ("source-database", var<string>()->required(), "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].") 22 25 ("schedule-database", var<string>(), "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].") 23 26 ("max-current", var<double>(90), "Global maximum current limit in uA") 24 27 ("max-zd", var<double>(75), "Global zenith distance limit in degree") 28 ("min-moon", var<double>(5), "Minimum angular separation to moon in degree") 29 ("exponent-visibility-ratio", var<double>(1), "Global exponent for visibility ratio") 30 ("exponent-time-shift", var<double>(1), "Global exponent for time shift visibility ratio") 25 31 ("source", vars<string>(), "List of all TeV sources to be included, names according to the database") 26 32 ("setup.*", var<double>(), "Setup for the sources to be observed") … … 29 35 ("data-taking.start", var<double>(-12), "Begin of data-taking in degree of sun below horizon") 30 36 ("data-taking.end", var<double>(-13.75), "End of data-taking in degree of sun below horizon") 37 ("use-visibility-ratio", po_bool(), "Use visibility ratio for scaling of threshold)") 38 ("use-lst-shadow", po_bool(), "Respect shadow from LST for scheduling") 39 ("print-hist", var<bool>(), "Print histogram to console; default: false)") 31 40 ("enter-schedule-into-database", var<bool>(), "Enter schedule into database (required schedule-database, false: dry-run)") 32 41 ; … … 34 43 po::positional_options_description p; 35 44 p.add("date", 1); // The first positional options 45 p.add("loop", 1); // The first positional options 36 46 37 47 conf.AddOptions(control); … … 44 54 "makeschedule - Creates an automatic schedule\n" 45 55 "\n" 46 "Usage: makeschedule [yyyy-mm-dd ]\n";56 "Usage: makeschedule [yyyy-mm-dd [N]]\n"; 47 57 cout << endl; 48 58 } … … 50 60 void PrintHelp() 51 61 { 52 #ifdef HAVE_ROOT53 62 cout << 54 63 "First for each minute of the night a list is calculated of all " 55 "selected sources fulfil ing all global and all source specific "64 "selected sources fulfilling all global and all source specific " 56 65 "constraints, e.g. on the zenith distance or the current.\n" 57 66 "\n" … … 105 114 "\n"; 106 115 cout << endl; 107 #endif108 116 } 109 117 … … 141 149 static double max_current; 142 150 static double max_zd; 143 144 // Source descrition 151 static double min_moon; 152 static double vis_exponent; 153 static double ts_exponent; 154 static bool vis_ratio; 155 static bool lst_shadow; 156 157 // Source description 145 158 string name; 146 159 uint16_t key; … … 151 164 MyDouble maxcurrent; 152 165 double penalty; 166 double time_shift; 167 double visexponent; 168 double tsexponent; 169 170 // Derived values 171 double visratio; 172 double timeshift; 153 173 154 174 // Possible observation time … … 168 188 //bool IsSpecial() const { return threshold==std::numeric_limits<double>::max(); } 169 189 190 //distribution of zenith distance 191 int zddistr[90]; 192 170 193 double zd(const double &jd) const 171 194 { … … 173 196 } 174 197 198 /* 199 // See https://www.fact-project.org/logbook/showthread.php?tid=6601 200 201 75 -135 outmost knot left side 202 70 -134 tower left side 203 65 -128 dish left side 204 58 -125 intermediate position of dish on left side 205 52 -118 intermediate position of dish on left side 206 50 ~ -110 top of dish - to have the rail outside the FoV zd has to be 49 (see attached picture taken at zd50, az-113.5) 207 52 -100 intermediate position of dish on right side 208 58 -92 intermediate position of dish on right side 209 65 -89 dish right side 210 71 -85 tower right side 211 75 -84 outmost knot right side 212 */ 213 214 bool HasShadowFromLST(const ZdAzPosn &hrz) const 215 { 216 return (hrz.zd > 0.0379373*pow(hrz.az+109.23, 2) + 49.1224) 217 || (hrz.az>-79 && hrz.az<=-59 && hrz.zd>=70); 218 } 219 175 220 bool valid(const SolarObjects &so) const 176 221 { 177 const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD); 178 const double current = FACT::PredictI(so, equ); 222 const ZdAzPosn hrz = GetHrzFromEqu(equ, so.fJD+timeshift); 223 const ZdAzPosn true_hrz = timeshift==0 ? hrz : ZdAzPosn(GetHrzFromEqu(equ, so.fJD)); 224 const double current = FACT::PredictI(so, equ); 179 225 180 226 if (current>max_current) 181 227 return false; 182 228 183 if ( hrz.alt<=0 || 90-hrz.alt>max_zd)229 if (true_hrz.zd>max_zd) 184 230 return false; 185 231 186 if (maxzd.valid && 90-hrz.alt>maxzd.val)232 if (maxzd.valid && true_hrz.zd>maxzd.val) 187 233 return false; 188 234 235 if (lst_shadow && HasShadowFromLST(true_hrz)) 236 return false; 237 189 238 if (maxcurrent.valid && current>maxcurrent.val) 239 return false; 240 241 if (Nova::GetAngularSeparation(so.fMoonEqu, equ)<min_moon) 190 242 return false; 191 243 … … 205 257 double getThreshold(const SolarObjects &so) const 206 258 { 207 const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD); 208 const double current = FACT::PredictI(so, equ); 209 210 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664); 211 212 return penalty*ratio*pow(current/6.2, 0.394); 259 const ZdAzPosn hrz = GetHrzFromEqu(equ, so.fJD+timeshift); 260 const double current = FACT::PredictI(so, equ); 261 262 const double ratio = pow(cos(hrz.zd*M_PI/180), -2.664); 263 const double ratio2 = pow(current/6.2, 0.394); 264 265 return penalty*visratio*ratio*ratio2; 213 266 } 214 267 215 268 bool calcThreshold(const SolarObjects &so) 216 269 { 217 const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD); 218 const double current = FACT::PredictI(so, equ); 270 const ZdAzPosn hrz = GetHrzFromEqu(equ, so.fJD+timeshift); 271 const ZdAzPosn true_hrz = timeshift==0 ? hrz : ZdAzPosn(GetHrzFromEqu(equ, so.fJD)); 272 const double current = FACT::PredictI(so, equ); 219 273 220 274 if (current>max_current) 221 275 return false; 222 276 223 if ( hrz.alt<=0 || 90-hrz.alt>max_zd)277 if (true_hrz.zd>max_zd) 224 278 return false; 225 279 226 if (maxzd.valid && 90-hrz.alt>maxzd.val) 280 if (maxzd.valid && true_hrz.zd>maxzd.val) 281 return false; 282 283 if (lst_shadow && HasShadowFromLST(true_hrz)) 227 284 return false; 228 285 … … 230 287 return false; 231 288 232 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664); 233 threshold = penalty*ratio*pow(current/6.2, 0.394); 289 if (Nova::GetAngularSeparation(so.fMoonEqu, equ)<min_moon) 290 return false; 291 292 const double ratio = pow(cos(hrz.zd*M_PI/180), -2.664); 293 const double ratio2 = pow(current/6.2, 0.394); 294 295 threshold = penalty*visratio*ratio*ratio2; 234 296 235 297 return true; … … 239 301 double Source::max_zd; 240 302 double Source::max_current; 303 double Source::min_moon; 304 double Source::vis_exponent; 305 double Source::ts_exponent; 306 bool Source::vis_ratio; 307 bool Source::lst_shadow; 241 308 242 309 bool SortByThreshold(const Source &i, const Source &j) { return i.threshold<j.threshold; } … … 367 434 bool RescheduleIntermediateSources(vector<Source> &obs) 368 435 { 436 bool changed = false; 437 369 438 for (size_t i=1; i<obs.size()-1; i++) 370 439 { … … 413 482 cout << "Intermediate source [" << obs[i].name << "] removed" << endl; 414 483 415 const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;484 // const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40; 416 485 417 486 obs[i-1].end = intersection; … … 421 490 i--; 422 491 492 changed = true; 493 423 494 if (obs.size()>1 && obs[i].name==obs[i+1].name) 424 495 { … … 433 504 } 434 505 435 if (underflow) 436 cout << "WARNING - Neighbor source < 40min as well." << endl; 437 } 438 return false; 439 } 440 441 void RemoveMiniSources(vector<Source> &obs) 442 { 506 // if (underflow) 507 // cout << "WARNING - Neighbor source < 40min as well." << endl; 508 } 509 return changed; 510 } 511 512 // Remove mini source is able to violate limits 513 bool RemoveMiniSources(vector<Source> &obs) 514 { 515 bool changed = false; 516 443 517 for (size_t i=1; i<obs.size()-1; i++) 444 518 { … … 446 520 continue; 447 521 448 if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP") 522 // if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP") 523 // continue; 524 525 cout << "Mini source [" << obs[i].name << "] detected < 5min" << endl; 526 527 if (obs[i-1].name==obs[i+1].name) 528 { 529 cout << "Combined surrounding identical sources [" << obs[i-1].name << "]" << endl; 530 531 obs[i-1].end = obs[i+1].end; 532 obs.erase(obs.begin()+i, obs.begin()+i+1); 533 i -= 2; 534 535 changed = true; 536 449 537 continue; 450 451 cout << "Mini source [" << obs[i].name << "] detected < 5min" << endl; 452 538 } 539 540 /* 453 541 if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP") 454 542 { … … 460 548 i -= 2; 461 549 550 changed = true; 551 462 552 cout << "Combined two surrounding sleep into one" << endl; 463 553 464 554 continue; 465 } 555 }*/ 556 557 // Note that this (intentionally) violates our safety limits 466 558 467 559 if (obs[i-1].name=="SLEEP") 468 560 { 561 cout << "Preceding SLEEP detected... extended next source [" << obs[i+1].name << "]" << endl; 562 469 563 obs[i-1].end = obs[i+1].begin; 470 564 obs.erase(obs.begin()+i); 471 565 i--; 472 566 473 c out << "Extended previous sleep" << endl;567 changed = true; 474 568 475 569 continue; … … 478 572 if (obs[i+1].name=="SLEEP") 479 573 { 480 obs[i+1].begin = obs[i-1].end; 574 cout << "Following SLEEP detected... extended previous source [" << obs[i-1].name << "]" << endl; 575 576 obs[i-1].end = obs[i+1].begin; 481 577 obs.erase(obs.begin()+i); 482 483 cout << "Extended following sleep" << endl;484 485 578 i--; 579 580 changed = true; 581 486 582 continue; 487 583 } 488 } 584 585 // Algorithm based on RescheduleIntermediateSources 586 // Note that this (intentionally) violates our safety limits 587 588 const SolarObjects so1(obs[i-1].end); 589 const SolarObjects so2(obs[i+1].begin); 590 591 const double th1 = obs[i-1].getThreshold(so1); 592 const double th2 = obs[i+1].getThreshold(so2); 593 594 const double intersection = th1<th2 ? obs[i+1].begin : obs[i-1].end; 595 596 cout << "Intermediate mini source [" << obs[i].name << "] detected" << endl; 597 cout << (th1<th2?"Previous":"Next") << " source [" << (th1<th2?obs[i-1].name:obs[i+1].name) << "] extended" << endl; 598 599 obs[i-1].end = intersection; 600 obs[i+1].begin = intersection; 601 obs.erase(obs.begin()+i); 602 603 i--; 604 605 changed = true; 606 } 607 608 return changed; 489 609 } 490 610 … … 504 624 } 505 625 506 void Print(c onst vector<Source> &obs, double startup_offset)507 { 508 cout << Time(obs[0].begin-startup_offset).GetAsStr() << " STARTUP\n";626 void Print(char id, const vector<Source> &obs, double startup_offset) 627 { 628 cout << id << " " << Time(obs[0].begin-startup_offset).GetAsStr() << " STARTUP\n"; 509 629 for (const auto& src: obs) 510 630 { … … 514 634 for (const auto& pre: src.preobs) 515 635 { 516 cout << tm << " " << pre << "\n";636 cout << id << " " << tm << " " << pre << "\n"; 517 637 tm = " "; 518 638 } 519 639 } 520 640 521 cout << tm << " " << src.name << " [";641 cout << id << " " << tm << " " << src.name << " ["; 522 642 cout << src.duration()*24*60 << "'"; 523 643 if (src.name!="SLEEP") 524 cout << Tools::Form("; %.1f/%.1f ", src.zd(src.begin), src.zd(src.end));644 cout << Tools::Form("; %.1f/%.1f; %.2f; %.2fh", src.zd(src.begin), src.zd(src.end), src.visratio, src.timeshift*24); 525 645 cout << "]"; 526 646 … … 530 650 cout << "\n"; 531 651 } 532 cout << Time(obs.back().end).GetAsStr() << " SHUTDOWN" << endl; 652 cout << id << " " << Time(obs.back().end).GetAsStr() << " SHUTDOWN" << endl; 653 } 654 655 void PrintRect(ofstream &svg, const uint16_t &x, const int16_t &y, const uint16_t &w, const uint16_t &h, const Source &src) 656 { 657 static vector<string> colors = 658 { 659 "gray", // rgb(0,0,0) 660 "red", 661 "olive", 662 "blue", 663 "maroon", 664 "green", 665 "teal", 666 "navy", 667 "purple", 668 "silver", 669 "lime", 670 "fuchsia", 671 "yellow", 672 }; 673 674 static map<uint16_t, uint16_t> col_keys; 675 676 const auto ic = col_keys.find(src.key); 677 if (ic==col_keys.end()) 678 { 679 if (col_keys.size()>=colors.size()) 680 cerr << "WARNING: More different sources (" << col_keys.size() << ") than colors (" << colors.size() << ")." << endl; 681 682 const auto n = col_keys.size(); 683 684 col_keys[src.key] = n%colors.size(); 685 686 svg << "<g>\n"; 687 svg << " <rect x=\"" << 0 << "\" y=\"" << n*18 << "\" width=\"140\" height=\"18\" fill=\"" << colors[col_keys[src.key]] << "\"/>\n"; 688 svg << " <text x=\"" << 0 << "\" y=\"" << n*18+14 << "\" font-family=\"Verdana\" font-size=\"16\" fill=\"white\">" << src.name << "</text>\n"; 689 svg << "</g>\n"; 690 } 691 692 const string color = colors[col_keys[src.key]]; 693 694 svg << "<g>\n"; 695 svg << " <defs>\n"; 696 svg << " <linearGradient x1=\"0\" x2=\"0\" y1=\"0\" y2=\"1\" id=\"grad\">\n"; 697 698 for (int i=0; i<11; i++) 699 { 700 const double opa = 1-src.zd(src.begin+(src.end-src.begin)*0.10*i)/(src.maxzd.valid?src.maxzd.val:src.max_zd); 701 svg << " <stop style=\"stop-color:" << color << "; stop-opacity:" << opa << "\" offset=\"" << i*10 << "%\"/>\n"; 702 } 703 svg << " </linearGradient>\n"; 704 svg << " </defs>\n"; 705 706 svg << " <rect"; 707 svg << " x=\"" << x*w+150 << "\""; 708 svg << " y=\"" << y << "\""; 709 svg << " width=\"" << w << "\""; 710 svg << " height=\"" << h << "\""; 711 svg << " style=\"fill:url(#grad);\""; 712 svg << "/>\n"; 713 svg << "</g>\n"; 714 } 715 716 717 void PrintSVG(ofstream &svg, const uint16_t &scale, const Time &T0, const vector<Source> &obs) 718 { 719 const double jd = floor(T0.JD()); 720 721 for (const auto& src: obs) 722 { 723 if (src.preobs.size()>0 || src.key==0 || src.name=="SLEEP") 724 continue; 725 726 const double diff = src.begin-jd; 727 728 const uint32_t b = 24*60*diff; 729 const uint32_t d = 24*60*(src.end-src.begin)+0.5; 730 731 PrintRect(svg, diff, b%(24*60)-6*60, scale, d, src); 732 } 533 733 } 534 734 … … 594 794 } 595 795 596 int main(int argc, const char* argv[]) 597 { 598 // gROOT->SetBatch(); 599 600 Configuration conf(argv[0]); 601 conf.SetPrintUsage(PrintUsage); 602 SetupConfiguration(conf); 603 604 if (!conf.DoParse(argc, argv, PrintHelp)) 605 return 127; 606 607 // ------------------ Eval config --------------------- 608 609 const int enter = conf.Has("enter-schedule-into-database") ? (conf.Get<bool>("enter-schedule-into-database") ? 1 : -1) : 0; 610 if (enter && !conf.Has("schedule-database")) 611 throw runtime_error("enter-schedule-into-database required schedule-database."); 612 613 Time time; 614 if (conf.Has("date")) 615 time.SetFromStr(conf.Get<string>("date")+" 12:00:00"); 616 617 if (enter && floor(time.JD())<ceil(Time().JD())) 618 throw runtime_error("Only future schedules can be entered into the database."); 619 620 Source::max_current = conf.Get<double>("max-current"); 621 Source::max_zd = conf.Get<double>("max-zd"); 622 623 const double startup_offset = conf.Get<double>("startup.offset")/60/24; 624 625 const double angle_sun_set = conf.Get<double>("data-taking.start"); 626 const double angle_sun_rise = conf.Get<double>("data-taking.end"); 627 628 if (startup_offset<0 || startup_offset>120) 629 throw runtime_error("Only values [0;120] are allowed for startup.offset"); 630 631 if (angle_sun_set>-6) 632 throw runtime_error("datataking.start not allowed before sun at -6deg"); 633 634 if (angle_sun_rise>-6) 635 throw runtime_error("datataking.end not allowed after sun at -6deg"); 636 637 // -12: nautical 638 // Sun set with the same date than th provided date 639 // Sun rise on the following day 640 const RstTime sun_set = GetSolarRst(floor(time.JD())-0.5, angle_sun_set); 641 const RstTime sun_rise = GetSolarRst(floor(time.JD())+0.5, angle_sun_rise); 642 643 const double sunset = ceil(sun_set.set*24*60) /24/60 + 1e-9; 644 const double sunrise = floor(sun_rise.rise*24*60)/24/60 + 1e-9; 645 646 cout << "\n"; 647 648 cout << "Date: " << Time(floor(sunset)).GetAsStr() << "\n"; 649 cout << "Set: " << Time(sunset).GetAsStr() << " [" << Time(sun_set.set) << "]\n"; 650 cout << "Rise: " << Time(sunrise).GetAsStr() << " [" << Time(sun_rise.rise) << "]\n"; 651 652 cout << "\n"; 653 654 cout << "Global maximum current: " << Source::max_current << " uA/pix\n"; 655 cout << "Global zenith distance: " << Source::max_zd << " deg\n"; 656 657 cout << "\n"; 658 659 // ------------- Get Sources from databasse --------------------- 660 661 const vector<string> ourcenames = conf.Vec<string>("source"); 796 vector<Source> GetSources(Configuration &conf) 797 { 662 798 const vector<string> sourcenames = conf.Vec<string>("source"); 663 799 cout << "Nsources = " << sourcenames.size() << "\n"; … … 676 812 for (const auto &row: res) 677 813 { 814 //cout << endl; 678 815 const string name = string(row[0]); 679 816 … … 687 824 src.penalty = conf.Has("setup."+name+".penalty") ? 688 825 conf.Get<double>("setup."+name+".penalty") : 1; 826 src.visexponent = conf.Has("setup."+name+".exponent-visibility-ratio") ? 827 conf.Get<double>("setup."+name+".exponent-visibility-ratio") : 1; 828 src.tsexponent = conf.Has("setup."+name+".exponent-time-shift") ? 829 conf.Get<double>("setup."+name+".exponent-time-shift") : 1; 830 src.time_shift = conf.Has("setup."+name+".time-shift") ? 831 conf.Get<double>("setup."+name+".time-shift")/24 : 0; 689 832 690 833 src.preobs = conf.Vec<string>("preobs."+name); 691 834 692 693 cout << "[" << name << "]"; 835 cout << "[" << src.key << ":" << name << "]"; 694 836 695 837 if (src.maxzd.valid) … … 697 839 if (src.penalty!=1) 698 840 cout << " Penalty=" << src.penalty; 841 if (src.visexponent!=1) 842 cout << " Exponent [VR]=" << src.visexponent; 843 if (src.tsexponent!=1) 844 cout << " Exponent [TS]=" << src.tsexponent; 845 if (src.time_shift!=0) 846 cout << " Time-shift=" << src.time_shift*24 << "h"; 699 847 700 848 cout << " " << boost::algorithm::join(src.preobs, "+") << endl; 701 849 850 for (int k=0; k<90; k++) 851 src.zddistr[k]=0; 702 852 /* 703 853 RstTime t1 = GetObjectRst(floor(sunset)-1, src.equ); … … 713 863 cout << endl; 714 864 715 // ------------------------------------------------------------------------- 716 717 vector<Source> obs; 718 719 const uint32_t n = nearbyint((sunrise-sunset)*24*60); 720 for (uint32_t i=0; i<n; i++) 721 { 722 const double jd = sunset + i/24./60.; 723 724 const SolarObjects so(jd); 725 726 vector<Source> vis; 865 return sources; 866 } 867 868 int main(int argc, const char* argv[]) 869 { 870 // gROOT->SetBatch(); 871 872 Configuration conf(argv[0]); 873 conf.SetPrintUsage(PrintUsage); 874 SetupConfiguration(conf); 875 876 if (!conf.DoParse(argc, argv, PrintHelp)) 877 return 127; 878 879 // ------------------ Eval config --------------------- 880 881 const int enter = conf.Has("enter-schedule-into-database") ? (conf.Get<bool>("enter-schedule-into-database") ? 1 : -1) : 0; 882 if (enter && !conf.Has("schedule-database")) 883 throw runtime_error("enter-schedule-into-database requires schedule-database."); 884 885 Time time; 886 if (conf.Has("date")) 887 time.SetFromStr(conf.Get<string>("date")+" 12:00:00"); 888 else 889 time=floor(Time().JD()); 890 891 /* 892 cout << "-------> " << setprecision(12) << time.JD() << endl; 893 cout << " " << time.GetAsStr() 894 << " " << ceil(Time().JD()) 895 << " " << floor(Time().JD()) 896 << endl; 897 */ 898 const uint16_t loop = conf.Get<uint16_t>("loop"); 899 900 if (enter && floor(time.JD())<ceil(Time().JD())) 901 throw runtime_error("Only future schedules can be entered into the database."); 902 903 Source::max_current = conf.Get<double>("max-current"); 904 Source::max_zd = conf.Get<double>("max-zd"); 905 Source::min_moon = conf.Get<double>("min-moon"); 906 Source::vis_exponent= conf.Get<double>("exponent-visibility-ratio"); 907 Source::ts_exponent = conf.Get<double>("exponent-time-shift"); 908 Source::vis_ratio = conf.Get<bool>("use-visibility-ratio"); 909 Source::lst_shadow = conf.Get<bool>("use-lst-shadow"); 910 911 const double startup_offset = conf.Get<double>("startup.offset")/60/24; 912 913 const double angle_sun_set = conf.Get<double>("data-taking.start"); 914 const double angle_sun_rise = conf.Get<double>("data-taking.end"); 915 916 if (startup_offset<0 || startup_offset>120) 917 throw runtime_error("Only values [0;120] are allowed for startup.offset"); 918 919 if (angle_sun_set>-6) 920 throw runtime_error("datataking.start not allowed before sun at -6deg"); 921 922 if (angle_sun_rise>-6) 923 throw runtime_error("datataking.end not allowed after sun at -6deg"); 924 925 cout << "\n"; 926 927 cout << "Global maximum current allowed: " << Source::max_current << " uA/pix\n"; 928 cout << "Global maximum zenith distance: " << Source::max_zd << " deg\n"; 929 cout << "Min. angular separation to moon: " << Source::min_moon << " deg\n"; 930 931 cout << "\n"; 932 933 // ------------- SVG --------------------------------------------- 934 935 const uint16_t svgscale = loop>0 ? ::min((2*1080)/loop, 5) : 0; 936 937 ofstream svg; 938 if (conf.Has("svg") && loop>0) 939 { 940 svg.open(conf.Get<string>("svg")); 941 if (!svg) 942 { 943 cerr << "ERROR: Writing to '" << conf.Get<string>("svg") << "' failed: " << strerror(errno) << endl; 944 return 1; 945 } 946 947 svg << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"; 948 svg << "<svg width=\"" << loop*svgscale+150 << "\" height=\"840\" xmlns=\"http://www.w3.org/2000/svg\">\n"; 949 svg << "<rect width=\"" << loop*svgscale+150 << "\" height=\"840\" fill=\"black\"/>\n"; 950 951 //svg << "<svg width=\"" << loop*svgscale+150 << "\" height=\"840\">\n"; 952 //svg << "<rect x=\"0\" y=\"0\" width=\"" << loop*svgscale+150 << "\" height=\"840\" fill=\"black\"/>\n"; 953 } 954 955 // ------------- Get Sources from databasse --------------------- 956 957 vector<Source> sources = GetSources(conf); 958 959 //cout << "VR: " << Source::vis_ratio << endl; 960 961 const uint16_t N = ::max<uint16_t>(loop, 1); 962 for (int iday=0; iday<N; iday++) 963 { 964 // ----------- Define time slots for the day --------------------- 965 966 // -12: nautical 967 // Sun set with the same date than th provided date 968 // Sun rise on the following day 969 const RstTime sun_set = GetSolarRst(floor(time.JD()+iday)-0.5, angle_sun_set); 970 const RstTime sun_rise = GetSolarRst(floor(time.JD()+iday)+0.5, angle_sun_rise); 971 972 const double sunset = ceil(sun_set.set*24*60) /24/60 + 1e-9; 973 const double sunrise = floor(sun_rise.rise*24*60)/24/60 + 1e-9; 974 975 cout << "Date: " << Time(floor(sunset)).GetAsStr() << "\n"; 976 cout << "Set: " << Time(sunset).GetAsStr() << " [" << Time(sun_set.set) << "]\n"; 977 cout << "Rise: " << Time(sunrise).GetAsStr() << " [" << Time(sun_rise.rise) << "]\n"; 978 cout << "\n"; 979 980 const double vis_sun = sun_rise.rise - sun_set.set; 981 727 982 for (auto& src: sources) 728 983 { 729 if (src.calcThreshold(so)) 730 vis.emplace_back(src); 984 const double maxzd = min(Source::max_zd, src.maxzd.valid?src.maxzd.val:90); 985 986 int8_t sign = 0; 987 988 RstTime rst; 989 const int rc = GetObjectRst(rst, src.equ, sun_set.set, 90-maxzd); 990 switch (rc) 991 { 992 case -1: 993 // circumpolar below horizon - not visible 994 src.visratio = 0; 995 break; 996 997 case 0: 998 { 999 // Although, according to the documentation, everything should be alinged 1000 // between JD and JD+1, it is not... therefore, we align ourselves 1001 if (rst.rise-sun_set.set>1) 1002 rst.rise -= 1; 1003 if (rst.set-sun_set.set>1) 1004 rst.set -= 1; 1005 if (rst.rise-sun_set.set<0) 1006 rst.rise += 1; 1007 if (rst.set-sun_set.set<0) 1008 rst.set += 1; 1009 1010 // Now: rst.rise/rst.set is aligned between sunset and sunset+1: 1011 // 1012 // JD JD+1 1013 // 1) SUN.SET | | sun.rise RST.RISE rst.set ... SUN.SET | 1014 // 2) SUN.SET |++++++++++++++++++++++++| sun.rise rst.set RST.RISE ... SUN.SET | 1015 // 3) SUN.SET | RST.RISE+++++++++| sun.rise rst.set ... SUN.SET | 1016 // 4) SUN.SET |++++++++rst.set | sun.rise rst.rise ... SUN.SET | 1017 // 5) SUN.SET | RST.RISE+++rst.set | sun.rise ... SUN.SET | 1018 // 6) SUN.SET |+++rst.set RST.RISE+++| sun.rise ... SUN.SET | 1019 1020 if (rst.rise<sun_rise.rise && sun_rise.rise<rst.set) 1021 sign = 1; 1022 if (rst.set<sun_rise.rise && sun_rise.rise<rst.rise) 1023 sign = -1; 1024 1025 // Identify case 6 1026 const bool case6 = rst.set<rst.rise && rst.rise<sun_rise.rise; 1027 1028 // make sure that the visibility calculation of the source 1029 // yields positive valiues and the overlap calculation works 1030 if (rst.set<rst.rise) 1031 rst.rise -= 1; 1032 1033 // Now: rst.rise is always before rst.set 1034 // 1035 // JD JD+1 1036 // 1) ... SUN.SET | | sun.rise RST.RISE rst.set ... 1037 // 2) ... RST.RISE SUN.SET |++++++++++++++++++++++++| sun.rise rst.set ... 1038 // 3) ... SUN.SET | RST.RISE+++++++++| sun.rise rst.set ... 1039 // 4) ... RST.RISE SUN.SET |++++++++rst.set | sun.rise rst.rise ... 1040 // 5) ... SUN.SET | RST.RISE+++rst.set | sun.rise ... 1041 // 6) ... RST.RISE SUN.SET |+++rst.set [RST.RISE]++| sun.rise ... 1042 1043 // Total visibility of the source itself (this is now always positive) 1044 const double vis_rst = rst.set-rst.rise; 1045 1046 // Calculate overlap region 1047 // Visibility of the source between sunset and sunrise 1048 const double vis_night = case6 ? vis_sun - (1-vis_rst) : min(rst.set, sun_rise.rise) - max(rst.rise, sun_set.set); 1049 1050 //vis_sun-(1-vis_rst) = vis_sun-(1-(set-rise)) = vis_sun-(1-(set-(rise-1))) = vis_sun-(1-set+(rise-1)) = vis_sun-(rise-set) 1051 1052 // General visibility of the source (but maximum one night) 1053 const double vis_total = min(vis_rst, vis_sun); 1054 1055 // Special treatment of case 1 (vis_night<0 means that there is no overlap) 1056 src.visratio = vis_night<0 ? 0 : vis_night / vis_total; 1057 1058 // cout << src.name << " " << Time(rst.rise).GetAsStr() << " " << Time(rst.set).GetAsStr() << " " << vis_night << " " << vis_total << " " << vis_rst << " " << vis_sun << " " << case6 << " " << test << endl; 1059 break; 1060 } 1061 1062 case 1: 1063 // circumpolar above horizon -> always visible 1064 src.visratio = 1; 1065 break; 1066 } 1067 1068 src.timeshift = src.time_shift * sign*(1-pow(src.visratio, src.tsexponent*Source::ts_exponent)); 1069 1070 src.visratio = Source::vis_ratio ? pow(src.visratio, src.visexponent*Source::vis_exponent) : 1; 731 1071 } 732 1072 733 // In case no source was found, add a sleep source 734 Source src("SLEEP"); 735 vis.emplace_back(src); 736 737 // Source has higher priority if minimum observation time not yet fullfilled 738 sort(vis.begin(), vis.end(), SortByThreshold); 739 740 if (obs.size()>0 && obs.back().name==vis[0].name) 1073 // ------------------------------------------------------------------------- 1074 1075 vector<Source> obs; 1076 1077 const uint32_t n = nearbyint((sunrise-sunset)*24*60); 1078 for (uint32_t i=0; i<n; i++) 1079 { 1080 const double jd = sunset + i/24./60.; 1081 1082 const SolarObjects so(jd); 1083 1084 vector<Source> vis; 1085 for (auto& src: sources) 1086 { 1087 if (src.calcThreshold(so)) 1088 vis.emplace_back(src); 1089 } 1090 1091 // In case no source was found, add a sleep source 1092 Source src("SLEEP"); 1093 vis.emplace_back(src); 1094 1095 // Source has higher priority if minimum observation time not yet fullfilled 1096 sort(vis.begin(), vis.end(), SortByThreshold); 1097 1098 if (obs.size()>0 && obs.back().name==vis[0].name) 1099 continue; 1100 1101 vis[0].begin = jd; 1102 obs.emplace_back(vis[0]); 1103 } 1104 1105 if (obs.size()==0) 1106 { 1107 cout << "No source found." << endl; 1108 return 1; 1109 } 1110 1111 // ------------------------------------------------------------------------- 1112 1113 for (auto it=obs.begin(); it<obs.end()-1; it++) 1114 it[0].end = it[1].begin; 1115 obs.back().end = sunrise; 1116 1117 // ------------------------------------------------------------------------- 1118 1119 Print('-', obs, startup_offset); 1120 cout << endl; 1121 1122 // ------------------------------------------------------------------------- 1123 1124 while (RescheduleFirstSources(obs)); 1125 while (RescheduleLastSources(obs)); 1126 while (RescheduleIntermediateSources(obs)); 1127 while (RemoveMiniSources(obs)); 1128 1129 CheckStartupAndShutdown(obs); 1130 1131 // --------------------------------------------------------------------- 1132 1133 cout << endl; 1134 Print('*', obs, startup_offset); 1135 cout << endl; 1136 1137 // --------------------------------------------------------------------- 1138 1139 //fill zd-distribution 1140 for (auto& src2: sources) 1141 { 1142 //cout << "source: " << src2.name << " " << src2.key << endl; 1143 for (const auto& src: obs) 1144 { 1145 //cout << "source: " << src.name << " " << src.key << endl; 1146 const uint32_t n2 = nearbyint((src.end-src.begin)*24*60); 1147 for (uint32_t i=0; i<n2; i++) 1148 { 1149 const double jd = src.begin + i/24./60.; 1150 const int zd = (int)src.zd(jd); 1151 //cout << " zd:" << floor(zd) << " " << zd << endl; 1152 //cout << " zd: " << zd << " " << src2.zddistr[zd] << flush; 1153 if (src.key==src2.key) 1154 { 1155 //cout << " -> " << flush; 1156 src2.zddistr[zd]++; 1157 } 1158 //cout << " " << src2.zddistr[zd] << endl; 1159 } 1160 } 1161 } 1162 1163 // --------------------------------------------------------------------- 1164 1165 if (svg.is_open()) 1166 PrintSVG(svg, svgscale, time, obs); 1167 1168 // --------------------------------------------------------------------- 1169 1170 if (loop>0) 741 1171 continue; 742 1172 743 vis[0].begin = jd; 744 obs.emplace_back(vis[0]); 745 } 746 747 if (obs.size()==0) 748 { 749 cout << "No source found." << endl; 1173 if (!enter) 1174 return 0; 1175 1176 const string scheduledb = conf.Get<string>("schedule-database"); 1177 1178 Database db(scheduledb); 1179 1180 if (enter>0) 1181 db.query("LOCK TABLES Schedule WRITE"); 1182 1183 const int rc = FillSql(db, enter, obs, startup_offset); 1184 1185 if (enter>0) 1186 db.query("UNLOCK TABLES"); 1187 1188 return rc; 1189 } 1190 1191 if (svg.is_open()) 1192 svg << "</svg>" << endl; 1193 1194 // --------------------------------------------------------------------- 1195 if (conf.Has("print-hist") && conf.Get<bool>("print-hist")) 1196 { 1197 for (const auto& src: sources) 1198 { 1199 cout << "ZD-DISTRIBUTION [" << src.name << "]: " << flush; 1200 //print zd-distr 1201 for (int k=0; k<90; k++) 1202 cout << src.zddistr[k] << " " << flush; 1203 cout << endl; 1204 //plot zd-distr 1205 for (int k=0; k<90; k++) 1206 { 1207 if (src.zddistr[k]==0) 1208 continue; 1209 printf("%02d", k); 1210 for (int k2=0; k2<(int)(src.zddistr[k]/loop); k2++) 1211 cout << "*" << flush; 1212 cout << endl; 1213 } 1214 } 1215 } 1216 1217 if (!conf.Has("hist")) 1218 return 0; 1219 1220 ofstream hist(conf.Get<string>("hist").c_str()); 1221 if (!hist) 1222 { 1223 cerr << "ERROR: Writing to '" << conf.Get<string>("hist") << "' failed: " << strerror(errno) << endl; 750 1224 return 1; 751 1225 } 752 1226 753 // ------------------------------------------------------------------------- 754 755 for (auto it=obs.begin(); it<obs.end()-1; it++) 756 it[0].end = it[1].begin; 757 obs.back().end = sunrise; 758 759 // ------------------------------------------------------------------------- 760 761 Print(obs, startup_offset); 762 cout << endl; 763 764 // ------------------------------------------------------------------------- 765 766 while (RescheduleFirstSources(obs)); 767 while (RescheduleLastSources(obs)); 768 while (RescheduleIntermediateSources(obs)); 769 770 RemoveMiniSources(obs); 771 CheckStartupAndShutdown(obs); 772 773 // --------------------------------------------------------------------- 774 775 cout << endl; 776 Print(obs, startup_offset); 777 cout << endl; 778 779 // --------------------------------------------------------------------- 780 781 if (!enter) 782 return 0; 783 784 const string scheduledb = conf.Get<string>("schedule-database"); 785 786 Database db(scheduledb); 787 788 if (enter>0) 789 db.query("LOCK TABLES Schedule WRITE"); 790 791 const int rc = FillSql(db, enter, obs, startup_offset); 792 793 if (enter>0) 794 db.query("UNLOCK TABLES"); 795 796 // --------------------------------------------------------------------- 797 798 return rc; 799 } 1227 hist << "# " << flush; 1228 for (const auto& src: sources) 1229 hist << src.name << " " << flush; 1230 hist << endl; 1231 1232 hist << "# " << flush; 1233 for (const auto& src: sources) 1234 hist << src.key << " " << flush; 1235 hist << endl; 1236 1237 for (int k=0; k<90; k++) 1238 { 1239 hist << k << " " << flush; 1240 for (const auto& src: sources) 1241 hist << src.zddistr[k] << " " << flush; 1242 hist << endl; 1243 } 1244 1245 return 0; 1246 }
Note:
See TracChangeset
for help on using the changeset viewer.