| 1 | #include "Prediction.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("Makeplots");
|
|---|
| 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 |
|
|---|
| 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 sunset = sun_set.set;
|
|---|
| 123 | const double sunrise = sun_rise.rise;
|
|---|
| 124 |
|
|---|
| 125 | const string fDatabase = conf.Get<string>("source-database");
|
|---|
| 126 |
|
|---|
| 127 | // ------------- Get Sources from databasse ---------------------
|
|---|
| 128 |
|
|---|
| 129 | const mysqlpp::StoreQueryResult res =
|
|---|
| 130 | Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1").store();
|
|---|
| 131 |
|
|---|
| 132 | // ------------- Create canvases and frames ---------------------
|
|---|
| 133 |
|
|---|
| 134 | // It is important to use an offset which is larger than
|
|---|
| 135 | // 1970-01-01 00:00:00. This one will not work if your
|
|---|
| 136 | // local time zone is positive!
|
|---|
| 137 | TH1S hframe("", "", 1, Time(sunset).Mjd()*24*3600, Time(sunrise).Mjd()*24*3600);
|
|---|
| 138 | hframe.SetStats(kFALSE);
|
|---|
| 139 | hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
|
|---|
| 140 | hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
|
|---|
| 141 | hframe.GetXaxis()->CenterTitle();
|
|---|
| 142 | hframe.GetYaxis()->CenterTitle();
|
|---|
| 143 | hframe.GetXaxis()->SetTimeDisplay(true);
|
|---|
| 144 | hframe.GetYaxis()->SetTitleSize(0.040);
|
|---|
| 145 | hframe.GetXaxis()->SetTitleSize(0.040);
|
|---|
| 146 | hframe.GetXaxis()->SetTitleOffset(1.1);
|
|---|
| 147 | hframe.GetYaxis()->SetLabelSize(0.040);
|
|---|
| 148 | hframe.GetXaxis()->SetLabelSize(0.040);
|
|---|
| 149 |
|
|---|
| 150 | TCanvas c1;
|
|---|
| 151 | c1.SetFillColor(kWhite);
|
|---|
| 152 | c1.SetBorderMode(0);
|
|---|
| 153 | c1.SetFrameBorderMode(0);
|
|---|
| 154 | c1.SetLeftMargin(0.085);
|
|---|
| 155 | c1.SetRightMargin(0.01);
|
|---|
| 156 | c1.SetTopMargin(0.03);
|
|---|
| 157 | c1.SetGrid();
|
|---|
| 158 | hframe.GetYaxis()->SetTitle("Altitude [deg]");
|
|---|
| 159 | hframe.SetMinimum(15);
|
|---|
| 160 | hframe.SetMaximum(90);
|
|---|
| 161 | hframe.DrawCopy();
|
|---|
| 162 |
|
|---|
| 163 | TCanvas c2;
|
|---|
| 164 | c2.SetFillColor(kWhite);
|
|---|
| 165 | c2.SetBorderMode(0);
|
|---|
| 166 | c2.SetFrameBorderMode(0);
|
|---|
| 167 | c2.SetLeftMargin(0.085);
|
|---|
| 168 | c2.SetRightMargin(0.01);
|
|---|
| 169 | c2.SetTopMargin(0.03);
|
|---|
| 170 | c2.SetGrid();
|
|---|
| 171 | hframe.GetYaxis()->SetTitle("Predicted Current [#muA]");
|
|---|
| 172 | hframe.SetMinimum(0);
|
|---|
| 173 | hframe.SetMaximum(100);
|
|---|
| 174 | hframe.DrawCopy();
|
|---|
| 175 |
|
|---|
| 176 | TCanvas c3;
|
|---|
| 177 | c3.SetFillColor(kWhite);
|
|---|
| 178 | c3.SetBorderMode(0);
|
|---|
| 179 | c3.SetFrameBorderMode(0);
|
|---|
| 180 | c3.SetLeftMargin(0.085);
|
|---|
| 181 | c3.SetRightMargin(0.01);
|
|---|
| 182 | c3.SetTopMargin(0.03);
|
|---|
| 183 | c3.SetGrid();
|
|---|
| 184 | c3.SetLogy();
|
|---|
| 185 | hframe.GetYaxis()->SetTitle("Estimated relative threshold");
|
|---|
| 186 | hframe.GetYaxis()->SetMoreLogLabels();
|
|---|
| 187 | hframe.SetMinimum(0.9);
|
|---|
| 188 | hframe.SetMaximum(11);
|
|---|
| 189 | hframe.DrawCopy();
|
|---|
| 190 |
|
|---|
| 191 | TCanvas c4;
|
|---|
| 192 | c4.SetFillColor(kWhite);
|
|---|
| 193 | c4.SetBorderMode(0);
|
|---|
| 194 | c4.SetFrameBorderMode(0);
|
|---|
| 195 | c4.SetLeftMargin(0.085);
|
|---|
| 196 | c4.SetRightMargin(0.01);
|
|---|
| 197 | c4.SetTopMargin(0.03);
|
|---|
| 198 | c4.SetGrid();
|
|---|
| 199 | hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
|
|---|
| 200 | hframe.SetMinimum(0);
|
|---|
| 201 | hframe.SetMaximum(180);
|
|---|
| 202 | hframe.DrawCopy();
|
|---|
| 203 |
|
|---|
| 204 | Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
|
|---|
| 205 | Int_t style[] = { kSolid, kDashed, kDotted };
|
|---|
| 206 |
|
|---|
| 207 | TLegend leg(0, 0, 1, 1);
|
|---|
| 208 |
|
|---|
| 209 | // ------------- Loop over sources ---------------------
|
|---|
| 210 |
|
|---|
| 211 | Int_t cnt=0;
|
|---|
| 212 | for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
|
|---|
| 213 | {
|
|---|
| 214 | // Eval row
|
|---|
| 215 | const string name = (*v)[0].c_str();
|
|---|
| 216 |
|
|---|
| 217 | EquPosn pos;
|
|---|
| 218 | pos.ra = double((*v)[1])*15;
|
|---|
| 219 | pos.dec = double((*v)[2]);
|
|---|
| 220 |
|
|---|
| 221 | // Create graphs
|
|---|
| 222 | TGraph g1, g2, g3, g4, gr, gm;
|
|---|
| 223 | g1.SetName(name.data());
|
|---|
| 224 | g2.SetName(name.data());
|
|---|
| 225 | g3.SetName(name.data());
|
|---|
| 226 | g4.SetName(name.data());
|
|---|
| 227 | g1.SetLineWidth(2);
|
|---|
| 228 | g2.SetLineWidth(2);
|
|---|
| 229 | g3.SetLineWidth(2);
|
|---|
| 230 | g4.SetLineWidth(2);
|
|---|
| 231 | gm.SetLineWidth(1);
|
|---|
| 232 | g1.SetLineStyle(style[cnt/6]);
|
|---|
| 233 | g2.SetLineStyle(style[cnt/6]);
|
|---|
| 234 | g3.SetLineStyle(style[cnt/6]);
|
|---|
| 235 | g4.SetLineStyle(style[cnt/6]);
|
|---|
| 236 | g1.SetLineColor(color[cnt%6]);
|
|---|
| 237 | g2.SetLineColor(color[cnt%6]);
|
|---|
| 238 | g3.SetLineColor(color[cnt%6]);
|
|---|
| 239 | g4.SetLineColor(color[cnt%6]);
|
|---|
| 240 | gm.SetLineColor(kYellow);
|
|---|
| 241 |
|
|---|
| 242 | if (cnt==0)
|
|---|
| 243 | leg.AddEntry(gm.Clone(), "Moon", "l");
|
|---|
| 244 | leg.AddEntry(g1.Clone(), name.data(), "l");
|
|---|
| 245 |
|
|---|
| 246 | // Loop over 24 hours
|
|---|
| 247 | for (double h=0; h<1; h+=1./(24*12))
|
|---|
| 248 | {
|
|---|
| 249 | const SolarObjects so(jd+h);
|
|---|
| 250 |
|
|---|
| 251 | // get local position of source
|
|---|
| 252 | const HrzPosn hrz = GetHrzFromEqu(pos, so.fJD);
|
|---|
| 253 |
|
|---|
| 254 | if (v==res.begin())
|
|---|
| 255 | cout << Time(so.fJD) <<" " << 90-so.fMoonHrz.alt << endl;
|
|---|
| 256 |
|
|---|
| 257 | const double cur = FACT::PredictI(so, pos);
|
|---|
| 258 |
|
|---|
| 259 | // Relative energy threshold prediction
|
|---|
| 260 | const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
|
|---|
| 261 |
|
|---|
| 262 | // Add points to curve
|
|---|
| 263 | const double axis = Time(so.fJD).Mjd()*24*3600;
|
|---|
| 264 |
|
|---|
| 265 | // If there is a gap of more than one bin, start a new curve
|
|---|
| 266 | CheckForGap(c1, g1, axis);
|
|---|
| 267 | CheckForGap(c1, gm, axis);
|
|---|
| 268 | CheckForGap(c2, g2, axis);
|
|---|
| 269 | CheckForGap(c3, g3, axis);
|
|---|
| 270 | CheckForGap(c4, g4, axis);
|
|---|
| 271 |
|
|---|
| 272 | // Add data
|
|---|
| 273 | if (no_limits || cur<max_current)
|
|---|
| 274 | g1.SetPoint(g1.GetN(), axis, hrz.alt);
|
|---|
| 275 |
|
|---|
| 276 | if (no_limits || 90-hrz.alt<max_zd)
|
|---|
| 277 | g2.SetPoint(g2.GetN(), axis, cur);
|
|---|
| 278 |
|
|---|
| 279 | if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
|
|---|
| 280 | g3.SetPoint(g3.GetN(), axis, ratio*pow(cur/6.2, 0.394));
|
|---|
| 281 |
|
|---|
| 282 | if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
|
|---|
| 283 | {
|
|---|
| 284 | const double angle = GetAngularSeparation(so.fMoonEqu, pos);
|
|---|
| 285 | g4.SetPoint(g4.GetN(), axis, angle);
|
|---|
| 286 | }
|
|---|
| 287 |
|
|---|
| 288 | if (cnt==0)
|
|---|
| 289 | gm.SetPoint(gm.GetN(), axis, so.fMoonHrz.alt);
|
|---|
| 290 | }
|
|---|
| 291 |
|
|---|
| 292 | if (cnt==0)
|
|---|
| 293 | DrawClone(c1, gm);
|
|---|
| 294 |
|
|---|
| 295 | DrawClone(c1, g1);
|
|---|
| 296 | DrawClone(c2, g2);
|
|---|
| 297 | DrawClone(c3, g3);
|
|---|
| 298 | DrawClone(c4, g4);
|
|---|
| 299 | }
|
|---|
| 300 |
|
|---|
| 301 |
|
|---|
| 302 | // Save three plots
|
|---|
| 303 | TCanvas c5;
|
|---|
| 304 | c5.SetFillColor(kWhite);
|
|---|
| 305 | c5.SetBorderMode(0);
|
|---|
| 306 | c5.SetFrameBorderMode(0);
|
|---|
| 307 | leg.Draw();
|
|---|
| 308 |
|
|---|
| 309 | const string t = Time(jd).GetAsStr("%Y%m%d");
|
|---|
| 310 |
|
|---|
| 311 | c1.SaveAs((t+"-ZenithDistance.eps").c_str());
|
|---|
| 312 | c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
|
|---|
| 313 | c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
|
|---|
| 314 | c4.SaveAs((t+"-MoonDist.eps").c_str());
|
|---|
| 315 | c5.SaveAs((t+"-Legend.eps").c_str());
|
|---|
| 316 |
|
|---|
| 317 | c1.SaveAs((t+"-ZenithDistance.root").c_str());
|
|---|
| 318 | c2.SaveAs((t+"-PredictedCurrent.root").c_str());
|
|---|
| 319 | c3.SaveAs((t+"-RelativeThreshold.root").c_str());
|
|---|
| 320 | c4.SaveAs((t+"-MoonDist.root").c_str());
|
|---|
| 321 |
|
|---|
| 322 | c1.Print((t+".pdf(").c_str(), "pdf");
|
|---|
| 323 | c2.Print((t+".pdf" ).c_str(), "pdf");
|
|---|
| 324 | c3.Print((t+".pdf" ).c_str(), "pdf");
|
|---|
| 325 | c4.Print((t+".pdf" ).c_str(), "pdf");
|
|---|
| 326 | c5.Print((t+".pdf)").c_str(), "pdf");
|
|---|
| 327 |
|
|---|
| 328 | return 0;
|
|---|
| 329 | }
|
|---|
