Context Navigation

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

Visit:

Last change on this file since 17654 was 17654, checked in by tbretz, 11 years ago
Removed the obsolete ra/dec arguments; output pdf as well
File size: 11.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
14	using namespace std;
15	using namespace Nova;
16
17	// ------------------------------------------------------------------------
18
19	void 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
30	void 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
43	void 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>(100), "Maximum current to display in other plots.")
52	("max-zd", var<double>(50), "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
63	void 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\n";
78	// "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
79	cout << endl;
80	}
81
82	int main(int argc, const char* argv[])
83	{
84	gROOT->SetBatch();
85
86	Configuration conf(argv[0]);
87	conf.SetPrintUsage(PrintUsage);
88	SetupConfiguration(conf);
89
90	if (!conf.DoParse(argc, argv))
91	return 127;
92
93	// if (conf.Has("ra")^conf.Has("dec"))
94	// {
95	// cout << "ERROR - Either --ra or --dec missing." << endl;
96	// return 1;
97	// }
98
99	// ------------------ Eval config ---------------------
100
101	Time time;
102	if (conf.Has("date-time"))
103	time.SetFromStr(conf.Get<string>("date-time"));
104
105	const double max_current = conf.Get<double>("max-current");
106	const double max_zd = conf.Get<double>("max-zd");
107	const double no_limits = conf.Get<bool>("no-limits");
108
109	// -12: nautical
110	// Sun set with the same date than th provided date
111	// Sun rise on the following day
112	const RstTime sun_set = GetSolarRst(time.JD()-0.5, -12);
113	const RstTime sun_rise = GetSolarRst(time.JD()+0.5, -12);
114
115	const double jd = floor(time.Mjd())+2400001;
116	const double mjd = floor(time.Mjd())+49718+0.5; // root time
117
118	cout << "Time: " << time << endl;
119	cout << "Base: " << Time(jd-0.5) << endl;
120	cout << "Set: " << Time(sun_set.set) << endl;
121	cout << "Rise: " << Time(sun_rise.rise) << endl;
122
123	const double jd0 = fmod(sun_set.set, 1); // ~0.3
124	const double jd1 = fmod(sun_rise.rise, 1); // ~0.8
125
126	const string fDatabase = conf.Get<string>("source-database");
127
128	// ------------------ Precalc moon coord ---------------------
129
130	vector<tuple<EquPosn, double, double>> fMoonCoords;
131	vector<EquPosn> fSunCoords;
132
133	for (double h=0; h<1; h+=1./(24*12))
134	{
135	const EquPosn sun = GetSolarEquCoords(jd+h);
136	const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
137	const double disk = GetLunarDisk(jd+h);
138	const double edist = GetLunarEarthDist(jd+h)/384400;
139
140	fMoonCoords.emplace_back(moon, disk, edist);
141	fSunCoords.emplace_back(sun);
142	}
143
144	// ------------- Get Sources from databasse ---------------------
145
146	const mysqlpp::StoreQueryResult res =
147	Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1").store();
148
149	// ------------- Create canvases and frames ---------------------
150
151	// It is important to use an offset which is larger than
152	// 1970-01-01 00:00:00. This one will not work if your
153	// local time zone is positive!
154	TH1S hframe("", "", 1, (mjd+jd0)243600, (mjd+jd1)243600);
155	hframe.SetStats(kFALSE);
156	hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
157	hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
158	hframe.GetXaxis()->CenterTitle();
159	hframe.GetYaxis()->CenterTitle();
160	hframe.GetXaxis()->SetTimeDisplay(true);
161	hframe.GetYaxis()->SetTitleSize(0.040);
162	hframe.GetXaxis()->SetTitleSize(0.040);
163	hframe.GetXaxis()->SetTitleOffset(1.1);
164	hframe.GetYaxis()->SetLabelSize(0.040);
165	hframe.GetXaxis()->SetLabelSize(0.040);
166
167	TCanvas c1;
168	c1.SetFillColor(kWhite);
169	c1.SetBorderMode(0);
170	c1.SetFrameBorderMode(0);
171	c1.SetLeftMargin(0.085);
172	c1.SetRightMargin(0.01);
173	c1.SetTopMargin(0.03);
174	c1.SetGrid();
175	hframe.GetYaxis()->SetTitle("Altitude [deg]");
176	hframe.SetMinimum(15);
177	hframe.SetMaximum(90);
178	hframe.DrawCopy();
179
180	TCanvas c2;
181	c2.SetFillColor(kWhite);
182	c2.SetBorderMode(0);
183	c2.SetFrameBorderMode(0);
184	c2.SetLeftMargin(0.085);
185	c2.SetRightMargin(0.01);
186	c2.SetTopMargin(0.03);
187	c2.SetGrid();
188	hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
189	hframe.SetMinimum(0);
190	hframe.SetMaximum(100);
191	hframe.DrawCopy();
192
193	TCanvas c3;
194	c3.SetFillColor(kWhite);
195	c3.SetBorderMode(0);
196	c3.SetFrameBorderMode(0);
197	c3.SetLeftMargin(0.085);
198	c3.SetRightMargin(0.01);
199	c3.SetTopMargin(0.03);
200	c3.SetGrid();
201	c3.SetLogy();
202	hframe.GetYaxis()->SetTitle("Estimated relative threshold");
203	hframe.SetMinimum(0.9);
204	hframe.SetMaximum(180);
205	hframe.DrawCopy();
206
207	TCanvas c4;
208	c4.SetFillColor(kWhite);
209	c4.SetBorderMode(0);
210	c4.SetFrameBorderMode(0);
211	c4.SetLeftMargin(0.085);
212	c4.SetRightMargin(0.01);
213	c4.SetTopMargin(0.03);
214	c4.SetGrid();
215	hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
216	hframe.SetMinimum(0);
217	hframe.SetMaximum(180);
218	hframe.DrawCopy();
219
220	Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
221	Int_t style[] = { kSolid, kDashed, kDotted };
222
223	TLegend leg(0, 0, 1, 1);
224
225	// ------------- Loop over sources ---------------------
226
227	Int_t cnt=0;
228	for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
229	{
230	// Eval row
231	const string name = (*v)[0].c_str();
232
233	EquPosn pos;
234	pos.ra = double((v)[1])15;
235	pos.dec = double((*v)[2]);
236
237	// Create graphs
238	TGraph g1, g2, g3, g4, gr, gm;
239	g1.SetName(name.data());
240	g2.SetName(name.data());
241	g3.SetName(name.data());
242	g4.SetName(name.data());
243	g1.SetLineWidth(2);
244	g2.SetLineWidth(2);
245	g3.SetLineWidth(2);
246	g4.SetLineWidth(2);
247	gm.SetLineWidth(1);
248	g1.SetLineStyle(style[cnt/6]);
249	g2.SetLineStyle(style[cnt/6]);
250	g3.SetLineStyle(style[cnt/6]);
251	g4.SetLineStyle(style[cnt/6]);
252	g1.SetLineColor(color[cnt%6]);
253	g2.SetLineColor(color[cnt%6]);
254	g3.SetLineColor(color[cnt%6]);
255	g4.SetLineColor(color[cnt%6]);
256	gm.SetLineColor(kYellow);
257
258	if (cnt==0)
259	leg.AddEntry(gm.Clone(), "Moon", "l");
260	leg.AddEntry(g1.Clone(), name.data(), "l");
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	const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
272
273	// Get sun and moon properties and
274	const EquPosn sun = fSunCoords[i];
275	const EquPosn moon = get<0>(fMoonCoords[i]);
276	const double disk = get<1>(fMoonCoords[i]);
277	const double edist = get<2>(fMoonCoords[i]);
278	const HrzPosn hrzm = GetHrzFromEqu(moon, jd+h);
279	const ZdAzPosn hrzs = GetHrzFromEqu(sun, jd+h);
280
281	if (v==res.begin())
282	cout << Time(jd+h) <<" " << 90-hrzm.alt << endl;
283
284	// Distance between source and moon
285	const double angle = GetAngularSeparation(moon, pos);
286
287	// Current prediction
288	const double sin_malt = hrzm.alt<0 ? 0 : sin(hrzm.alt*M_PI/180);
289	const double cos_mdist = cos(angle*M_PI/180);
290	const double sin_szd = sin(hrzs.zd*M_PI/180);
291
292	const double k0 = pow(disk, 2.63);
293	const double k1 = pow(sin_malt, 0.60);
294	const double k2 = pow(edist, -2.00);
295	const double k3 = exp(0.67cos_mdistcos_mdistcos_mdistcos_mdist);
296	const double k4 = exp(-97.8+105.8sin_szdsin_szd);
297
298	const double cur = 6.2 + 95.7k0k1k2k3 + k4;
299
300	// Relative energy threshold prediction
301	//const double cs = cos((90+hrz.alt)*M_PI/180);
302	//const double ratio = (10.sqrt(409600.cscs+9009.) + 6400.cs - 60.)/10.;
303	const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
304
305	// Add points to curve
306	const double axis = (mjd+h)243600;
307
308	// If there is a gap of more than one bin, start a new curve
309	CheckForGap(c1, g1, axis);
310	CheckForGap(c1, gm, axis);
311	CheckForGap(c2, g2, axis);
312	CheckForGap(c3, g3, axis);
313	CheckForGap(c4, g4, axis);
314
315	// Add data
316	if (no_limits \|\| cur<max_current)
317	g1.SetPoint(g1.GetN(), axis, hrz.alt);
318
319	if (no_limits \|\| 90-hrz.alt<max_zd)
320	g2.SetPoint(g2.GetN(), axis, cur);
321
322	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
323	g3.SetPoint(g3.GetN(), axis, ratio*cur/6.2);
324
325	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
326	g4.SetPoint(g4.GetN(), axis, angle);
327
328	if (cnt==0)
329	gm.SetPoint(gm.GetN(), axis, hrzm.alt);
330	}
331
332	if (cnt==0)
333	DrawClone(c1, gm);
334
335	DrawClone(c1, g1);
336	DrawClone(c2, g2);
337	DrawClone(c3, g3);
338	DrawClone(c4, g4);
339	}
340
341
342	// Save three plots
343	TCanvas c5;
344	c5.SetFillColor(kWhite);
345	c5.SetBorderMode(0);
346	c5.SetFrameBorderMode(0);
347	leg.Draw();
348
349	const string t = Time(jd).GetAsStr("%Y%m%d");
350
351	c1.SaveAs((t+"-ZenithDistance.eps").c_str());
352	c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
353	c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
354	c4.SaveAs((t+"-MoonDist.eps").c_str());
355	c5.SaveAs((t+"-Legend.eps").c_str());
356
357	c1.SaveAs((t+"-ZenithDistance.root").c_str());
358	c2.SaveAs((t+"-PredictedCurrent.root").c_str());
359	c3.SaveAs((t+"-RelativeThreshold.root").c_str());
360	c4.SaveAs((t+"-MoonDist.root").c_str());
361
362	c1.Print((t+".pdf(").c_str(), "pdf");
363	c2.Print((t+".pdf" ).c_str(), "pdf");
364	c3.Print((t+".pdf" ).c_str(), "pdf");
365	c4.Print((t+".pdf" ).c_str(), "pdf");
366	c5.Print((t+".pdf)").c_str(), "pdf");
367
368	return 0;
369	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: