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

Last change on this file since 14711 was 14705, checked in by tbretz, 12 years ago
Added more command-line options; added a plot for the estimated relative energy threshold; apply zd and current limits on either plot; small improvements to the plot layout
File size: 9.4 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
44// ========================================================================
45// ========================================================================
46// ========================================================================
47
48void SetupConfiguration(Configuration &conf)
49{
50 po::options_description control("Smart FACT");
51 control.add_options()
52 ("ra", var<double>(), "Source right ascension")
53 ("dec", var<double>(), "Source declination")
54 ("date-time", var<string>(), "SQL time (UTC)")
55 ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
56 ("max-current", var<double>(75), "Maximum current to display in other plots.")
57 ("max-zd", var<double>(75), "Maximum zenith distance to display in other plots")
58 ("no-limits", po_switch(), "Switch off limits in plots")
59 ;
60
61 po::positional_options_description p;
62 p.add("date-time", 1); // The first positional options
63
64 conf.AddOptions(control);
65 conf.SetArgumentPositions(p);
66}
67
68void PrintUsage()
69{
70 cout <<
71 "makeplots - The astronomy plotter\n"
72 "\n"
73 "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
74 cout << endl;
75}
76
77int main(int argc, const char* argv[])
78{
79 gROOT->SetBatch();
80
81 Configuration conf(argv[0]);
82 conf.SetPrintUsage(PrintUsage);
83 SetupConfiguration(conf);
84
85 if (!conf.DoParse(argc, argv))
86 return 127;
87
88 if (conf.Has("ra")^conf.Has("dec"))
89 {
90 cout << "ERROR - Either --ra or --dec missing." << endl;
91 return 1;
92 }
93
94 // ------------------ Eval config ---------------------
95
96 const double lon = -(17.+53./60+26.525/3600);
97 const double lat = 28.+45./60+42.462/3600;
98
99 ln_lnlat_posn observer;
100 observer.lng = lon;
101 observer.lat = lat;
102
103 Time time;
104 if (conf.Has("date-time"))
105 time.SetFromStr(conf.Get<string>("date-time"));
106
107 const double max_current = conf.Get<double>("max-current");
108 const double max_zd = conf.Get<double>("max-zd");
109 const double no_limits = conf.Get<bool>("no-limits");
110
111 ln_rst_time sun_astronomical;
112 ln_get_solar_rst_horizon(time.JD(), &observer, -12, &sun_astronomical);
113
114 const double jd = floor(time.JD());
115 const double mjd = floor(time.Mjd())+49718-0.5;
116
117 const double jd0 = fmod(sun_astronomical.set, 1);
118 const double jd1 = fmod(sun_astronomical.rise, 1);
119
120 cout << "JD0=" << jd+jd0 << endl;
121 cout << "JD1=" << jd+jd1 << endl;
122
123 const string fDatabase = conf.Get<string>("source-database");
124
125 // ------------------ Precalc moon coord ---------------------
126
127 vector<pair<ln_equ_posn, double>> fMoonCoords;
128
129 for (double h=0; h<1; h+=1./(24*12))
130 {
131 ln_equ_posn moon;
132 ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);
133
134 const double disk = ln_get_lunar_disk(jd+h);
135
136 fMoonCoords.push_back(make_pair(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 TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
147 hframe.SetStats(kFALSE);
148 hframe.GetXaxis()->SetTimeFormat("%Hh%M'");
149 hframe.GetXaxis()->SetTitle("Time [UTC]");
150 hframe.GetXaxis()->CenterTitle();
151 hframe.GetYaxis()->CenterTitle();
152 hframe.GetXaxis()->SetTimeDisplay(true);
153 hframe.GetYaxis()->SetTitleSize(0.040);
154 hframe.GetXaxis()->SetTitleSize(0.040);
155 hframe.GetXaxis()->SetTitleOffset(1.1);
156 hframe.GetYaxis()->SetLabelSize(0.040);
157 hframe.GetXaxis()->SetLabelSize(0.040);
158
159 TCanvas c1;
160 gPad->SetLeftMargin(0.085);
161 gPad->SetRightMargin(0.01);
162 gPad->SetTopMargin(0.03);
163 gPad->SetGrid();
164 hframe.GetYaxis()->SetTitle("Altitude [deg]");
165 hframe.SetMinimum(15);
166 hframe.SetMaximum(90);
167 hframe.DrawCopy();
168
169 TCanvas c2;
170 gPad->SetLeftMargin(0.085);
171 gPad->SetRightMargin(0.01);
172 gPad->SetTopMargin(0.03);
173 gPad->SetGrid();
174 hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
175 hframe.SetMinimum(0);
176 hframe.SetMaximum(100);
177 hframe.DrawCopy();
178
179 TCanvas c3;
180 gPad->SetLeftMargin(0.085);
181 gPad->SetRightMargin(0.01);
182 gPad->SetTopMargin(0.03);
183 gPad->SetGrid();
184 gPad->SetLogy();
185 hframe.GetYaxis()->SetTitle("Estimated relative threshold");
186 hframe.SetMinimum(0.9);
187 hframe.SetMaximum(180);
188 hframe.DrawCopy();
189
190 Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
191 Int_t style[] = { kSolid, kDashed, kDotted };
192
193 TLegend leg(0, 0, 1, 1);
194
195 // ------------- Loop over sources ---------------------
196
197 Int_t cnt=0;
198 for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
199 {
200 // Eval row
201 const string name = (*v)[0].c_str();
202
203 ln_equ_posn pos;
204 pos.ra = double((*v)[1])*15;
205 pos.dec = double((*v)[2]);
206
207 // Create graphs
208 TGraph g1, g2, g3;
209 g1.SetName(name.data());
210 g2.SetName(name.data());
211 g3.SetName(name.data());
212 g1.SetLineWidth(2);
213 g2.SetLineWidth(2);
214 g3.SetLineWidth(2);
215 g1.SetLineStyle(style[cnt/6]);
216 g2.SetLineStyle(style[cnt/6]);
217 g3.SetLineStyle(style[cnt/6]);
218 g1.SetLineColor(color[cnt%6]);
219 g2.SetLineColor(color[cnt%6]);
220 g3.SetLineColor(color[cnt%6]);
221
222 // Loop over 24 hours
223 int i=0;
224 for (double h=0; h<1; h+=1./(24*12), i++)
225 {
226 // check if it is betwene sun-rise and sun-set
227 if (h<jd0 || h>jd1)
228 continue;
229
230 // get local position of source
231 ln_hrz_posn hrz;
232 ln_get_hrz_from_equ(&pos, &observer, jd+h, &hrz);
233
234 // Get moon properties and
235 ln_equ_posn moon = fMoonCoords[i].first;
236 const double disk = fMoonCoords[i].second;
237
238 ln_hrz_posn hrzm;
239 ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrzm);
240
241 // Distance between source and moon
242 const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec);
243
244 // Current prediction
245 const double lc = angle*hrzm.alt*pow(disk, 6)/360/360;
246 const double cur = lc>0 ? 7.7+4942*lc : 7.7;
247
248 // Relative energy threshold prediction
249 const double cs = cos((90+hrz.alt)*M_PI/180);
250 const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;
251
252 // Add points to curve
253 const double axis = (mjd+h)*24*3600;
254
255 // If there is a gap of more than one bin, start a new curve
256 if (g1.GetN()>0 && axis-g1.GetX()[g1.GetN()-1]>450)
257 {
258 c1.cd();
259 ((TGraph*)g1.DrawClone("C"))->SetBit(kCanDelete);
260 while (g1.GetN())
261 g1.RemovePoint(0);
262 }
263
264 if (g2.GetN()>0 && axis-g2.GetX()[g2.GetN()-1]>450)
265 {
266 c2.cd();
267 ((TGraph*)g2.DrawClone("C"))->SetBit(kCanDelete);
268 while (g2.GetN())
269 g2.RemovePoint(0);
270 }
271
272 if (g3.GetN()>0 && axis-g3.GetX()[g3.GetN()-1]>450)
273 {
274 c3.cd();
275 ((TGraph*)g3.DrawClone("C"))->SetBit(kCanDelete);
276 while (g3.GetN())
277 g3.RemovePoint(0);
278 }
279
280 // Add data
281 if (no_limits || cur<max_current)
282 g1.SetPoint(g1.GetN(), axis, hrz.alt);
283
284 if (no_limits || 90-hrz.alt<max_zd)
285 g2.SetPoint(g2.GetN(), axis, cur);
286
287 if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
288 g3.SetPoint(g3.GetN(), axis, ratio*cur/7.7);
289 }
290
291 // Add graphs to canvases and add corresponding entry to legend
292 TGraph *g = 0;
293
294 c1.cd();
295 g = (TGraph*)g1.DrawClone("C");
296 g->SetBit(kCanDelete);
297
298 c2.cd();
299 g = (TGraph*)g2.DrawClone("C");
300 g->SetBit(kCanDelete);
301
302 c3.cd();
303 g = (TGraph*)g3.DrawClone("C");
304 g->SetBit(kCanDelete);
305
306 leg.AddEntry(g, name.data(), "l");
307 }
308
309 // Save three plots
310 TCanvas c4;
311 leg.Draw();
312
313 c1.SaveAs("test1.eps");
314 c2.SaveAs("test2.eps");
315 c3.SaveAs("test3.eps");
316 c4.SaveAs("legend.eps");
317
318 c1.SaveAs("test1.root");
319 c2.SaveAs("test2.root");
320 c3.SaveAs("test3.root");
321
322 return 0;
323}
Note: See TracBrowser for help on using the repository browser.