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