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