Changeset 17957

Timestamp:

08/13/14 12:50:36 (10 years ago)

Author:

tbretz

Message:

Make use of new central current prediction.

File:

• trunk/FACT++/src/makeplots.cc (modified) (5 diffs)

trunk/FACT++/src/makeplots.cc

@@ -1 @@
-#include "externals/nova.h"
+#include "externals/Prediction.h"
 
 #include "Database.h"
@@ -121 +121 @@
     cout << "Rise: " << Time(sun_rise.rise) << endl;
 
-    const double jd0 = fmod(sun_set.set,   1);   // ~0.3
-    const double jd1 = fmod(sun_rise.rise, 1);   // ~0.8
+    const double sunset  = sun_set.set;
+    const double sunrise = sun_rise.rise;
 
     const string fDatabase = conf.Get<string>("source-database");
 
-    // ------------------ Precalc moon coord ---------------------
-
-    vector<tuple<EquPosn, double, double>> fMoonCoords;
-    vector<EquPosn> fSunCoords;
+    // ------------------ Precalc coord ---------------------
+
+    vector<SolarObjects> fCoordinates;
 
     for (double h=0; h<1; h+=1./(24*12))
-    {
-        const EquPosn sun  = GetSolarEquCoords(jd+h);
-        const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
-        const double  disk = GetLunarDisk(jd+h);
-        const double edist = GetLunarEarthDist(jd+h)/384400;
-
-        fMoonCoords.emplace_back(moon, disk, edist);
-        fSunCoords.emplace_back(sun);
-    }
+        if (jd+h>sunset && jd+h<sunrise)
+            fCoordinates.emplace_back(jd+h);
 
     // ------------- Get Sources from databasse ---------------------
@@ -152 +144 @@
     // 1970-01-01 00:00:00. This one will not work if your
     // local time zone is positive!
-    TH1S hframe("", "", 1, (mjd+jd0)*24*3600, (mjd+jd1)*24*3600);
+    TH1S hframe("", "", 1, Time(sunset).Mjd()*24*3600, Time(sunrise).Mjd()*24*3600);
     hframe.SetStats(kFALSE);
     hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
@@ -262 +254 @@
 
         // Loop over 24 hours
-        int i=0;
-        for (double h=0; h<1; h+=1./(24*12), i++)
+        for (auto it=fCoordinates.begin(); it!=fCoordinates.end(); it++)
         {
-            // check if it is betwene sun-rise and sun-set
-            if (h<jd0 || h>jd1)
-                continue;
-
             // get local position of source
-            const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
-
-            // Get sun and moon properties and
-            const EquPosn  sun   = fSunCoords[i];
-            const EquPosn  moon  = get<0>(fMoonCoords[i]);
-            const double   disk  = get<1>(fMoonCoords[i]);
-            const double   edist = get<2>(fMoonCoords[i]);
-            const HrzPosn  hrzm  = GetHrzFromEqu(moon, jd+h);
-            const ZdAzPosn hrzs  = GetHrzFromEqu(sun,  jd+h);
+            const HrzPosn hrz = GetHrzFromEqu(pos, it->fJD);
 
             if (v==res.begin())
-                cout << Time(jd+h) <<" " << 90-hrzm.alt <<  endl;
-
-            // Distance between source and moon
-            const double angle = GetAngularSeparation(moon, pos);
-
-            // Current prediction
-            const double sin_malt  = hrzm.alt<0 ? 0 : sin(hrzm.alt*M_PI/180);
-            const double cos_mdist = cos(angle*M_PI/180);
-            const double sin_szd   = sin(hrzs.zd*M_PI/180);
-
-            const double k0 = pow(disk,      2.63);
-            const double k1 = pow(sin_malt,  0.60);
-            const double k2 = pow(edist,    -2.00);
-            const double k3 = exp(0.67*cos_mdist*cos_mdist*cos_mdist*cos_mdist);
-            const double k4 = exp(-97.8+105.8*sin_szd*sin_szd);
-
-            const double cur = 6.2 + 95.7*k0*k1*k2*k3 + k4;
+                cout << Time(it->fJD) <<" " << 90-it->fMoonHrz.alt <<  endl;
+
+            const double cur = FACT::PredictI(*it, pos);
 
             // 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.;
             const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
 
             // Add points to curve
-            const double axis = (mjd+h)*24*3600;
+            const double axis = Time(it->fJD).Mjd()*24*3600;
 
             // If there is a gap of more than one bin, start a new curve
@@ -325 +288 @@
 
             if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
+            {
+                const double angle = GetAngularSeparation(it->fMoonEqu, pos);
                 g4.SetPoint(g4.GetN(), axis, angle);
+            }
 
             if (cnt==0)
-                gm.SetPoint(gm.GetN(), axis, hrzm.alt);
+                gm.SetPoint(gm.GetN(), axis, it->fMoonHrz.alt);
         }