source: trunk/FACT++/src/makeplots.cc@ 16991

Last change on this file since 16991 was 16969, checked in by tbretz, 11 years ago
Fixed a stupid bug: ln__equ_posn stores RS in degree not in hours, therefore the result of Angle() was wrong
File size: 10.8 KB
Line 
1#include "externals/nova.h"
2
3#include "Database.h"
4
5#include "Time.h"
6#include "Configuration.h"
7
8#include <TROOT.h>
9#include <TH1.h>
10#include <TGraph.h>
11#include <TCanvas.h>
12#include <TLegend.h>
13
14using namespace std;
15using namespace Nova;
16
17// ------------------------------------------------------------------------
18
19void CheckForGap(TCanvas &c, TGraph &g, double axis)
20{
21 if (g.GetN()==0 || axis-g.GetX()[g.GetN()-1]<450)
22 return;
23
24 c.cd();
25 ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
26 while (g.GetN())
27 g.RemovePoint(0);
28}
29
30void DrawClone(TCanvas &c, TGraph &g)
31{
32 if (g.GetN()==0)
33 return;
34
35 c.cd();
36 ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
37}
38
39// ========================================================================
40// ========================================================================
41// ========================================================================
42
43void SetupConfiguration(Configuration &conf)
44{
45 po::options_description control("Smart FACT");
46 control.add_options()
47 ("ra", var<double>(), "Source right ascension")
48 ("dec", var<double>(), "Source declination")
49 ("date-time", var<string>(), "SQL time (UTC)")
50 ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
51 ("max-current", var<double>(75), "Maximum current to display in other plots.")
52 ("max-zd", var<double>(75), "Maximum zenith distance to display in other plots")
53 ("no-limits", po_switch(), "Switch off limits in plots")
54 ;
55
56 po::positional_options_description p;
57 p.add("date-time", 1); // The first positional options
58
59 conf.AddOptions(control);
60 conf.SetArgumentPositions(p);
61}
62
63void PrintUsage()
64{
65 cout <<
66 "makeplots - The astronomy plotter\n"
67 "\n"
68 "Calculates several plots for the sources in the database\n"
69 "helpful or needed for scheduling. The Plot is always calculated\n"
70 "for the night which starts at the same so. So no matter if\n"
71 "you specify '1974-09-09 00:00:00' or '1974-09-09 21:56:00'\n"
72 "the plots will refer to the night 1974-09-09/1974-09-10.\n"
73 "The advantage is that specification of the date as in\n"
74 "1974-09-09 is enough. Time axis starts and ends at nautical\n"
75 "twilight which is 12deg below horizon.\n"
76 "\n"
77 "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
78 cout << endl;
79}
80
81int main(int argc, const char* argv[])
82{
83 gROOT->SetBatch();
84
85 Configuration conf(argv[0]);
86 conf.SetPrintUsage(PrintUsage);
87 SetupConfiguration(conf);
88
89 if (!conf.DoParse(argc, argv))
90 return 127;
91
92 if (conf.Has("ra")^conf.Has("dec"))
93 {
94 cout << "ERROR - Either --ra or --dec missing." << endl;
95 return 1;
96 }
97
98 // ------------------ Eval config ---------------------
99
100 Time time;
101 if (conf.Has("date-time"))
102 time.SetFromStr(conf.Get<string>("date-time"));
103
104 const double max_current = conf.Get<double>("max-current");
105 const double max_zd = conf.Get<double>("max-zd");
106 const double no_limits = conf.Get<bool>("no-limits");
107
108 // -12: nautical
109 // Sun set with the same date than th provided date
110 // Sun rise on the following day
111 const RstTime sun_set = GetSolarRst(time.JD()-0.5, -12);
112 const RstTime sun_rise = GetSolarRst(time.JD()+0.5, -12);
113
114 const double jd = floor(time.Mjd())+2400001;
115 const double mjd = floor(time.Mjd())+49718+0.5; // root time
116
117 cout << "Time: " << time << endl;
118 cout << "Base: " << Time(jd-0.5) << endl;
119 cout << "Set: " << Time(sun_set.set) << endl;
120 cout << "Rise: " << Time(sun_rise.rise) << endl;
121
122 const double jd0 = fmod(sun_set.set, 1); // ~0.3
123 const double jd1 = fmod(sun_rise.rise, 1); // ~0.8
124
125 const string fDatabase = conf.Get<string>("source-database");
126
127 // ------------------ Precalc moon coord ---------------------
128
129 vector<pair<EquPosn, double>> fMoonCoords;
130
131 for (double h=0; h<1; h+=1./(24*12))
132 {
133 const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
134 const double disk = GetLunarDisk(jd+h);
135
136 fMoonCoords.emplace_back(moon, disk);
137 }
138
139 // ------------- Get Sources from databasse ---------------------
140
141 const mysqlpp::StoreQueryResult res =
142 Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM source").store();
143
144 // ------------- Create canvases and frames ---------------------
145
146 // It is important to use an offset which is larger than
147 // 1970-01-01 00:00:00. This one will not work if your
148 // local time zone is positive!
149 TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
150 hframe.SetStats(kFALSE);
151 hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
152 hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
153 hframe.GetXaxis()->CenterTitle();
154 hframe.GetYaxis()->CenterTitle();
155 hframe.GetXaxis()->SetTimeDisplay(true);
156 hframe.GetYaxis()->SetTitleSize(0.040);
157 hframe.GetXaxis()->SetTitleSize(0.040);
158 hframe.GetXaxis()->SetTitleOffset(1.1);
159 hframe.GetYaxis()->SetLabelSize(0.040);
160 hframe.GetXaxis()->SetLabelSize(0.040);
161
162 TCanvas c1;
163 c1.SetFillColor(kWhite);
164 c1.SetBorderMode(0);
165 c1.SetFrameBorderMode(0);
166 c1.SetLeftMargin(0.085);
167 c1.SetRightMargin(0.01);
168 c1.SetTopMargin(0.03);
169 c1.SetGrid();
170 hframe.GetYaxis()->SetTitle("Altitude [deg]");
171 hframe.SetMinimum(15);
172 hframe.SetMaximum(90);
173 hframe.DrawCopy();
174
175 TCanvas c2;
176 c2.SetFillColor(kWhite);
177 c2.SetBorderMode(0);
178 c2.SetFrameBorderMode(0);
179 c2.SetLeftMargin(0.085);
180 c2.SetRightMargin(0.01);
181 c2.SetTopMargin(0.03);
182 c2.SetGrid();
183 hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
184 hframe.SetMinimum(0);
185 hframe.SetMaximum(100);
186 hframe.DrawCopy();
187
188 TCanvas c3;
189 c3.SetFillColor(kWhite);
190 c3.SetBorderMode(0);
191 c3.SetFrameBorderMode(0);
192 c3.SetLeftMargin(0.085);
193 c3.SetRightMargin(0.01);
194 c3.SetTopMargin(0.03);
195 c3.SetGrid();
196 c3.SetLogy();
197 hframe.GetYaxis()->SetTitle("Estimated relative threshold");
198 hframe.SetMinimum(0.9);
199 hframe.SetMaximum(180);
200 hframe.DrawCopy();
201
202 TCanvas c4;
203 c4.SetFillColor(kWhite);
204 c4.SetBorderMode(0);
205 c4.SetFrameBorderMode(0);
206 c4.SetLeftMargin(0.085);
207 c4.SetRightMargin(0.01);
208 c4.SetTopMargin(0.03);
209 c4.SetGrid();
210 hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
211 hframe.SetMinimum(0);
212 hframe.SetMaximum(180);
213 hframe.DrawCopy();
214
215 Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
216 Int_t style[] = { kSolid, kDashed, kDotted };
217
218 TLegend leg(0, 0, 1, 1);
219
220 // ------------- Loop over sources ---------------------
221
222 Int_t cnt=0;
223 for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
224 {
225 // Eval row
226 const string name = (*v)[0].c_str();
227
228 EquPosn pos;
229 pos.ra = double((*v)[1])*15;
230 pos.dec = double((*v)[2]);
231
232 // Create graphs
233 TGraph g1, g2, g3, g4, gr, gm;
234 g1.SetName(name.data());
235 g2.SetName(name.data());
236 g3.SetName(name.data());
237 g4.SetName(name.data());
238 g1.SetLineWidth(2);
239 g2.SetLineWidth(2);
240 g3.SetLineWidth(2);
241 g4.SetLineWidth(2);
242 gm.SetLineWidth(1);
243 g1.SetLineStyle(style[cnt/6]);
244 g2.SetLineStyle(style[cnt/6]);
245 g3.SetLineStyle(style[cnt/6]);
246 g4.SetLineStyle(style[cnt/6]);
247 g1.SetLineColor(color[cnt%6]);
248 g2.SetLineColor(color[cnt%6]);
249 g3.SetLineColor(color[cnt%6]);
250 g4.SetLineColor(color[cnt%6]);
251 gm.SetLineColor(kYellow);
252
253 if (cnt==0)
254 leg.AddEntry(gm.Clone(), "Moon", "l");
255 leg.AddEntry(g1.Clone(), name.data(), "l");
256
257 // Loop over 24 hours
258 int i=0;
259 for (double h=0; h<1; h+=1./(24*12), i++)
260 {
261 // check if it is betwene sun-rise and sun-set
262 if (h<jd0 || h>jd1)
263 continue;
264
265 // get local position of source
266 const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
267
268 // Get moon properties and
269 const EquPosn moon = fMoonCoords[i].first;
270 const double disk = fMoonCoords[i].second;
271 const HrzPosn hrzm = GetHrzFromEqu(moon, jd+h);
272
273 if (v==res.begin())
274 cout << Time(jd+h) <<" " << 90-hrzm.alt << endl;
275
276 // Distance between source and moon
277 const double angle = GetAngularSeparation(moon, pos);
278
279 // Current prediction
280 const double cang = sin(angle *M_PI/180);
281 const double calt = sin(hrzm.alt*M_PI/180);
282
283 const double lc = cang==0 ? -1 : calt*pow(disk, 3)/sqrt(cang);
284 const double cur = lc>0 ? 8.1+94.6*lc : 8.1;
285
286 // Relative energy threshold prediction
287 //const double cs = cos((90+hrz.alt)*M_PI/180);
288 //const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;
289 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
290
291 // Add points to curve
292 const double axis = (mjd+h)*24*3600;
293
294 // If there is a gap of more than one bin, start a new curve
295 CheckForGap(c1, g1, axis);
296 CheckForGap(c1, gm, axis);
297 CheckForGap(c2, g2, axis);
298 CheckForGap(c3, g3, axis);
299 CheckForGap(c4, g4, axis);
300
301 // Add data
302 if (no_limits || cur<max_current)
303 g1.SetPoint(g1.GetN(), axis, hrz.alt);
304
305 if (no_limits || 90-hrz.alt<max_zd)
306 g2.SetPoint(g2.GetN(), axis, cur);
307
308 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
309 g3.SetPoint(g3.GetN(), axis, ratio*cur/8.1);
310
311 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
312 g4.SetPoint(g4.GetN(), axis, angle);
313
314 if (cnt==0)
315 gm.SetPoint(gm.GetN(), axis, hrzm.alt);
316 }
317
318 if (cnt==0)
319 DrawClone(c1, gm);
320
321 DrawClone(c1, g1);
322 DrawClone(c2, g2);
323 DrawClone(c3, g3);
324 DrawClone(c4, g4);
325 }
326
327
328 // Save three plots
329 TCanvas c5;
330 c5.SetFillColor(kWhite);
331 c5.SetBorderMode(0);
332 c5.SetFrameBorderMode(0);
333 leg.Draw();
334
335 const string t = Time(jd).GetAsStr("%Y%m%d");
336
337 c1.SaveAs((t+"-ZenithDistance.eps").c_str());
338 c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
339 c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
340 c4.SaveAs((t+"-MoonDist.eps").c_str());
341 c5.SaveAs((t+"-Legend.eps").c_str());
342
343 c1.SaveAs((t+"-ZenithDistance.root").c_str());
344 c2.SaveAs((t+"-PredictedCurrent.root").c_str());
345 c3.SaveAs((t+"-RelativeThreshold.root").c_str());
346 c4.SaveAs((t+"-MoonDist.root").c_str());
347
348 return 0;
349}
Note: See TracBrowser for help on using the repository browser.