| 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 | // ========================================================================
|
|---|
| 45 | // ========================================================================
|
|---|
| 46 | // ========================================================================
|
|---|
| 47 |
|
|---|
| 48 | void 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 |
|
|---|
| 65 | void 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 |
|
|---|
| 74 | int 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 | }
|
|---|
