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

Last change on this file since 17624 was 17367, checked in by tbretz, 11 years ago
Missed the correct calculation of the coefficient derived from the angle to the moon.
File size: 11.1 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 // Distance between earth and moon relative to major semi-axis
280 const double dist = GetLunarEarthDist(jd+h)/384400;
281
282 // Current prediction
283 const double calt = sin(hrzm.alt*M_PI/180);
284 const double cang = 1-sin(angle*M_PI/180);
285
286 //const double lc = cang==0 ? -1 : calt*pow(disk, 2.2)*pow(dist, -2);
287 //const double cur = lc>0 ? 4+103*lc : 4;
288
289 const double lc = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2);
290 const double cur = lc>0 ? 7.2 + 69*lc : 7.2;
291
292 // Relative energy threshold prediction
293 //const double cs = cos((90+hrz.alt)*M_PI/180);
294 //const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;
295 const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
296
297 // Add points to curve
298 const double axis = (mjd+h)*24*3600;
299
300 // If there is a gap of more than one bin, start a new curve
301 CheckForGap(c1, g1, axis);
302 CheckForGap(c1, gm, axis);
303 CheckForGap(c2, g2, axis);
304 CheckForGap(c3, g3, axis);
305 CheckForGap(c4, g4, axis);
306
307 // Add data
308 if (no_limits || cur<max_current)
309 g1.SetPoint(g1.GetN(), axis, hrz.alt);
310
311 if (no_limits || 90-hrz.alt<max_zd)
312 g2.SetPoint(g2.GetN(), axis, cur);
313
314 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
315 g3.SetPoint(g3.GetN(), axis, ratio*cur/8.1);
316
317 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
318 g4.SetPoint(g4.GetN(), axis, angle);
319
320 if (cnt==0)
321 gm.SetPoint(gm.GetN(), axis, hrzm.alt);
322 }
323
324 if (cnt==0)
325 DrawClone(c1, gm);
326
327 DrawClone(c1, g1);
328 DrawClone(c2, g2);
329 DrawClone(c3, g3);
330 DrawClone(c4, g4);
331 }
332
333
334 // Save three plots
335 TCanvas c5;
336 c5.SetFillColor(kWhite);
337 c5.SetBorderMode(0);
338 c5.SetFrameBorderMode(0);
339 leg.Draw();
340
341 const string t = Time(jd).GetAsStr("%Y%m%d");
342
343 c1.SaveAs((t+"-ZenithDistance.eps").c_str());
344 c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
345 c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
346 c4.SaveAs((t+"-MoonDist.eps").c_str());
347 c5.SaveAs((t+"-Legend.eps").c_str());
348
349 c1.SaveAs((t+"-ZenithDistance.root").c_str());
350 c2.SaveAs((t+"-PredictedCurrent.root").c_str());
351 c3.SaveAs((t+"-RelativeThreshold.root").c_str());
352 c4.SaveAs((t+"-MoonDist.root").c_str());
353
354 return 0;
355}
Note: See TracBrowser for help on using the repository browser.