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

Last change on this file since 14476 was 14450, checked in by tbretz, 12 years ago
File size: 6.9 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 ;
57
58 po::positional_options_description p;
59 p.add("date-time", 1); // The first positional options
60
61 conf.AddOptions(control);
62 conf.SetArgumentPositions(p);
63}
64
65void PrintUsage()
66{
67 cout <<
68 "makeplots - The astronomy plotter\n"
69 "\n"
70 "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
71 cout << endl;
72}
73
74int main(int argc, const char* argv[])
75{
76 gROOT->SetBatch();
77
78 Configuration conf(argv[0]);
79 conf.SetPrintUsage(PrintUsage);
80 SetupConfiguration(conf);
81
82 if (!conf.DoParse(argc, argv))
83 return 127;
84
85 if (conf.Has("ra")^conf.Has("dec"))
86 {
87 cout << "ERROR - Either --ra or --dec missing." << endl;
88 return 1;
89 }
90
91 // ------------------ Eval config ---------------------
92
93 const double lon = -(17.+53./60+26.525/3600);
94 const double lat = 28.+45./60+42.462/3600;
95
96 ln_lnlat_posn observer;
97 observer.lng = lon;
98 observer.lat = lat;
99
100 Time time;
101 if (conf.Has("date-time"))
102 time.SetFromStr(conf.Get<string>("date-time"));
103
104 ln_rst_time sun_astronomical;
105 ln_get_solar_rst_horizon(time.JD(), &observer, -12, &sun_astronomical);
106
107 const double jd = floor(time.JD());
108 const double mjd = floor(time.Mjd())+49718-0.5;
109
110 const double jd0 = fmod(sun_astronomical.set, 1);
111 const double jd1 = fmod(sun_astronomical.rise, 1);
112
113 cout << "JD0=" << jd+jd0 << endl;
114 cout << "JD1=" << jd+jd1 << endl;
115
116 const string fDatabase = conf.Get<string>("source-database");
117
118 // ------------------ Precalc moon coord ---------------------
119
120 vector<pair<ln_equ_posn, double>> fMoonCoords;
121
122 for (double h=0; h<1; h+=1./(24*12))
123 {
124 ln_equ_posn moon;
125 ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01);
126
127 const double disk = ln_get_lunar_disk(jd+h);
128
129 fMoonCoords.push_back(make_pair(moon, disk));
130 }
131
132 // ------------- Get Sources from databasse ---------------------
133
134 const mysqlpp::StoreQueryResult res =
135 Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM source").store();
136
137 // ------------- Create canvases and frames ---------------------
138
139 TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
140 hframe.SetStats(kFALSE);
141 hframe.GetXaxis()->SetTimeFormat("%Hh%M'");
142 hframe.GetXaxis()->SetTitle("Time");
143 hframe.GetXaxis()->CenterTitle();
144 hframe.GetYaxis()->CenterTitle();
145 hframe.GetXaxis()->SetTimeDisplay(true);
146 hframe.GetXaxis()->SetLabelSize(0.025);
147
148 TCanvas c1;
149 gPad->SetRightMargin(0.01);
150 gPad->SetTopMargin(0.03);
151 gPad->SetGrid();
152 hframe.GetYaxis()->SetTitle("Altitude [deg]");
153 hframe.SetMinimum(15);
154 hframe.SetMaximum(90);
155 hframe.DrawCopy();
156
157 TCanvas c2;
158 gPad->SetRightMargin(0.01);
159 gPad->SetTopMargin(0.03);
160 gPad->SetGrid();
161 hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
162 hframe.SetMinimum(0);
163 hframe.SetMaximum(100);
164 hframe.DrawCopy();
165
166 Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
167 Int_t style[] = { kSolid, kDashed, kDotted };
168
169 TLegend leg(0, 0, 1, 1);
170
171 // ------------- Loop over sources ---------------------
172
173 Int_t cnt=0;
174 for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
175 {
176 // Eval row
177 const string name = (*v)[0].c_str();
178
179 ln_equ_posn pos;
180 pos.ra = double((*v)[1])*15;
181 pos.dec = double((*v)[2]);
182
183 // Create graphs
184 TGraph g1, g2;
185 g1.SetName(name.data());
186 g2.SetName(name.data());
187 g1.SetLineWidth(2);
188 g2.SetLineWidth(2);
189 g1.SetLineStyle(style[cnt/6]);
190 g2.SetLineStyle(style[cnt/6]);
191 g1.SetLineColor(color[cnt%6]);
192 g2.SetLineColor(color[cnt%6]);
193
194 // Loop over 24 hours
195 int i=0;
196 for (double h=0; h<1; h+=1./(24*12), i++)
197 {
198 // check if it is betwene sun-rise and sun-set
199 if (h<jd0 || h>jd1)
200 continue;
201
202 // get local position of source
203 ln_hrz_posn hrz;
204 ln_get_hrz_from_equ(&pos, &observer, jd+h, &hrz);
205
206 // Check if source is visible
207 if (hrz.alt<15)
208 continue;
209
210 // Add point to curve
211 const double axis = (mjd+h)*24*3600;
212 g1.SetPoint(g1.GetN(), axis, hrz.alt);
213
214 // Get moon properties and estimate current
215 ln_equ_posn moon = fMoonCoords[i].first;
216 const double disk = fMoonCoords[i].second;
217
218 ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrz);
219
220 const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec);
221
222 const double lc = angle*hrz.alt*pow(disk, 6)/360/360;
223
224 // Add point to curve
225 g2.SetPoint(g2.GetN(), axis, lc>0 ? 7.7+4942*lc : 7.7);
226 }
227
228 // Add graphs to canvases and add corresponding entry to legend
229 TGraph *g = 0;
230
231 c1.cd();
232 g = (TGraph*)g1.DrawClone("C");
233 g->SetBit(kCanDelete);
234
235 c2.cd();
236 g = (TGraph*)g2.DrawClone("C");
237 g->SetBit(kCanDelete);
238
239 leg.AddEntry(g, name.data(), "l");
240 }
241
242 // Save three plots
243 TCanvas c3;
244 leg.Draw();
245
246 c1.SaveAs("test1.eps");
247 c2.SaveAs("test2.eps");
248 c3.SaveAs("legend.eps");
249
250 c1.SaveAs("test1.root");
251 c2.SaveAs("test2.root");
252
253 return 0;
254}
Note: See TracBrowser for help on using the repository browser.