source: trunk/FACT++/src/makedata.cc@ 15034

#include <libnova/solar.h>
#include <libnova/lunar.h>
#include <libnova/rise_set.h>
#include <libnova/transform.h>

#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<string>(), "Source name")
        ("date-time", var<string>(), "SQL time (UTC)")
        ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
        ("max-current", var<double>(75), "Maximum current to display in other plots.")
        ("max-zd", var<double>(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<string>("date-time"));

    const double max_current = conf.Get<double>("max-current");
    const double max_zd      = conf.Get<double>("max-zd");
    const double no_limits   = conf.Get<bool>("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";

    if (!conf.Has("source-name"))
        return 1;

    const string source_name = conf.Get<string>("source-name");
    const string fDatabase   = conf.Get<string>("source-database");

    // ------------------ Precalc moon coord ---------------------

    vector<pair<ln_equ_posn, double>> 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<mysqlpp::Row>::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 (h<jd0 || h>jd1)
            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 cang = sin(angle   *M_PI/180);
        const double calt = sin(hrzm.alt*M_PI/180);

        const double lc = cang==0 ? -1 : calt*pow(disk, 3)/sqrt(cang);
        const double cur = lc>0 ? 8.1+94.6*lc : 8.1;

        // 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<max_current)
            cout << hrz.alt;
        cout << ", ";

        if (no_limits || 90-hrz.alt<max_zd)
            cout << cur;
        cout << ", ";

        if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
            cout << ratio*cur/8.1;
        cout << ", ";

        if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
            cout << angle;
        cout << "\n";
    }

    cout << flush;

    return 0;
}