Changeset 19684 for trunk/FACT++


Ignore:
Timestamp:
09/27/19 16:27:20 (5 years ago)
Author:
Daniela Dorner
Message:
added new algorithms (visiblity ratio, time-shift) and featues (plotting, zd-histogram)
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/FACT++/src/makeschedule.cc

    r19139 r19684  
    1919    control.add_options()
    2020        ("date", var<string>(), "SQL time (UTC), e.g. '2016-12-24' (equiv. '2016-12-24 12:00:00' to '2016-12-25 11:59:59')")
     21        ("svg", var<string>(), "File name to write a SVG file. Use together with --loop")
     22        ("hist", var<string>(), "File name to write a histogram file. Use together with --loop; default: false")
     23        ("loop", var<uint16_t>(uint16_t(0)), "Number of days to loop over (0:default, >=1:turns off database entry)")
    2124        ("source-database", var<string>()->required(), "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
    2225        ("schedule-database", var<string>(), "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
    2326        ("max-current", var<double>(90), "Global maximum current limit in uA")
    2427        ("max-zd", var<double>(75), "Global zenith distance limit in degree")
     28        ("min-moon", var<double>(5), "Minimum angular separation to moon in degree")
     29        ("exponent-visibility-ratio", var<double>(1), "Global exponent for visibility ratio")
     30        ("exponent-time-shift", var<double>(1), "Global exponent for time shift visibility ratio")
    2531        ("source", vars<string>(), "List of all TeV sources to be included, names according to the database")
    2632        ("setup.*", var<double>(), "Setup for the sources to be observed")
     
    2935        ("data-taking.start", var<double>(-12), "Begin of data-taking in degree of sun below horizon")
    3036        ("data-taking.end", var<double>(-13.75), "End of data-taking in degree of sun below horizon")
     37        ("use-visibility-ratio", po_bool(), "Use visibility ratio for scaling of threshold)")
     38        ("use-lst-shadow", po_bool(), "Respect shadow from LST for scheduling")
     39        ("print-hist", var<bool>(), "Print histogram to console; default: false)")
    3140        ("enter-schedule-into-database", var<bool>(), "Enter schedule into database (required schedule-database, false: dry-run)")
    3241        ;
     
    3443    po::positional_options_description p;
    3544    p.add("date", 1); // The first positional options
     45    p.add("loop", 1); // The first positional options
    3646
    3747    conf.AddOptions(control);
     
    4454        "makeschedule - Creates an automatic schedule\n"
    4555        "\n"
    46         "Usage: makeschedule [yyyy-mm-dd]\n";
     56        "Usage: makeschedule [yyyy-mm-dd [N]]\n";
    4757    cout << endl;
    4858}
     
    5060void PrintHelp()
    5161{
    52 #ifdef HAVE_ROOT
    5362    cout <<
    5463        "First for each minute of the night a list is calculated of all "
    55         "selected sources fulfiling all global and all source specific "
     64        "selected sources fulfilling all global and all source specific "
    5665        "constraints, e.g. on the zenith distance or the current.\n"
    5766        "\n"
     
    105114        "\n";
    106115    cout << endl;
    107 #endif
    108116}
    109117
     
    141149    static double max_current;
    142150    static double max_zd;
    143 
    144     // Source descrition
     151    static double min_moon;
     152    static double vis_exponent;
     153    static double ts_exponent;
     154    static bool   vis_ratio;
     155    static bool   lst_shadow;
     156
     157    // Source description
    145158    string name;
    146159    uint16_t key;
     
    151164    MyDouble maxcurrent;
    152165    double penalty;
     166    double time_shift;
     167    double visexponent;
     168    double tsexponent;
     169
     170    // Derived values
     171    double visratio;
     172    double timeshift;
    153173
    154174    // Possible observation time
     
    168188    //bool IsSpecial() const { return threshold==std::numeric_limits<double>::max(); }
    169189
     190    //distribution of zenith distance
     191    int zddistr[90];
     192
    170193    double zd(const double &jd) const
    171194    {
     
    173196    }
    174197
     198    /*
     199     // See https://www.fact-project.org/logbook/showthread.php?tid=6601
     200
     201     75 -135    outmost knot left side
     202     70 -134    tower left side
     203     65 -128    dish left side
     204     58 -125    intermediate position of dish on left side
     205     52 -118    intermediate position of dish on left side
     206     50 ~ -110  top of dish - to have the rail outside the FoV zd has to be 49 (see attached picture taken at zd50, az-113.5)
     207     52 -100    intermediate position of dish on right side
     208     58 -92     intermediate position of dish on right side
     209     65 -89     dish right side
     210     71 -85     tower right side
     211     75 -84     outmost knot right side
     212     */
     213
     214    bool HasShadowFromLST(const ZdAzPosn &hrz) const
     215    {
     216        return (hrz.zd > 0.0379373*pow(hrz.az+109.23, 2) + 49.1224)
     217            || (hrz.az>-79 && hrz.az<=-59 && hrz.zd>=70);
     218    }
     219
    175220    bool valid(const SolarObjects &so) const
    176221    {
    177         const HrzPosn hrz     = GetHrzFromEqu(equ, so.fJD);
    178         const double  current = FACT::PredictI(so, equ);
     222        const ZdAzPosn hrz      = GetHrzFromEqu(equ, so.fJD+timeshift);
     223        const ZdAzPosn true_hrz = timeshift==0 ? hrz : ZdAzPosn(GetHrzFromEqu(equ, so.fJD));
     224        const double   current  = FACT::PredictI(so, equ);
    179225
    180226        if (current>max_current)
    181227            return false;
    182228
    183         if (hrz.alt<=0 || 90-hrz.alt>max_zd)
     229        if (true_hrz.zd>max_zd)
    184230            return false;
    185231
    186         if (maxzd.valid && 90-hrz.alt>maxzd.val)
     232        if (maxzd.valid && true_hrz.zd>maxzd.val)
    187233            return false;
    188234
     235        if (lst_shadow && HasShadowFromLST(true_hrz))
     236            return false;
     237
    189238        if (maxcurrent.valid && current>maxcurrent.val)
     239            return false;
     240
     241        if (Nova::GetAngularSeparation(so.fMoonEqu, equ)<min_moon)
    190242            return false;
    191243
     
    205257    double getThreshold(const SolarObjects &so) const
    206258    {
    207         const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD);
    208         const double  current = FACT::PredictI(so, equ);
    209 
    210         const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
    211 
    212         return penalty*ratio*pow(current/6.2, 0.394);
     259        const ZdAzPosn hrz     = GetHrzFromEqu(equ, so.fJD+timeshift);
     260        const double   current = FACT::PredictI(so, equ);
     261
     262        const double ratio = pow(cos(hrz.zd*M_PI/180), -2.664);
     263        const double ratio2 = pow(current/6.2, 0.394);
     264
     265        return penalty*visratio*ratio*ratio2;
    213266    }
    214267
    215268    bool calcThreshold(const SolarObjects &so)
    216269    {
    217         const HrzPosn hrz     = GetHrzFromEqu(equ, so.fJD);
    218         const double  current = FACT::PredictI(so, equ);
     270        const ZdAzPosn hrz      = GetHrzFromEqu(equ, so.fJD+timeshift);
     271        const ZdAzPosn true_hrz = timeshift==0 ? hrz : ZdAzPosn(GetHrzFromEqu(equ, so.fJD));
     272        const double   current  = FACT::PredictI(so, equ);
    219273
    220274        if (current>max_current)
    221275            return false;
    222276
    223         if (hrz.alt<=0 || 90-hrz.alt>max_zd)
     277        if (true_hrz.zd>max_zd)
    224278            return false;
    225279
    226         if (maxzd.valid && 90-hrz.alt>maxzd.val)
     280        if (maxzd.valid && true_hrz.zd>maxzd.val)
     281            return false;
     282
     283        if (lst_shadow && HasShadowFromLST(true_hrz))
    227284            return false;
    228285
     
    230287            return false;
    231288
    232         const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
    233         threshold = penalty*ratio*pow(current/6.2, 0.394);
     289        if (Nova::GetAngularSeparation(so.fMoonEqu, equ)<min_moon)
     290            return false;
     291
     292        const double ratio = pow(cos(hrz.zd*M_PI/180), -2.664);
     293        const double ratio2 = pow(current/6.2, 0.394);
     294
     295        threshold = penalty*visratio*ratio*ratio2;
    234296
    235297        return true;
     
    239301double Source::max_zd;
    240302double Source::max_current;
     303double Source::min_moon;
     304double Source::vis_exponent;
     305double Source::ts_exponent;
     306bool Source::vis_ratio;
     307bool Source::lst_shadow;
    241308
    242309bool SortByThreshold(const Source &i, const Source &j) { return i.threshold<j.threshold; }
     
    367434bool RescheduleIntermediateSources(vector<Source> &obs)
    368435{
     436    bool changed = false;
     437
    369438    for (size_t i=1; i<obs.size()-1; i++)
    370439    {
     
    413482        cout << "Intermediate source [" << obs[i].name << "] removed" << endl;
    414483
    415         const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;
     484        // const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;
    416485
    417486        obs[i-1].end   = intersection;
     
    421490        i--;
    422491
     492        changed = true;
     493
    423494        if (obs.size()>1 && obs[i].name==obs[i+1].name)
    424495        {
     
    433504        }
    434505
    435         if (underflow)
    436             cout << "WARNING - Neighbor source < 40min as well." << endl;
    437     }
    438     return false;
    439 }
    440 
    441 void RemoveMiniSources(vector<Source> &obs)
    442 {
     506        // if (underflow)
     507        //     cout << "WARNING - Neighbor source < 40min as well." << endl;
     508    }
     509    return changed;
     510}
     511
     512// Remove mini source is able to violate limits
     513bool RemoveMiniSources(vector<Source> &obs)
     514{
     515    bool changed = false;
     516
    443517    for (size_t i=1; i<obs.size()-1; i++)
    444518    {
     
    446520            continue;
    447521
    448         if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
     522        // if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
     523        //     continue;
     524
     525        cout << "Mini source [" << obs[i].name << "] detected < 5min" << endl;
     526
     527        if (obs[i-1].name==obs[i+1].name)
     528        {
     529            cout << "Combined surrounding identical sources [" << obs[i-1].name << "]" << endl;
     530
     531            obs[i-1].end = obs[i+1].end;
     532            obs.erase(obs.begin()+i, obs.begin()+i+1);
     533            i -= 2;
     534
     535            changed = true;
     536
    449537            continue;
    450 
    451         cout << "Mini source [" << obs[i].name << "] detected < 5min" << endl;
    452 
     538        }
     539
     540        /*
    453541        if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
    454542        {
     
    460548            i -= 2;
    461549
     550            changed = true;
     551
    462552            cout << "Combined two surrounding sleep into one" << endl;
    463553
    464554            continue;
    465         }
     555        }*/
     556
     557        // Note that this (intentionally) violates our safety limits
    466558
    467559        if (obs[i-1].name=="SLEEP")
    468560        {
     561            cout << "Preceding SLEEP detected... extended next source [" << obs[i+1].name << "]" << endl;
     562
    469563            obs[i-1].end = obs[i+1].begin;
    470564            obs.erase(obs.begin()+i);
    471565            i--;
    472566
    473             cout << "Extended previous sleep" << endl;
     567            changed = true;
    474568
    475569            continue;
     
    478572        if (obs[i+1].name=="SLEEP")
    479573        {
    480             obs[i+1].begin = obs[i-1].end;
     574            cout << "Following SLEEP detected... extended previous source [" << obs[i-1].name << "]" << endl;
     575
     576            obs[i-1].end = obs[i+1].begin;
    481577            obs.erase(obs.begin()+i);
    482 
    483             cout << "Extended following sleep" << endl;
    484 
    485578            i--;
     579
     580            changed = true;
     581
    486582            continue;
    487583        }
    488     }
     584
     585        // Algorithm based on RescheduleIntermediateSources
     586        // Note that this (intentionally) violates our safety limits
     587
     588        const SolarObjects so1(obs[i-1].end);
     589        const SolarObjects so2(obs[i+1].begin);
     590
     591        const double th1 = obs[i-1].getThreshold(so1);
     592        const double th2 = obs[i+1].getThreshold(so2);
     593
     594        const double intersection = th1<th2 ? obs[i+1].begin : obs[i-1].end;
     595
     596        cout << "Intermediate mini source [" << obs[i].name << "] detected" << endl;
     597        cout << (th1<th2?"Previous":"Next") << " source [" << (th1<th2?obs[i-1].name:obs[i+1].name) << "] extended" << endl;
     598
     599        obs[i-1].end   = intersection;
     600        obs[i+1].begin = intersection;
     601        obs.erase(obs.begin()+i);
     602
     603        i--;
     604
     605        changed = true;
     606    }
     607
     608    return changed;
    489609}
    490610
     
    504624}
    505625
    506 void Print(const vector<Source> &obs, double startup_offset)
    507 {
    508     cout << Time(obs[0].begin-startup_offset).GetAsStr() << "  STARTUP\n";
     626void Print(char id, const vector<Source> &obs, double startup_offset)
     627{
     628    cout << id << " " << Time(obs[0].begin-startup_offset).GetAsStr() << "  STARTUP\n";
    509629    for (const auto& src: obs)
    510630    {
     
    514634            for (const auto& pre: src.preobs)
    515635            {
    516                 cout << tm << "  " << pre << "\n";
     636                cout << id << " " << tm << "  " << pre << "\n";
    517637                tm = "                   ";
    518638            }
    519639        }
    520640
    521         cout << tm << "  " << src.name << " [";
     641        cout << id << " " << tm << "  " << src.name << " [";
    522642        cout << src.duration()*24*60 << "'";
    523643        if (src.name!="SLEEP")
    524             cout << Tools::Form("; %.1f/%.1f", src.zd(src.begin), src.zd(src.end));
     644            cout << Tools::Form("; %.1f/%.1f; %.2f; %.2fh", src.zd(src.begin), src.zd(src.end), src.visratio, src.timeshift*24);
    525645        cout << "]";
    526646
     
    530650        cout << "\n";
    531651    }
    532     cout << Time(obs.back().end).GetAsStr() << "  SHUTDOWN" << endl;
     652    cout << id << " " << Time(obs.back().end).GetAsStr() << "  SHUTDOWN" << endl;
     653}
     654
     655void PrintRect(ofstream &svg, const uint16_t &x, const int16_t &y, const uint16_t &w, const uint16_t &h, const Source &src)
     656{
     657    static vector<string> colors =
     658    {
     659        "gray",     // rgb(0,0,0)
     660        "red",
     661        "olive",
     662        "blue",
     663        "maroon",
     664        "green",
     665        "teal",
     666        "navy",
     667        "purple",
     668        "silver",
     669        "lime",
     670        "fuchsia",
     671        "yellow",
     672    };
     673
     674    static map<uint16_t, uint16_t> col_keys;
     675
     676    const auto ic = col_keys.find(src.key);
     677    if (ic==col_keys.end())
     678    {
     679        if (col_keys.size()>=colors.size())
     680            cerr << "WARNING: More different sources (" << col_keys.size() << ") than colors (" << colors.size() << ")." << endl;
     681
     682        const auto n = col_keys.size();
     683
     684        col_keys[src.key] = n%colors.size();
     685
     686        svg << "<g>\n";
     687        svg << " <rect x=\"" << 0 << "\" y=\"" << n*18    << "\" width=\"140\" height=\"18\" fill=\"" << colors[col_keys[src.key]] << "\"/>\n";
     688        svg << " <text x=\"" << 0 << "\" y=\"" << n*18+14 << "\" font-family=\"Verdana\" font-size=\"16\" fill=\"white\">" << src.name << "</text>\n";
     689        svg << "</g>\n";
     690    }
     691
     692    const string color = colors[col_keys[src.key]];
     693
     694    svg << "<g>\n";
     695    svg << " <defs>\n";
     696    svg << "  <linearGradient x1=\"0\" x2=\"0\" y1=\"0\" y2=\"1\" id=\"grad\">\n";
     697
     698    for (int i=0; i<11; i++)
     699    {
     700        const double opa = 1-src.zd(src.begin+(src.end-src.begin)*0.10*i)/(src.maxzd.valid?src.maxzd.val:src.max_zd);
     701        svg << "   <stop style=\"stop-color:" << color << "; stop-opacity:" << opa << "\" offset=\"" << i*10 << "%\"/>\n";
     702    }
     703    svg << "  </linearGradient>\n";
     704    svg << " </defs>\n";
     705
     706    svg << " <rect";
     707    svg << " x=\"" << x*w+150 << "\"";
     708    svg << " y=\"" << y << "\"";
     709    svg << " width=\""  << w << "\"";
     710    svg << " height=\"" << h << "\"";
     711    svg << " style=\"fill:url(#grad);\"";
     712    svg << "/>\n";
     713    svg << "</g>\n";
     714}
     715
     716
     717void PrintSVG(ofstream &svg, const uint16_t &scale, const Time &T0, const vector<Source> &obs)
     718{
     719    const double jd = floor(T0.JD());
     720
     721    for (const auto& src: obs)
     722    {
     723        if (src.preobs.size()>0 || src.key==0 || src.name=="SLEEP")
     724            continue;
     725
     726        const double diff = src.begin-jd;
     727
     728        const uint32_t b = 24*60*diff;
     729        const uint32_t d = 24*60*(src.end-src.begin)+0.5;
     730
     731        PrintRect(svg, diff, b%(24*60)-6*60, scale, d, src);
     732    }
    533733}
    534734
     
    594794}
    595795
    596 int main(int argc, const char* argv[])
    597 {
    598 //    gROOT->SetBatch();
    599 
    600     Configuration conf(argv[0]);
    601     conf.SetPrintUsage(PrintUsage);
    602     SetupConfiguration(conf);
    603 
    604     if (!conf.DoParse(argc, argv, PrintHelp))
    605         return 127;
    606 
    607     // ------------------ Eval config ---------------------
    608 
    609     const int enter = conf.Has("enter-schedule-into-database") ? (conf.Get<bool>("enter-schedule-into-database") ? 1 : -1) : 0;
    610     if (enter && !conf.Has("schedule-database"))
    611         throw runtime_error("enter-schedule-into-database required schedule-database.");
    612 
    613     Time time;
    614     if (conf.Has("date"))
    615         time.SetFromStr(conf.Get<string>("date")+" 12:00:00");
    616 
    617     if (enter && floor(time.JD())<ceil(Time().JD()))
    618         throw runtime_error("Only future schedules can be entered into the database.");
    619 
    620     Source::max_current = conf.Get<double>("max-current");
    621     Source::max_zd      = conf.Get<double>("max-zd");
    622 
    623     const double startup_offset = conf.Get<double>("startup.offset")/60/24;
    624 
    625     const double angle_sun_set  = conf.Get<double>("data-taking.start");
    626     const double angle_sun_rise = conf.Get<double>("data-taking.end");
    627 
    628     if (startup_offset<0 || startup_offset>120)
    629         throw runtime_error("Only values [0;120] are allowed for startup.offset");
    630 
    631     if (angle_sun_set>-6)
    632         throw runtime_error("datataking.start not allowed before sun at -6deg");
    633 
    634     if (angle_sun_rise>-6)
    635         throw runtime_error("datataking.end not allowed after sun at -6deg");
    636 
    637     // -12: nautical
    638     // Sun set with the same date than th provided date
    639     // Sun rise on the following day
    640     const RstTime sun_set  = GetSolarRst(floor(time.JD())-0.5, angle_sun_set);
    641     const RstTime sun_rise = GetSolarRst(floor(time.JD())+0.5, angle_sun_rise);
    642 
    643     const double sunset  = ceil(sun_set.set*24*60)   /24/60 + 1e-9;
    644     const double sunrise = floor(sun_rise.rise*24*60)/24/60 + 1e-9;
    645 
    646     cout << "\n";
    647 
    648     cout << "Date: " << Time(floor(sunset)).GetAsStr() << "\n";
    649     cout << "Set:  " << Time(sunset).GetAsStr()  << "  [" << Time(sun_set.set)   << "]\n";
    650     cout << "Rise: " << Time(sunrise).GetAsStr() << "  [" << Time(sun_rise.rise) << "]\n";
    651 
    652     cout << "\n";
    653 
    654     cout << "Global maximum current:  " << Source::max_current << " uA/pix\n";
    655     cout << "Global zenith distance:  " << Source::max_zd << " deg\n";
    656 
    657     cout << "\n";
    658 
    659     // ------------- Get Sources from databasse ---------------------
    660 
    661     const vector<string> ourcenames = conf.Vec<string>("source");
     796vector<Source> GetSources(Configuration &conf)
     797{
    662798    const vector<string> sourcenames = conf.Vec<string>("source");
    663799    cout << "Nsources = " << sourcenames.size() << "\n";
     
    676812    for (const auto &row: res)
    677813    {
     814        //cout << endl;
    678815        const string name = string(row[0]);
    679816
     
    687824        src.penalty = conf.Has("setup."+name+".penalty") ?
    688825            conf.Get<double>("setup."+name+".penalty") : 1;
     826        src.visexponent = conf.Has("setup."+name+".exponent-visibility-ratio") ?
     827            conf.Get<double>("setup."+name+".exponent-visibility-ratio") : 1;
     828        src.tsexponent = conf.Has("setup."+name+".exponent-time-shift") ?
     829            conf.Get<double>("setup."+name+".exponent-time-shift") : 1;
     830        src.time_shift = conf.Has("setup."+name+".time-shift") ?
     831            conf.Get<double>("setup."+name+".time-shift")/24 : 0;
    689832
    690833        src.preobs = conf.Vec<string>("preobs."+name);
    691834
    692 
    693         cout << "[" << name << "]";
     835        cout << "[" << src.key << ":" << name << "]";
    694836
    695837        if (src.maxzd.valid)
     
    697839        if (src.penalty!=1)
    698840            cout << " Penalty=" << src.penalty;
     841        if (src.visexponent!=1)
     842            cout << " Exponent [VR]=" << src.visexponent;
     843        if (src.tsexponent!=1)
     844            cout << " Exponent [TS]=" << src.tsexponent;
     845        if (src.time_shift!=0)
     846            cout << " Time-shift=" << src.time_shift*24 << "h";
    699847
    700848        cout << " " << boost::algorithm::join(src.preobs, "+") << endl;
    701849
     850        for (int k=0; k<90; k++)
     851            src.zddistr[k]=0;
    702852        /*
    703853         RstTime t1 = GetObjectRst(floor(sunset)-1, src.equ);
     
    713863    cout << endl;
    714864
    715     // -------------------------------------------------------------------------
    716 
    717     vector<Source> obs;
    718 
    719     const uint32_t n = nearbyint((sunrise-sunset)*24*60);
    720     for (uint32_t i=0; i<n; i++)
    721     {
    722         const double jd = sunset + i/24./60.;
    723 
    724         const SolarObjects so(jd);
    725 
    726         vector<Source> vis;
     865    return sources;
     866}
     867
     868int main(int argc, const char* argv[])
     869{
     870//    gROOT->SetBatch();
     871
     872    Configuration conf(argv[0]);
     873    conf.SetPrintUsage(PrintUsage);
     874    SetupConfiguration(conf);
     875
     876    if (!conf.DoParse(argc, argv, PrintHelp))
     877        return 127;
     878
     879    // ------------------ Eval config ---------------------
     880
     881    const int enter = conf.Has("enter-schedule-into-database") ? (conf.Get<bool>("enter-schedule-into-database") ? 1 : -1) : 0;
     882    if (enter && !conf.Has("schedule-database"))
     883        throw runtime_error("enter-schedule-into-database requires schedule-database.");
     884
     885    Time time;
     886    if (conf.Has("date"))
     887        time.SetFromStr(conf.Get<string>("date")+" 12:00:00");
     888    else
     889        time=floor(Time().JD());
     890
     891    /*
     892    cout << "-------> " << setprecision(12) << time.JD() << endl;
     893    cout << "         " << time.GetAsStr()
     894        << " " << ceil(Time().JD())
     895        << " " << floor(Time().JD())
     896        << endl;
     897    */
     898    const uint16_t loop = conf.Get<uint16_t>("loop");
     899
     900    if (enter && floor(time.JD())<ceil(Time().JD()))
     901        throw runtime_error("Only future schedules can be entered into the database.");
     902
     903    Source::max_current = conf.Get<double>("max-current");
     904    Source::max_zd      = conf.Get<double>("max-zd");
     905    Source::min_moon    = conf.Get<double>("min-moon");
     906    Source::vis_exponent= conf.Get<double>("exponent-visibility-ratio");
     907    Source::ts_exponent = conf.Get<double>("exponent-time-shift");
     908    Source::vis_ratio   = conf.Get<bool>("use-visibility-ratio");
     909    Source::lst_shadow  = conf.Get<bool>("use-lst-shadow");
     910
     911    const double startup_offset = conf.Get<double>("startup.offset")/60/24;
     912
     913    const double angle_sun_set  = conf.Get<double>("data-taking.start");
     914    const double angle_sun_rise = conf.Get<double>("data-taking.end");
     915
     916    if (startup_offset<0 || startup_offset>120)
     917        throw runtime_error("Only values [0;120] are allowed for startup.offset");
     918
     919    if (angle_sun_set>-6)
     920        throw runtime_error("datataking.start not allowed before sun at -6deg");
     921
     922    if (angle_sun_rise>-6)
     923        throw runtime_error("datataking.end not allowed after sun at -6deg");
     924
     925    cout << "\n";
     926
     927    cout << "Global maximum current allowed:  " << Source::max_current << " uA/pix\n";
     928    cout << "Global maximum zenith distance:  " << Source::max_zd      << " deg\n";
     929    cout << "Min. angular separation to moon: " << Source::min_moon    << " deg\n";
     930
     931    cout << "\n";
     932
     933    // ------------- SVG ---------------------------------------------
     934
     935    const uint16_t svgscale = loop>0 ? ::min((2*1080)/loop, 5) : 0;
     936
     937    ofstream svg;
     938    if (conf.Has("svg") && loop>0)
     939    {
     940        svg.open(conf.Get<string>("svg"));
     941        if (!svg)
     942        {
     943            cerr << "ERROR: Writing to '" << conf.Get<string>("svg") << "' failed: " << strerror(errno) << endl;
     944            return 1;
     945        }
     946
     947        svg << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
     948        svg << "<svg width=\""  << loop*svgscale+150 << "\" height=\"840\" xmlns=\"http://www.w3.org/2000/svg\">\n";
     949        svg << "<rect width=\"" << loop*svgscale+150 << "\" height=\"840\" fill=\"black\"/>\n";
     950
     951        //svg << "<svg width=\"" << loop*svgscale+150 << "\" height=\"840\">\n";
     952        //svg << "<rect x=\"0\" y=\"0\" width=\"" << loop*svgscale+150 << "\" height=\"840\" fill=\"black\"/>\n";
     953    }
     954
     955    // ------------- Get Sources from databasse ---------------------
     956
     957    vector<Source> sources = GetSources(conf);
     958
     959    //cout << "VR: " << Source::vis_ratio << endl;
     960
     961    const uint16_t N = ::max<uint16_t>(loop, 1);
     962    for (int iday=0; iday<N; iday++)
     963    {
     964        // ----------- Define time slots for the day ---------------------
     965
     966        // -12: nautical
     967        // Sun set with the same date than th provided date
     968        // Sun rise on the following day
     969        const RstTime sun_set  = GetSolarRst(floor(time.JD()+iday)-0.5, angle_sun_set);
     970        const RstTime sun_rise = GetSolarRst(floor(time.JD()+iday)+0.5, angle_sun_rise);
     971
     972        const double sunset  = ceil(sun_set.set*24*60)   /24/60 + 1e-9;
     973        const double sunrise = floor(sun_rise.rise*24*60)/24/60 + 1e-9;
     974
     975        cout << "Date: " << Time(floor(sunset)).GetAsStr() << "\n";
     976        cout << "Set:  " << Time(sunset).GetAsStr()  << "  [" << Time(sun_set.set)   << "]\n";
     977        cout << "Rise: " << Time(sunrise).GetAsStr() << "  [" << Time(sun_rise.rise) << "]\n";
     978        cout << "\n";
     979
     980        const double vis_sun = sun_rise.rise - sun_set.set;
     981
    727982        for (auto& src: sources)
    728983        {
    729             if (src.calcThreshold(so))
    730                 vis.emplace_back(src);
     984            const double maxzd = min(Source::max_zd, src.maxzd.valid?src.maxzd.val:90);
     985
     986            int8_t sign = 0;
     987
     988            RstTime rst;
     989            const int rc = GetObjectRst(rst, src.equ, sun_set.set, 90-maxzd);
     990            switch (rc)
     991            {
     992            case -1:
     993                // circumpolar below horizon - not visible
     994                src.visratio = 0;
     995                break;
     996
     997            case 0:
     998                {
     999                    // Although, according to the documentation, everything should be alinged
     1000                    // between JD and JD+1, it is not... therefore, we align ourselves
     1001                    if (rst.rise-sun_set.set>1)
     1002                        rst.rise -= 1;
     1003                    if (rst.set-sun_set.set>1)
     1004                        rst.set -= 1;
     1005                    if (rst.rise-sun_set.set<0)
     1006                        rst.rise += 1;
     1007                    if (rst.set-sun_set.set<0)
     1008                        rst.set += 1;
     1009
     1010                    // Now: rst.rise/rst.set is aligned between sunset and sunset+1:
     1011                    //
     1012                    //              JD                                 JD+1
     1013                    // 1) SUN.SET |                        | sun.rise  RST.RISE rst.set  ... SUN.SET |
     1014                    // 2) SUN.SET |++++++++++++++++++++++++| sun.rise  rst.set  RST.RISE ... SUN.SET |
     1015                    // 3) SUN.SET |       RST.RISE+++++++++| sun.rise  rst.set           ... SUN.SET |
     1016                    // 4) SUN.SET |++++++++rst.set         | sun.rise  rst.rise          ... SUN.SET |
     1017                    // 5) SUN.SET |  RST.RISE+++rst.set    | sun.rise                    ... SUN.SET |
     1018                    // 6) SUN.SET |+++rst.set   RST.RISE+++| sun.rise                    ... SUN.SET |
     1019
     1020                    if (rst.rise<sun_rise.rise && sun_rise.rise<rst.set)
     1021                        sign = 1;
     1022                    if (rst.set<sun_rise.rise && sun_rise.rise<rst.rise)
     1023                        sign = -1;
     1024
     1025                    // Identify case 6
     1026                    const bool case6 = rst.set<rst.rise && rst.rise<sun_rise.rise;
     1027
     1028                    // make sure that the visibility calculation of the source
     1029                    // yields positive valiues and the overlap calculation works
     1030                    if (rst.set<rst.rise)
     1031                        rst.rise -= 1;
     1032
     1033                    // Now: rst.rise is always before rst.set
     1034                    //
     1035                    //                    JD                               JD+1
     1036                    // 1) ...          SUN.SET |                        | sun.rise  RST.RISE rst.set ...
     1037                    // 2) ... RST.RISE SUN.SET |++++++++++++++++++++++++| sun.rise  rst.set          ...
     1038                    // 3) ...          SUN.SET |       RST.RISE+++++++++| sun.rise  rst.set          ...
     1039                    // 4) ... RST.RISE SUN.SET |++++++++rst.set         | sun.rise  rst.rise         ...
     1040                    // 5) ...          SUN.SET |  RST.RISE+++rst.set    | sun.rise                   ...
     1041                    // 6) ... RST.RISE SUN.SET |+++rst.set  [RST.RISE]++| sun.rise                   ...
     1042
     1043                    // Total visibility of the source itself (this is now always positive)
     1044                    const double vis_rst = rst.set-rst.rise;
     1045
     1046                    // Calculate overlap region
     1047                    // Visibility of the source between sunset and sunrise
     1048                    const double vis_night = case6 ? vis_sun - (1-vis_rst) : min(rst.set, sun_rise.rise) - max(rst.rise, sun_set.set);
     1049
     1050                    //vis_sun-(1-vis_rst) = vis_sun-(1-(set-rise)) = vis_sun-(1-(set-(rise-1))) = vis_sun-(1-set+(rise-1)) = vis_sun-(rise-set)
     1051
     1052                    // General visibility of the source (but maximum one night)
     1053                    const double vis_total = min(vis_rst, vis_sun);
     1054
     1055                    // Special treatment of case 1 (vis_night<0 means that there is no overlap)
     1056                    src.visratio = vis_night<0 ? 0 : vis_night / vis_total;
     1057
     1058                    // cout << src.name << "   " << Time(rst.rise).GetAsStr() << "   " << Time(rst.set).GetAsStr() << "   " << vis_night << " " << vis_total << " " << vis_rst << " " << vis_sun << " " << case6 << " " << test << endl;
     1059                    break;
     1060                }
     1061
     1062            case 1:
     1063                // circumpolar above horizon -> always visible
     1064                src.visratio = 1;
     1065                break;
     1066            }
     1067
     1068            src.timeshift = src.time_shift * sign*(1-pow(src.visratio, src.tsexponent*Source::ts_exponent));
     1069
     1070            src.visratio = Source::vis_ratio ? pow(src.visratio, src.visexponent*Source::vis_exponent) : 1;
    7311071        }
    7321072
    733         // In case no source was found, add a sleep source
    734         Source src("SLEEP");
    735         vis.emplace_back(src);
    736 
    737         // Source has higher priority if minimum observation time not yet fullfilled
    738         sort(vis.begin(), vis.end(), SortByThreshold);
    739 
    740         if (obs.size()>0 && obs.back().name==vis[0].name)
     1073        // -------------------------------------------------------------------------
     1074
     1075        vector<Source> obs;
     1076
     1077        const uint32_t n = nearbyint((sunrise-sunset)*24*60);
     1078        for (uint32_t i=0; i<n; i++)
     1079        {
     1080            const double jd = sunset + i/24./60.;
     1081
     1082            const SolarObjects so(jd);
     1083
     1084            vector<Source> vis;
     1085            for (auto& src: sources)
     1086            {
     1087                if (src.calcThreshold(so))
     1088                    vis.emplace_back(src);
     1089            }
     1090
     1091            // In case no source was found, add a sleep source
     1092            Source src("SLEEP");
     1093            vis.emplace_back(src);
     1094
     1095            // Source has higher priority if minimum observation time not yet fullfilled
     1096            sort(vis.begin(), vis.end(), SortByThreshold);
     1097
     1098            if (obs.size()>0 && obs.back().name==vis[0].name)
     1099                continue;
     1100
     1101            vis[0].begin = jd;
     1102            obs.emplace_back(vis[0]);
     1103        }
     1104
     1105        if (obs.size()==0)
     1106        {
     1107            cout << "No source found." << endl;
     1108            return 1;
     1109        }
     1110
     1111        // -------------------------------------------------------------------------
     1112
     1113        for (auto it=obs.begin(); it<obs.end()-1; it++)
     1114            it[0].end = it[1].begin;
     1115        obs.back().end = sunrise;
     1116
     1117        // -------------------------------------------------------------------------
     1118
     1119        Print('-', obs, startup_offset);
     1120        cout << endl;
     1121
     1122        // -------------------------------------------------------------------------
     1123
     1124        while (RescheduleFirstSources(obs));
     1125        while (RescheduleLastSources(obs));
     1126        while (RescheduleIntermediateSources(obs));
     1127        while (RemoveMiniSources(obs));
     1128
     1129        CheckStartupAndShutdown(obs);
     1130
     1131        // ---------------------------------------------------------------------
     1132
     1133        cout << endl;
     1134        Print('*', obs, startup_offset);
     1135        cout << endl;
     1136
     1137        // ---------------------------------------------------------------------
     1138
     1139        //fill zd-distribution
     1140        for (auto& src2: sources)
     1141        {
     1142            //cout << "source: " << src2.name << " " << src2.key << endl;
     1143            for (const auto& src: obs)
     1144            {
     1145                //cout << "source: " << src.name << " " << src.key << endl;
     1146                const uint32_t n2 = nearbyint((src.end-src.begin)*24*60);
     1147                for (uint32_t i=0; i<n2; i++)
     1148                {
     1149                    const double jd = src.begin + i/24./60.;
     1150                    const int zd = (int)src.zd(jd);
     1151                    //cout << " zd:" << floor(zd) << " " << zd << endl;
     1152                    //cout << " zd: " << zd << " " << src2.zddistr[zd] << flush;
     1153                    if (src.key==src2.key)
     1154                    {
     1155                        //cout << " -> " << flush;
     1156                        src2.zddistr[zd]++;
     1157                    }
     1158                    //cout << " " << src2.zddistr[zd] << endl;
     1159                }
     1160            }
     1161        }
     1162
     1163        // ---------------------------------------------------------------------
     1164
     1165        if (svg.is_open())
     1166            PrintSVG(svg, svgscale, time, obs);
     1167
     1168        // ---------------------------------------------------------------------
     1169
     1170        if (loop>0)
    7411171            continue;
    7421172
    743         vis[0].begin = jd;
    744         obs.emplace_back(vis[0]);
    745     }
    746 
    747     if (obs.size()==0)
    748     {
    749         cout << "No source found." << endl;
     1173        if (!enter)
     1174            return 0;
     1175
     1176        const string scheduledb = conf.Get<string>("schedule-database");
     1177
     1178        Database db(scheduledb);
     1179
     1180        if (enter>0)
     1181            db.query("LOCK TABLES Schedule WRITE");
     1182
     1183        const int rc = FillSql(db, enter, obs, startup_offset);
     1184
     1185        if (enter>0)
     1186            db.query("UNLOCK TABLES");
     1187
     1188        return rc;
     1189    }
     1190
     1191    if (svg.is_open())
     1192        svg << "</svg>" << endl;
     1193
     1194    // ---------------------------------------------------------------------
     1195    if (conf.Has("print-hist") && conf.Get<bool>("print-hist"))
     1196    {
     1197        for (const auto& src: sources)
     1198        {
     1199            cout << "ZD-DISTRIBUTION [" << src.name << "]: " << flush;
     1200            //print zd-distr
     1201            for (int k=0; k<90; k++)
     1202                cout << src.zddistr[k] << " " << flush;
     1203            cout << endl;
     1204            //plot zd-distr
     1205            for (int k=0; k<90; k++)
     1206            {
     1207                if (src.zddistr[k]==0)
     1208                    continue;
     1209                printf("%02d", k);
     1210                for (int k2=0; k2<(int)(src.zddistr[k]/loop); k2++)
     1211                    cout << "*" << flush;
     1212                cout << endl;
     1213            }
     1214        }
     1215    }
     1216
     1217    if (!conf.Has("hist"))
     1218        return 0;
     1219
     1220    ofstream hist(conf.Get<string>("hist").c_str());
     1221    if (!hist)
     1222    {
     1223        cerr << "ERROR: Writing to '" << conf.Get<string>("hist") << "' failed: " << strerror(errno) << endl;
    7501224        return 1;
    7511225    }
    7521226
    753     // -------------------------------------------------------------------------
    754 
    755     for (auto it=obs.begin(); it<obs.end()-1; it++)
    756         it[0].end = it[1].begin;
    757     obs.back().end = sunrise;
    758 
    759     // -------------------------------------------------------------------------
    760 
    761     Print(obs, startup_offset);
    762     cout << endl;
    763 
    764     // -------------------------------------------------------------------------
    765 
    766     while (RescheduleFirstSources(obs));
    767     while (RescheduleLastSources(obs));
    768     while (RescheduleIntermediateSources(obs));
    769 
    770     RemoveMiniSources(obs);
    771     CheckStartupAndShutdown(obs);
    772 
    773     // ---------------------------------------------------------------------
    774 
    775     cout << endl;
    776     Print(obs, startup_offset);
    777     cout << endl;
    778 
    779     // ---------------------------------------------------------------------
    780 
    781     if (!enter)
    782         return 0;
    783 
    784     const string scheduledb = conf.Get<string>("schedule-database");
    785 
    786     Database db(scheduledb);
    787 
    788     if (enter>0)
    789         db.query("LOCK TABLES Schedule WRITE");
    790 
    791     const int rc = FillSql(db, enter, obs, startup_offset);
    792 
    793     if (enter>0)
    794         db.query("UNLOCK TABLES");
    795 
    796     // ---------------------------------------------------------------------
    797 
    798     return rc;
    799 }
     1227    hist << "# " << flush;
     1228    for (const auto& src: sources)
     1229        hist << src.name << " " << flush;
     1230    hist << endl;
     1231
     1232    hist << "# " << flush;
     1233    for (const auto& src: sources)
     1234        hist << src.key << " " << flush;
     1235    hist << endl;
     1236
     1237    for (int k=0; k<90; k++)
     1238    {
     1239        hist << k << " " << flush;
     1240        for (const auto& src: sources)
     1241            hist << src.zddistr[k] << " " << flush;
     1242        hist << endl;
     1243    }
     1244
     1245    return 0;
     1246}
Note: See TracChangeset for help on using the changeset viewer.