source: trunk/FACT++/src/makeschedule.cc@ 18411

Last change on this file since 18411 was 18411, checked in by tbretz, 9 years ago
Added new program makeschedule.
File size: 22.6 KB
Line 
1#include "externals/Prediction.h"
2
3//#include <unordered_map>
4#include <boost/algorithm/string/join.hpp>
5
6#include "Database.h"
7
8#include "tools.h"
9#include "Time.h"
10#include "Configuration.h"
11
12/*
13#include <TROOT.h>
14#include <TH1.h>
15#include <TGraph.h>
16#include <TCanvas.h>
17#include <TLegend.h>
18*/
19
20using namespace std;
21using namespace Nova;
22
23// -----------------------------------------------------------------------
24
25/*
26void CheckForGap(TCanvas &c, TGraph &g, double axis)
27{
28 if (g.GetN()==0 || axis-g.GetX()[g.GetN()-1]<450)
29 return;
30
31 c.cd();
32 ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
33 while (g.GetN())
34 g.RemovePoint(0);
35}
36
37void DrawClone(TCanvas &c, TGraph &g)
38{
39 if (g.GetN()==0)
40 return;
41
42 c.cd();
43 ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
44}
45*/
46
47// ========================================================================
48// ========================================================================
49// ========================================================================
50
51void SetupConfiguration(Configuration &conf)
52{
53 po::options_description control("Makeschedule");
54 control.add_options()
55 ("date", var<string>(), "SQL time (UTC), e.g. '2016-12-24'")
56 ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
57 ("max-current", var<double>(90), "Global maximum current limit in uA")
58 ("max-zd", var<double>(75), "Global zenith distance limit in degree")
59 ("source", vars<string>(), "List of all TeV sources to be included, names according to the database")
60 ("setup.*", var<double>(), "Setup for the sources to be observed")
61 ("preobs.*", vars<string>(), "Prescheduled observations")
62 ("startup.offset", var<double>(15), "Determines how many minutes the startup is scheduled before data-taking.start [0;120]")
63 ("data-taking.start", var<double>(-12), "Begin of data-taking in degree of sun below horizon")
64 ("data-taking.end", var<double>(-13.75), "End of data-taking in degree of sun below horizon")
65 ;
66
67 po::positional_options_description p;
68 p.add("date", 1); // The first positional options
69
70 conf.AddOptions(control);
71 conf.SetArgumentPositions(p);
72}
73
74void PrintUsage()
75{
76 cout <<
77 "makeschedule - Creates an automatic schedule\n"
78 "\n"
79 "Usage: makeschedule [yyyy-mm-dd]\n";
80 cout << endl;
81}
82
83void PrintHelp()
84{
85#ifdef HAVE_ROOT
86 cout <<
87 "\n"
88 "Examples:\n"
89 "\n"
90 " makeschedule 2016-12-24\n"
91 "\n"
92 "Calculate the Christmas schedule for 2016 using all TeV sources from the\n"
93 "database. If the date is omitted the current date is used.\n"
94 "\n"
95 " makeschedule --source='Mrk 421' --source='Mrk 501' --source='Crab'\n"
96 "\n"
97 "Use only the mentioned sources to calculate the schedule.\n"
98 "\n"
99 " makeschedule --source=Crab --setup.Crab.max-zd=30\n"
100 "\n"
101 "Limit the zenith distance of Crab into the range [0;30]deg.\n"
102 "\n"
103 " makeschedule --source=Crab --setup.Crab.max-current=50\n"
104 "\n"
105 "Limit the maximum estimated current of Crab at 50uA.\n"
106 "\n"
107 " makeschedule --source='IC 310' '--setup.IC 310.penalty=1.2'\n"
108 "\n"
109 "Multiply IC310's estimated relative threshold by a factor 1.2\n"
110 "\n";
111 cout << endl;
112#endif
113}
114
115
116struct MyDouble
117{
118 double val;
119 bool valid;
120 MyDouble(Configuration &conf, const string &str) : val(0)
121 {
122 valid = conf.Has(str);
123 if (valid)
124 val = conf.Get<double>(str);
125 }
126 MyDouble() : val(0), valid(false) {}
127};
128
129/*
130struct MinMax
131{
132 MyDouble min;
133 MyDouble max;
134 MinMax(Configuration &conf, const string &str)
135 {
136 min = MyDouble(conf, str+".min");
137 max = MyDouble(conf, str+".max");
138 }
139 MinMax() {}
140};
141*/
142
143struct Source
144{
145 // Global limits
146 static double max_current;
147 static double max_zd;
148
149 // Source descrition
150 string name;
151 EquPosn equ;
152
153 // Source specific limits
154 MyDouble maxzd;
155 MyDouble maxcurrent;
156 double penalty;
157
158 // Possible observation time
159 double begin;
160 double end;
161
162 // Threshold (sorting reference)
163 double threshold;
164
165 double duration() const { return end-begin; };
166
167 // Pre-observations (e.g. ratescan)
168 vector<string> preobs;
169
170 Source(const string &n="") : name(n), begin(0), threshold(std::numeric_limits<double>::max()) { }
171
172 //bool IsSpecial() const { return threshold==std::numeric_limits<double>::max(); }
173
174 double zd(const double &jd) const
175 {
176 return 90-GetHrzFromEqu(equ, jd).alt;
177 }
178
179 bool valid(const SolarObjects &so) const
180 {
181 const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD);
182 const double current = FACT::PredictI(so, equ);
183
184 if (current>max_current)
185 return false;
186
187 if (hrz.alt<=0 || 90-hrz.alt>max_zd)
188 return false;
189
190 if (maxzd.valid && 90-hrz.alt>maxzd.val)
191 return false;
192
193 if (maxcurrent.valid && current>maxcurrent.val)
194 return false;
195
196 return true;
197 }
198
199 bool IsRangeValid(const double &jd_begin, const double &jd_end) const
200 {
201 const double step = 1./24/60;
202 for (double jd=jd_begin; jd<jd_end+step/2; jd++)
203 if (!valid(SolarObjects(jd)))
204 return false;
205
206 return true;
207 }
208
209 double getThreshold(const SolarObjects &so) const
210 {
211 const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD);
212 const double current = FACT::PredictI(so, equ);
213
214 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
215
216 return penalty*ratio*pow(current/6.2, 0.394);
217 }
218};
219
220double Source::max_zd;
221double Source::max_current;
222
223bool SortByThreshold(const Source &i, const Source &j) { return i.threshold<j.threshold; }
224
225bool RescheduleFirstSources(vector<Source> &obs)
226{
227 if (obs.size()<2 || obs[0].duration()>=40./24/60 || obs[0].name=="SLEEP" || obs[1].name=="SLEEP")
228 return false;
229
230 cout << "First source [" << obs[0].name << "] detected < 40min" << endl;
231
232 const double obs1_duration = obs[1].end - obs[0].begin - 40./24/60;
233 const double obs0_end = obs[0].begin + 40./24/60;
234
235 // Check that:
236 // - the duration for the shrunken source obs[1] is still above 40min
237 // - obs[0] does not exceed 60deg at the end of its new window
238 // - obs[0] does not exceed any limit within the new window
239
240 if (obs1_duration>=40./24/60 && obs[0].IsRangeValid(obs[0].end, obs0_end))
241 {
242 obs[0].end = obs0_end;
243 obs[1].begin = obs0_end;
244
245 cout << "First source [" << obs[0].name << "] extended to 40min" << endl;
246
247 return false;
248 }
249
250 // Try to remove the first source, check if second source fullfills all limits
251 if (obs[1].IsRangeValid(obs[0].begin, obs[0].end))
252 {
253 cout << "First source [" << obs[0].name << "] removed" << endl;
254
255 obs[1].begin = obs[0].begin;
256 obs.erase(obs.begin());
257
258 return true;
259 }
260
261 // Try to remove the second source, check if the first source fullfills all limits
262 if (obs[0].IsRangeValid(obs[1].begin, obs[1].end))
263 {
264 cout << "Second source [" << obs[1].name << "] removed" << endl;
265
266 obs[0].end = obs[1].end;
267 obs.erase(obs.begin()+1);
268
269 if (obs.size()==0 || obs[0].name!=obs[1].name)
270 return true;
271
272 obs[0].end = obs[1].end;
273 obs.erase(obs.begin()+1);
274
275 cout << "Combined first two indentical sources [" << obs[0].name << "] into one observation" << endl;
276
277 return true;
278 }
279
280 cout << "No reschedule possible within limit." << endl;
281
282 return false;
283}
284
285bool RescheduleLastSources(vector<Source> &obs)
286{
287 // If observation time is smaller than 40min for the first source
288 // extend it to 40min if zenith angle will not go above 60deg.
289 const int last = obs.size()-1;
290 if (obs.size()<2 || obs[last].duration()>=40./24/60 || obs[last].name=="SLEEP" || obs[last-1].name=="SLEEP")
291 return false;
292
293 cout << "Last source [" << obs[last].name << "] detected < 40min" << endl;
294
295 const double obs1_duration = obs[last].end - 40./24/60 - obs[last-1].begin;
296 const double obs0_begin = obs[last].end - 40./24/60;
297
298 // Check that:
299 // - the duration for the shrunken source obs[1] is still above 40min
300 // - obs[0] does not exceed 60deg at the end of its new window
301 // - obs[0] does not exceed any limit within the new window
302
303 if (obs1_duration>=40./24/60 && obs[last].IsRangeValid(obs0_begin, obs[last].begin))
304 {
305 obs[last].begin = obs0_begin;
306 obs[last-1].end = obs0_begin;
307
308 cout << "Last source [" << obs[last].name << "] extended to 40min" << endl;
309
310 return false;
311 }
312
313 // Try to remove the last source, check if second source fullfills all limits
314 if (obs[last-1].IsRangeValid(obs[last].begin, obs[last].end))
315 {
316 cout << "Last source [" << obs[last].name << "] removed" << endl;
317
318 obs[last-1].end = obs[last].end;
319 obs.pop_back();
320
321 return true;
322 }
323
324 // Try to remove the second last source, check if the first source fullfills all limits
325 if (obs[last].IsRangeValid(obs[last-1].begin, obs[last-1].end))
326 {
327 cout << "Second last source [" << obs[last-1].name << "] removed" << endl;
328
329 obs[last].begin = obs[last-1].begin;
330 obs.erase(obs.begin()+obs.size()-2);
331
332 if (obs.size()==0 || obs[last-1].name!=obs[last-2].name)
333 return true;
334
335 obs[last-2].end = obs[last-1].end;
336 obs.pop_back();
337
338 cout << "Combined last two indentical sources [" << obs[last-1].name << "] into one observation" << endl;
339
340 return true;
341 }
342
343 cout << "No reschedule possible within limit." << endl;
344
345 return false;
346}
347
348bool RescheduleIntermediateSources(vector<Source> &obs)
349{
350 for (size_t i=1; i<obs.size()-1; i++)
351 {
352 if (obs[i].duration()>=40./24/60)
353 continue;
354
355 if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
356 continue;
357
358 cout << "Intermediate source [" << obs[i].name << "] detected < 40min" << endl;
359
360 double intersection = -1;
361
362 if (obs[i-1].name=="SLEEP")
363 intersection = obs[i].begin;
364
365 if (obs[i+1].name=="SLEEP")
366 intersection = obs[i].end;
367
368 if (obs[i-1].name==obs[i+1].name)
369 intersection = obs[i].begin;
370
371 if (intersection<0)
372 {
373 const double step = 1./24/60;
374 for (double jd=obs[i].begin; jd<obs[i].end-step/2; jd+=step)
375 {
376 const SolarObjects so(jd);
377 if (obs[i-1].getThreshold(so)>=obs[i+1].getThreshold(so))
378 {
379 intersection = jd;
380 break;
381 }
382 }
383 }
384
385 if ((obs[i-1].name!="SLEEP" && !obs[i-1].IsRangeValid(obs[i-1].end, intersection)) ||
386 (obs[i+1].name!="SLEEP" && !obs[i+1].IsRangeValid(intersection, obs[i+1].begin)))
387 {
388 cout << "No reschedule possible within limits." << endl;
389 continue;
390 }
391
392 cout << "Intermediate source [" << obs[i].name << "] removed" << endl;
393
394 const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;
395
396 obs[i-1].end = intersection;
397 obs[i+1].begin = intersection;
398 obs.erase(obs.begin()+i);
399
400 i--;
401
402 if (obs.size()>1 && obs[i].name==obs[i+1].name)
403 {
404 obs[i].end = obs[i+1].end;
405 obs.erase(obs.begin()+i+1);
406
407 cout << "Combined two surrounding indentical sources [" << obs[i].name << "] into one observation" << endl;
408
409 i--;
410
411 continue;
412 }
413
414 if (underflow)
415 cout << "WARNING - Neighbor source < 40min as well." << endl;
416 }
417 return false;
418}
419
420void Print(const vector<Source> &obs, double startup_offset)
421{
422 cout << Time(obs[0].begin-startup_offset).GetAsStr() << " STARTUP\n";
423 for (const auto& src: obs)
424 {
425 string tm = Time(src.begin).GetAsStr();
426
427 if (src.preobs.size()>0)
428 {
429 for (const auto& pre: src.preobs)
430 {
431 cout << tm << " " << pre << "\n";
432 tm = " ";
433 }
434 }
435
436 cout << tm << " " << src.name << " [";
437 cout << src.duration()*24*60 << "'";
438 if (src.name!="SLEEP")
439 cout << Tools::Form("; %.1f/%.1f", src.zd(src.begin), src.zd(src.end));
440 cout << "]";
441
442 if (src.duration()*24*60<40)
443 cout << " (!)";
444
445 cout << "\n";
446 }
447 cout << Time(obs.back().end).GetAsStr() << " SHUTDOWN" << endl;
448}
449
450int main(int argc, const char* argv[])
451{
452// gROOT->SetBatch();
453
454 Configuration conf(argv[0]);
455 conf.SetPrintUsage(PrintUsage);
456 SetupConfiguration(conf);
457
458 if (!conf.DoParse(argc, argv, PrintHelp))
459 return 127;
460
461 // ------------------ Eval config ---------------------
462
463 Time time;
464 if (conf.Has("date"))
465 time.SetFromStr(conf.Get<string>("date")+" 12:00:00");
466
467 Source::max_current = conf.Get<double>("max-current");
468 Source::max_zd = conf.Get<double>("max-zd");
469
470 const double startup_offset = conf.Get<double>("startup.offset")/60/24;
471
472 const double angle_sun_set = conf.Get<double>("data-taking.start");
473 const double angle_sun_rise = conf.Get<double>("data-taking.end");
474
475 if (startup_offset<0 || startup_offset>120)
476 throw runtime_error("Only values [0;120] are allowed for startup.offset");
477
478 if (angle_sun_set>-6)
479 throw runtime_error("datataking.start not allowed before sun at -6deg");
480
481 if (angle_sun_rise>-6)
482 throw runtime_error("datataking.end not allowed after sun at -6deg");
483
484 // -12: nautical
485 // Sun set with the same date than th provided date
486 // Sun rise on the following day
487 const RstTime sun_set = GetSolarRst(time.JD()-0.5, angle_sun_set);
488 const RstTime sun_rise = GetSolarRst(time.JD()+0.5, angle_sun_rise);
489
490 const double sunset = ceil(sun_set.set*24*60)/24/60;
491 const double sunrise = floor(sun_rise.rise*24*60)/24/60;
492
493 cout << "\n";
494 cout << "Date: " << Time(floor(sunset)) << "\n";
495 cout << "Set: " << Time(sun_set.set) << "\n";
496 cout << "Rise: " << Time(sun_rise.rise) << "\n";
497 cout << endl;
498
499 // ------------- Get Sources from databasse ---------------------
500
501 const vector<string> ourcenames = conf.Vec<string>("source");
502 const vector<string> sourcenames = conf.Vec<string>("source");
503 cout << "Nsources = " << sourcenames.size() << "\n";
504
505 string query = "SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1";
506 if (sourcenames.size()>0)
507 query += " AND fSourceName IN ('" + boost::algorithm::join(sourcenames, "', '")+"')";
508
509 const string fDatabase = conf.Get<string>("source-database");
510 const mysqlpp::StoreQueryResult res =
511 Database(fDatabase).query(query).store();
512
513 // ------------------ Eval config ---------------------
514
515 vector<Source> sources;
516 for (const auto &row: res)
517 {
518 const string name = string(row[0]);
519
520 Source src(name);
521
522 src.equ.ra = double(row[1])*15;
523 src.equ.dec = double(row[2]);
524
525 src.maxzd = MyDouble(conf, "setup."+name+".max-zd");
526 src.maxcurrent = MyDouble(conf, "setup."+name+".max-current");
527 src.penalty = conf.Has("setup."+name+".penalty") ?
528 conf.Get<double>("setup."+name+".penalty") : 1;
529
530 src.preobs = conf.Vec<string>("preobs."+name);
531
532
533 cout << "[" << name << "]";
534
535 if (src.maxzd.valid)
536 cout << " Zd<" << src.maxzd.val;
537 if (src.penalty!=1)
538 cout << " Penalty=" << src.penalty;
539
540 cout << " " << boost::algorithm::join(src.preobs, "+") << endl;
541
542 /*
543 RstTime t1 = GetObjectRst(floor(sunset)-1, src.equ);
544 RstTime t2 = GetObjectRst(floor(sunset), src.equ);
545
546 src.rst.transit = t1.transit<floor(sunset) ? t2.transit : t1.transit;
547 src.rst.rise = t1.rise>src.rst.transit ? t2.rise : t1.rise;
548 src.rst.set = t1.set <src.rst.transit ? t2.set : t1.set;
549 */
550
551 sources.emplace_back(src);
552 }
553 cout << endl;
554
555 // -------------------------------------------------------------------------
556
557 vector<Source> obs;
558
559 const double step = 1./24/60;
560 for (double jd=sunset; jd<sunrise-step/2; jd+=step)
561 {
562 const SolarObjects so(jd);
563
564 vector<Source> vis;
565 for (auto& src: sources)
566 {
567 if (src.valid(so))
568 vis.emplace_back(src);
569 }
570
571 // In case no source was found, add a sleep source
572 Source src("SLEEP");
573 vis.emplace_back(src);
574
575 // Source has higher priority if minimum observation time not yet fullfilled
576 sort(vis.begin(), vis.end(), SortByThreshold);
577
578 if (obs.size()>0 && obs.back().name==vis[0].name)
579 continue;
580
581 vis[0].begin = jd;
582 obs.emplace_back(vis[0]);
583 }
584
585 if (obs.size()==0)
586 {
587 cout << "No source found." << endl;
588 return 0;
589 }
590
591 // -------------------------------------------------------------------------
592
593 for (auto it=obs.begin(); it<obs.end()-1; it++)
594 it[0].end = it[1].begin;
595 obs.back().end = sunrise;
596
597 // -------------------------------------------------------------------------
598
599 Print(obs, startup_offset);
600 cout << endl;
601
602 // -------------------------------------------------------------------------
603
604 while (RescheduleFirstSources(obs));
605 while (RescheduleLastSources(obs));
606 while (RescheduleIntermediateSources(obs));
607
608 // ---------------------------------------------------------------------
609
610 Print(obs, startup_offset);
611 cout << endl;
612
613 // ---------------------------------------------------------------------
614
615 return 1;
616
617 // ------------- Create canvases and frames ---------------------
618/*
619 // It is important to use an offset which is larger than
620 // 1970-01-01 00:00:00. This one will not work if your
621 // local time zone is positive!
622 TH1S hframe("", "", 1, Time(sunset).Mjd()*24*3600, Time(sunrise).Mjd()*24*3600);
623 hframe.SetStats(kFALSE);
624 hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
625 hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
626 hframe.GetXaxis()->CenterTitle();
627 hframe.GetYaxis()->CenterTitle();
628 hframe.GetXaxis()->SetTimeDisplay(true);
629 hframe.GetYaxis()->SetTitleSize(0.040);
630 hframe.GetXaxis()->SetTitleSize(0.040);
631 hframe.GetXaxis()->SetTitleOffset(1.1);
632 hframe.GetYaxis()->SetLabelSize(0.040);
633 hframe.GetXaxis()->SetLabelSize(0.040);
634
635 TCanvas c1;
636 c1.SetFillColor(kWhite);
637 c1.SetBorderMode(0);
638 c1.SetFrameBorderMode(0);
639 c1.SetLeftMargin(0.085);
640 c1.SetRightMargin(0.01);
641 c1.SetTopMargin(0.03);
642 c1.SetGrid();
643 hframe.GetYaxis()->SetTitle("Altitude [deg]");
644 hframe.SetMinimum(15);
645 hframe.SetMaximum(90);
646 hframe.DrawCopy();
647
648 TCanvas c2;
649 c2.SetFillColor(kWhite);
650 c2.SetBorderMode(0);
651 c2.SetFrameBorderMode(0);
652 c2.SetLeftMargin(0.085);
653 c2.SetRightMargin(0.01);
654 c2.SetTopMargin(0.03);
655 c2.SetGrid();
656 hframe.GetYaxis()->SetTitle("Predicted Current [#muA]");
657 hframe.SetMinimum(0);
658 hframe.SetMaximum(100);
659 hframe.DrawCopy();
660
661 TCanvas c3;
662 c3.SetFillColor(kWhite);
663 c3.SetBorderMode(0);
664 c3.SetFrameBorderMode(0);
665 c3.SetLeftMargin(0.085);
666 c3.SetRightMargin(0.01);
667 c3.SetTopMargin(0.03);
668 c3.SetGrid();
669 c3.SetLogy();
670 hframe.GetYaxis()->SetTitle("Estimated relative threshold");
671 hframe.GetYaxis()->SetMoreLogLabels();
672 hframe.SetMinimum(0.9);
673 hframe.SetMaximum(11);
674 hframe.DrawCopy();
675
676 TCanvas c4;
677 c4.SetFillColor(kWhite);
678 c4.SetBorderMode(0);
679 c4.SetFrameBorderMode(0);
680 c4.SetLeftMargin(0.085);
681 c4.SetRightMargin(0.01);
682 c4.SetTopMargin(0.03);
683 c4.SetGrid();
684 hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
685 hframe.SetMinimum(0);
686 hframe.SetMaximum(180);
687 hframe.DrawCopy();
688 Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
689 Int_t style[] = { kSolid, kDashed, kDotted };
690
691 // TLegend leg(0, 0, 1, 1);
692
693 // ------------- Loop over sources ---------------------
694
695 for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
696 {
697 // Eval row
698 const string name = (*v)[0].c_str();
699
700 EquPosn pos;
701 pos.ra = double((*v)[1])*15;
702 pos.dec = double((*v)[2]);
703
704 // Loop over 24 hours
705 for (double h=0; h<1; h+=1./(24*12))
706 {
707 const SolarObjects so(jd+h);
708
709 // get local position of source
710 const HrzPosn hrz = GetHrzFromEqu(pos, so.fJD);
711
712 if (v==res.begin())
713 cout << Time(so.fJD) <<" " << 90-so.fMoonHrz.alt << endl;
714
715 const double cur = FACT::PredictI(so, pos);
716
717 // Relative energy threshold prediction
718 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
719
720 // Add points to curve
721 const double axis = Time(so.fJD).Mjd()*24*3600;
722
723 // If there is a gap of more than one bin, start a new curve
724
725 // Add data
726 if (no_limits || cur<max_current)
727 g1.SetPoint(g1.GetN(), axis, hrz.alt);
728
729 if (no_limits || 90-hrz.alt<max_zd)
730 g2.SetPoint(g2.GetN(), axis, cur);
731
732 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
733 g3.SetPoint(g3.GetN(), axis, ratio*pow(cur/6.2, 0.394));
734
735 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
736 {
737 const double angle = GetAngularSeparation(so.fMoonEqu, pos);
738 g4.SetPoint(g4.GetN(), axis, angle);
739 }
740
741 if (cnt==0)
742 gm.SetPoint(gm.GetN(), axis, so.fMoonHrz.alt);
743 }
744 }
745
746 // Save three plots
747 TCanvas c5;
748 c5.SetFillColor(kWhite);
749 c5.SetBorderMode(0);
750 c5.SetFrameBorderMode(0);
751 leg.Draw();
752
753 const string t = Time(jd).GetAsStr("%Y%m%d");
754
755 c1.SaveAs((t+"-ZenithDistance.eps").c_str());
756 c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
757 c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
758 c4.SaveAs((t+"-MoonDist.eps").c_str());
759 c5.SaveAs((t+"-Legend.eps").c_str());
760
761 c1.SaveAs((t+"-ZenithDistance.root").c_str());
762 c2.SaveAs((t+"-PredictedCurrent.root").c_str());
763 c3.SaveAs((t+"-RelativeThreshold.root").c_str());
764 c4.SaveAs((t+"-MoonDist.root").c_str());
765
766 c1.Print((t+".pdf(").c_str(), "pdf");
767 c2.Print((t+".pdf" ).c_str(), "pdf");
768 c3.Print((t+".pdf" ).c_str(), "pdf");
769 c4.Print((t+".pdf" ).c_str(), "pdf");
770 c5.Print((t+".pdf)").c_str(), "pdf");
771*/
772}
Note: See TracBrowser for help on using the repository browser.