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

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