Index: trunk/FACT++/src/makedata.cc =================================================================== --- trunk/FACT++/src/makedata.cc (revision 14784) +++ trunk/FACT++/src/makedata.cc (revision 14784) @@ -0,0 +1,218 @@ +#include +#include +#include +#include + +#include "Database.h" + +#include "Time.h" +#include "Configuration.h" + +using namespace std; + +// ------------------------------------------------------------------------ + +double Angle(double ra, double dec, double r, double d) +{ + const double theta0 = M_PI/2-d*M_PI/180; + const double phi0 = r*M_PI/12; + + const double theta1 = M_PI/2-dec*M_PI/180; + const double phi1 = ra*M_PI/12; + + const double x0 = sin(theta0) * cos(phi0); + const double y0 = sin(theta0) * sin(phi0); + const double z0 = cos(theta0); + + const double x1 = sin(theta1) * cos(phi1); + const double y1 = sin(theta1) * sin(phi1); + const double z1 = cos(theta1); + + double arg = x0*x1 + y0*y1 + z0*z1; + if(arg > 1.0) arg = 1.0; + if(arg < -1.0) arg = -1.0; + + return acos(arg) * 180/M_PI; +} + +// ======================================================================== +// ======================================================================== +// ======================================================================== + +void SetupConfiguration(Configuration &conf) +{ + po::options_description control("Smart FACT"); + control.add_options() + ("source-name", var(), "Source name") + ("date-time", var(), "SQL time (UTC)") + ("source-database", var(""), "Database link as in\n\tuser:password@server[:port]/database.") + ("max-current", var(75), "Maximum current to display in other plots.") + ("max-zd", var(75), "Maximum zenith distance to display in other plots") + ("no-limits", po_switch(), "Switch off limits in plots") + ; + + po::positional_options_description p; + p.add("source-name", 1); // The 1st positional options + p.add("date-time", 2); // The 2nd positional options + + conf.AddOptions(control); + conf.SetArgumentPositions(p); +} + +void PrintUsage() +{ + cout << + "makedata - The astronomy data listing\n" + "\n" +// "Calculates several plots for the sources in the database\n" +// "helpful or needed for scheduling. The Plot is always calculated\n" +// "for the night which starts at the same so. So no matter if\n" +// "you specify '1974-09-09 00:00:00' or '1974-09-09 21:56:00'\n" +// "the plots will refer to the night 1974-09-09/1974-09-10.\n" +// "The advantage is that specification of the date as in\n" +// "1974-09-09 is enough. Time axis starts and ends at nautical\n" +// "twilight which is 12deg below horizon.\n" +// "\n" + "Usage: makedata sql-datetime [--ra={ra} --dec={dec}]\n"; + cout << endl; +} + +int main(int argc, const char* argv[]) +{ + Configuration conf(argv[0]); + conf.SetPrintUsage(PrintUsage); + SetupConfiguration(conf); + + if (!conf.DoParse(argc, argv)) + return 127; + + if (!conf.Has("source-name")) + { + cout << "ERROR - --source-name missing." << endl; + return 1; + } + + // ------------------ Eval config --------------------- + + const double lon = -(17.+53./60+26.525/3600); + const double lat = 28.+45./60+42.462/3600; + + ln_lnlat_posn observer; + observer.lng = lon; + observer.lat = lat; + + Time time; + if (conf.Has("date-time")) + time.SetFromStr(conf.Get("date-time")); + + const string source_name = conf.Get("source-name"); + const double max_current = conf.Get("max-current"); + const double max_zd = conf.Get("max-zd"); + const double no_limits = conf.Get("no-limits"); + + // -12: nautical + ln_rst_time sun_set; // Sun set with the same date than th provided date + ln_rst_time sun_rise; // Sun rise on the following day + ln_get_solar_rst_horizon(time.JD()-0.5, &observer, -12, &sun_set); + ln_get_solar_rst_horizon(time.JD()+0.5, &observer, -12, &sun_rise); + + const double jd = floor(time.Mjd())+2400001; + const double mjd = floor(time.Mjd())+49718+0.5; + + + const double jd0 = fmod(sun_set.set, 1); // ~0.3 + const double jd1 = fmod(sun_rise.rise, 1); // ~0.8 + + cout << Time::iso << time << ", " << mjd-49718 << ", "; + cout << jd0 << ", "; + cout << jd1 << "\n"; + + const string fDatabase = conf.Get("source-database"); + + // ------------------ Precalc moon coord --------------------- + + vector> fMoonCoords; + + for (double h=0; h<1; h+=1./(24*12)) + { + ln_equ_posn moon; + ln_get_lunar_equ_coords_prec(jd+h, &moon, 0.01); + + const double disk = ln_get_lunar_disk(jd+h); + + fMoonCoords.push_back(make_pair(moon, disk)); + } + + // ------------- Get Sources from databasse --------------------- + + const mysqlpp::StoreQueryResult res = + Database(fDatabase).query("SELECT fRightAscension, fDeclination FROM source WHERE fSourceName='"+source_name+"'").store(); + + // ------------- Create canvases and frames --------------------- + + vector::const_iterator row=res.begin(); + + ln_equ_posn pos; + pos.ra = double((*row)[0])*15; + pos.dec = double((*row)[1]); + + // Loop over 24 hours + for (int i=0; i<24*12; i++) + { + const double h = double(i)/(24*12); + + // check if it is betwene sun-rise and sun-set + if (hjd1) + continue; + + // get local position of source + ln_hrz_posn hrz; + ln_get_hrz_from_equ(&pos, &observer, jd+h, &hrz); + + // Get moon properties and + ln_equ_posn moon = fMoonCoords[i].first; + const double disk = fMoonCoords[i].second; + + ln_hrz_posn hrzm; + ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrzm); + + // Distance between source and moon + const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec); + + // Current prediction + const double lc = angle*hrzm.alt*pow(disk, 6)/360/360; + const double cur = lc>0 ? 7.7+4942*lc : 7.7; + + // Relative energy threshold prediction + const double cs = cos((90+hrz.alt)*M_PI/180); + const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.; + + // Add points to curve + // const double axis = (mjd+h)*24*3600; + + Time t(mjd-49718); + t += boost::posix_time::minutes(i*5); + + cout << t << ", " << h << ", "; + + if (no_limits || cur