| 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>(75), "Maximum current to display in other plots.")
|
|---|
| 52 | ("max-zd", var<double>(75), "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 [--ra={ra} --dec={dec}]\n";
|
|---|
| 78 | cout << endl;
|
|---|
| 79 | }
|
|---|
| 80 |
|
|---|
| 81 | int main(int argc, const char* argv[])
|
|---|
| 82 | {
|
|---|
| 83 | gROOT->SetBatch();
|
|---|
| 84 |
|
|---|
| 85 | Configuration conf(argv[0]);
|
|---|
| 86 | conf.SetPrintUsage(PrintUsage);
|
|---|
| 87 | SetupConfiguration(conf);
|
|---|
| 88 |
|
|---|
| 89 | if (!conf.DoParse(argc, argv))
|
|---|
| 90 | return 127;
|
|---|
| 91 |
|
|---|
| 92 | if (conf.Has("ra")^conf.Has("dec"))
|
|---|
| 93 | {
|
|---|
| 94 | cout << "ERROR - Either --ra or --dec missing." << endl;
|
|---|
| 95 | return 1;
|
|---|
| 96 | }
|
|---|
| 97 |
|
|---|
| 98 | // ------------------ Eval config ---------------------
|
|---|
| 99 |
|
|---|
| 100 | Time time;
|
|---|
| 101 | if (conf.Has("date-time"))
|
|---|
| 102 | time.SetFromStr(conf.Get<string>("date-time"));
|
|---|
| 103 |
|
|---|
| 104 | const double max_current = conf.Get<double>("max-current");
|
|---|
| 105 | const double max_zd = conf.Get<double>("max-zd");
|
|---|
| 106 | const double no_limits = conf.Get<bool>("no-limits");
|
|---|
| 107 |
|
|---|
| 108 | // -12: nautical
|
|---|
| 109 | // Sun set with the same date than th provided date
|
|---|
| 110 | // Sun rise on the following day
|
|---|
| 111 | const RstTime sun_set = GetSolarRst(time.JD()-0.5, -12);
|
|---|
| 112 | const RstTime sun_rise = GetSolarRst(time.JD()+0.5, -12);
|
|---|
| 113 |
|
|---|
| 114 | const double jd = floor(time.Mjd())+2400001;
|
|---|
| 115 | const double mjd = floor(time.Mjd())+49718+0.5; // root time
|
|---|
| 116 |
|
|---|
| 117 | cout << "Time: " << time << endl;
|
|---|
| 118 | cout << "Base: " << Time(jd-0.5) << endl;
|
|---|
| 119 | cout << "Set: " << Time(sun_set.set) << endl;
|
|---|
| 120 | cout << "Rise: " << Time(sun_rise.rise) << endl;
|
|---|
| 121 |
|
|---|
| 122 | const double jd0 = fmod(sun_set.set, 1); // ~0.3
|
|---|
| 123 | const double jd1 = fmod(sun_rise.rise, 1); // ~0.8
|
|---|
| 124 |
|
|---|
| 125 | const string fDatabase = conf.Get<string>("source-database");
|
|---|
| 126 |
|
|---|
| 127 | // ------------------ Precalc moon coord ---------------------
|
|---|
| 128 |
|
|---|
| 129 | vector<pair<EquPosn, double>> fMoonCoords;
|
|---|
| 130 |
|
|---|
| 131 | for (double h=0; h<1; h+=1./(24*12))
|
|---|
| 132 | {
|
|---|
| 133 | const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
|
|---|
| 134 | const double disk = GetLunarDisk(jd+h);
|
|---|
| 135 |
|
|---|
| 136 | fMoonCoords.emplace_back(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 | // It is important to use an offset which is larger than
|
|---|
| 147 | // 1970-01-01 00:00:00. This one will not work if your
|
|---|
| 148 | // local time zone is positive!
|
|---|
| 149 | TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
|
|---|
| 150 | hframe.SetStats(kFALSE);
|
|---|
| 151 | hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
|
|---|
| 152 | hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
|
|---|
| 153 | hframe.GetXaxis()->CenterTitle();
|
|---|
| 154 | hframe.GetYaxis()->CenterTitle();
|
|---|
| 155 | hframe.GetXaxis()->SetTimeDisplay(true);
|
|---|
| 156 | hframe.GetYaxis()->SetTitleSize(0.040);
|
|---|
| 157 | hframe.GetXaxis()->SetTitleSize(0.040);
|
|---|
| 158 | hframe.GetXaxis()->SetTitleOffset(1.1);
|
|---|
| 159 | hframe.GetYaxis()->SetLabelSize(0.040);
|
|---|
| 160 | hframe.GetXaxis()->SetLabelSize(0.040);
|
|---|
| 161 |
|
|---|
| 162 | TCanvas c1;
|
|---|
| 163 | c1.SetFillColor(kWhite);
|
|---|
| 164 | c1.SetBorderMode(0);
|
|---|
| 165 | c1.SetFrameBorderMode(0);
|
|---|
| 166 | c1.SetLeftMargin(0.085);
|
|---|
| 167 | c1.SetRightMargin(0.01);
|
|---|
| 168 | c1.SetTopMargin(0.03);
|
|---|
| 169 | c1.SetGrid();
|
|---|
| 170 | hframe.GetYaxis()->SetTitle("Altitude [deg]");
|
|---|
| 171 | hframe.SetMinimum(15);
|
|---|
| 172 | hframe.SetMaximum(90);
|
|---|
| 173 | hframe.DrawCopy();
|
|---|
| 174 |
|
|---|
| 175 | TCanvas c2;
|
|---|
| 176 | c2.SetFillColor(kWhite);
|
|---|
| 177 | c2.SetBorderMode(0);
|
|---|
| 178 | c2.SetFrameBorderMode(0);
|
|---|
| 179 | c2.SetLeftMargin(0.085);
|
|---|
| 180 | c2.SetRightMargin(0.01);
|
|---|
| 181 | c2.SetTopMargin(0.03);
|
|---|
| 182 | c2.SetGrid();
|
|---|
| 183 | hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
|
|---|
| 184 | hframe.SetMinimum(0);
|
|---|
| 185 | hframe.SetMaximum(100);
|
|---|
| 186 | hframe.DrawCopy();
|
|---|
| 187 |
|
|---|
| 188 | TCanvas c3;
|
|---|
| 189 | c3.SetFillColor(kWhite);
|
|---|
| 190 | c3.SetBorderMode(0);
|
|---|
| 191 | c3.SetFrameBorderMode(0);
|
|---|
| 192 | c3.SetLeftMargin(0.085);
|
|---|
| 193 | c3.SetRightMargin(0.01);
|
|---|
| 194 | c3.SetTopMargin(0.03);
|
|---|
| 195 | c3.SetGrid();
|
|---|
| 196 | c3.SetLogy();
|
|---|
| 197 | hframe.GetYaxis()->SetTitle("Estimated relative threshold");
|
|---|
| 198 | hframe.SetMinimum(0.9);
|
|---|
| 199 | hframe.SetMaximum(180);
|
|---|
| 200 | hframe.DrawCopy();
|
|---|
| 201 |
|
|---|
| 202 | TCanvas c4;
|
|---|
| 203 | c4.SetFillColor(kWhite);
|
|---|
| 204 | c4.SetBorderMode(0);
|
|---|
| 205 | c4.SetFrameBorderMode(0);
|
|---|
| 206 | c4.SetLeftMargin(0.085);
|
|---|
| 207 | c4.SetRightMargin(0.01);
|
|---|
| 208 | c4.SetTopMargin(0.03);
|
|---|
| 209 | c4.SetGrid();
|
|---|
| 210 | hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
|
|---|
| 211 | hframe.SetMinimum(0);
|
|---|
| 212 | hframe.SetMaximum(180);
|
|---|
| 213 | hframe.DrawCopy();
|
|---|
| 214 |
|
|---|
| 215 | Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
|
|---|
| 216 | Int_t style[] = { kSolid, kDashed, kDotted };
|
|---|
| 217 |
|
|---|
| 218 | TLegend leg(0, 0, 1, 1);
|
|---|
| 219 |
|
|---|
| 220 | // ------------- Loop over sources ---------------------
|
|---|
| 221 |
|
|---|
| 222 | Int_t cnt=0;
|
|---|
| 223 | for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
|
|---|
| 224 | {
|
|---|
| 225 | // Eval row
|
|---|
| 226 | const string name = (*v)[0].c_str();
|
|---|
| 227 |
|
|---|
| 228 | EquPosn pos;
|
|---|
| 229 | pos.ra = double((*v)[1])*15;
|
|---|
| 230 | pos.dec = double((*v)[2]);
|
|---|
| 231 |
|
|---|
| 232 | // Create graphs
|
|---|
| 233 | TGraph g1, g2, g3, g4, gr, gm;
|
|---|
| 234 | g1.SetName(name.data());
|
|---|
| 235 | g2.SetName(name.data());
|
|---|
| 236 | g3.SetName(name.data());
|
|---|
| 237 | g4.SetName(name.data());
|
|---|
| 238 | g1.SetLineWidth(2);
|
|---|
| 239 | g2.SetLineWidth(2);
|
|---|
| 240 | g3.SetLineWidth(2);
|
|---|
| 241 | g4.SetLineWidth(2);
|
|---|
| 242 | gm.SetLineWidth(1);
|
|---|
| 243 | g1.SetLineStyle(style[cnt/6]);
|
|---|
| 244 | g2.SetLineStyle(style[cnt/6]);
|
|---|
| 245 | g3.SetLineStyle(style[cnt/6]);
|
|---|
| 246 | g4.SetLineStyle(style[cnt/6]);
|
|---|
| 247 | g1.SetLineColor(color[cnt%6]);
|
|---|
| 248 | g2.SetLineColor(color[cnt%6]);
|
|---|
| 249 | g3.SetLineColor(color[cnt%6]);
|
|---|
| 250 | g4.SetLineColor(color[cnt%6]);
|
|---|
| 251 | gm.SetLineColor(kYellow);
|
|---|
| 252 |
|
|---|
| 253 | if (cnt==0)
|
|---|
| 254 | leg.AddEntry(gm.Clone(), "Moon", "l");
|
|---|
| 255 | leg.AddEntry(g1.Clone(), name.data(), "l");
|
|---|
| 256 |
|
|---|
| 257 | // Loop over 24 hours
|
|---|
| 258 | int i=0;
|
|---|
| 259 | for (double h=0; h<1; h+=1./(24*12), i++)
|
|---|
| 260 | {
|
|---|
| 261 | // check if it is betwene sun-rise and sun-set
|
|---|
| 262 | if (h<jd0 || h>jd1)
|
|---|
| 263 | continue;
|
|---|
| 264 |
|
|---|
| 265 | // get local position of source
|
|---|
| 266 | const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
|
|---|
| 267 |
|
|---|
| 268 | // Get moon properties and
|
|---|
| 269 | const EquPosn moon = fMoonCoords[i].first;
|
|---|
| 270 | const double disk = fMoonCoords[i].second;
|
|---|
| 271 | const HrzPosn hrzm = GetHrzFromEqu(moon, jd+h);
|
|---|
| 272 |
|
|---|
| 273 | if (v==res.begin())
|
|---|
| 274 | cout << Time(jd+h) <<" " << 90-hrzm.alt << endl;
|
|---|
| 275 |
|
|---|
| 276 | // Distance between source and moon
|
|---|
| 277 | const double angle = GetAngularSeparation(moon, pos);
|
|---|
| 278 |
|
|---|
| 279 | // Distance between earth and moon relative to major semi-axis
|
|---|
| 280 | const double dist = GetLunarEarthDist(jd+h)/384400;
|
|---|
| 281 |
|
|---|
| 282 | // Current prediction
|
|---|
| 283 | const double cang = sin(angle *M_PI/180);
|
|---|
| 284 | const double calt = sin(hrzm.alt*M_PI/180);
|
|---|
| 285 |
|
|---|
| 286 | const double lc = cang==0 ? -1 : calt*pow(disk, 2.2)*pow(dist, -2);
|
|---|
| 287 | const double cur = lc>0 ? 4+103*lc : 4;
|
|---|
| 288 |
|
|---|
| 289 | // Relative energy threshold prediction
|
|---|
| 290 | //const double cs = cos((90+hrz.alt)*M_PI/180);
|
|---|
| 291 | //const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;
|
|---|
| 292 | const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
|
|---|
| 293 |
|
|---|
| 294 | // Add points to curve
|
|---|
| 295 | const double axis = (mjd+h)*24*3600;
|
|---|
| 296 |
|
|---|
| 297 | // If there is a gap of more than one bin, start a new curve
|
|---|
| 298 | CheckForGap(c1, g1, axis);
|
|---|
| 299 | CheckForGap(c1, gm, axis);
|
|---|
| 300 | CheckForGap(c2, g2, axis);
|
|---|
| 301 | CheckForGap(c3, g3, axis);
|
|---|
| 302 | CheckForGap(c4, g4, axis);
|
|---|
| 303 |
|
|---|
| 304 | // Add data
|
|---|
| 305 | if (no_limits || cur<max_current)
|
|---|
| 306 | g1.SetPoint(g1.GetN(), axis, hrz.alt);
|
|---|
| 307 |
|
|---|
| 308 | if (no_limits || 90-hrz.alt<max_zd)
|
|---|
| 309 | g2.SetPoint(g2.GetN(), axis, cur);
|
|---|
| 310 |
|
|---|
| 311 | if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
|
|---|
| 312 | g3.SetPoint(g3.GetN(), axis, ratio*cur/8.1);
|
|---|
| 313 |
|
|---|
| 314 | if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
|
|---|
| 315 | g4.SetPoint(g4.GetN(), axis, angle);
|
|---|
| 316 |
|
|---|
| 317 | if (cnt==0)
|
|---|
| 318 | gm.SetPoint(gm.GetN(), axis, hrzm.alt);
|
|---|
| 319 | }
|
|---|
| 320 |
|
|---|
| 321 | if (cnt==0)
|
|---|
| 322 | DrawClone(c1, gm);
|
|---|
| 323 |
|
|---|
| 324 | DrawClone(c1, g1);
|
|---|
| 325 | DrawClone(c2, g2);
|
|---|
| 326 | DrawClone(c3, g3);
|
|---|
| 327 | DrawClone(c4, g4);
|
|---|
| 328 | }
|
|---|
| 329 |
|
|---|
| 330 |
|
|---|
| 331 | // Save three plots
|
|---|
| 332 | TCanvas c5;
|
|---|
| 333 | c5.SetFillColor(kWhite);
|
|---|
| 334 | c5.SetBorderMode(0);
|
|---|
| 335 | c5.SetFrameBorderMode(0);
|
|---|
| 336 | leg.Draw();
|
|---|
| 337 |
|
|---|
| 338 | const string t = Time(jd).GetAsStr("%Y%m%d");
|
|---|
| 339 |
|
|---|
| 340 | c1.SaveAs((t+"-ZenithDistance.eps").c_str());
|
|---|
| 341 | c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
|
|---|
| 342 | c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
|
|---|
| 343 | c4.SaveAs((t+"-MoonDist.eps").c_str());
|
|---|
| 344 | c5.SaveAs((t+"-Legend.eps").c_str());
|
|---|
| 345 |
|
|---|
| 346 | c1.SaveAs((t+"-ZenithDistance.root").c_str());
|
|---|
| 347 | c2.SaveAs((t+"-PredictedCurrent.root").c_str());
|
|---|
| 348 | c3.SaveAs((t+"-RelativeThreshold.root").c_str());
|
|---|
| 349 | c4.SaveAs((t+"-MoonDist.root").c_str());
|
|---|
| 350 |
|
|---|
| 351 | return 0;
|
|---|
| 352 | }
|
|---|
