Changeset 17651 for trunk/FACT++/src/makedata.cc

Timestamp:

04/04/14 11:02:41 (10 years ago)

Author:

tbretz

Message:

Updated with the current prediction from the calibration paper; make plots only for VHE sources.

File:

• trunk/FACT++/src/makedata.cc (modified) (2 diffs)

trunk/FACT++/src/makedata.cc

@@ -99 +99 @@
     // ------------------ Precalc moon coord ---------------------
 
-    vector<pair<EquPosn, double>> fMoonCoords;
+    vector<tuple<EquPosn, double, double>> fMoonCoords;
+    vector<EquPosn> fSunCoords;
 
     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);
+        fMoonCoords.emplace_back(moon, disk, edist);
+        fSunCoords.emplace_back(sun);
     }
 
@@ -136 +140 @@
         const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
 
-        // Get moon properties and
-        const EquPosn moon = fMoonCoords[i].first;
-        const double  disk = fMoonCoords[i].second;
-        const HrzPosn hrzm = GetHrzFromEqu(moon, 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);
 
-        // Angular distance between source and moon
+        // Distance between source and moon
         const double angle = GetAngularSeparation(moon, pos);
 
-        // Distance between earth and moon relative to major semi-axis
-        const double dist  = GetLunarEarthDist(jd+h)/384400;
+        // 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);
 
-        // Current prediction
-        const double calt = sin(hrzm.alt*M_PI/180);
-        const double cang = 1-sin(angle*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 lc = cang==0 ? -1 : calt*pow(disk, 2.2)*pow(dist, -2);
-        //const double cur = lc>0 ? 4+103*lc : 4;
-
-        const double lc  = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2);
-        const double cur = lc>0 ? 7.2 + 69*lc : 7.2;
+        const double cur = 6.2 + 95.7*k0*k1*k2*k3 + k4;
 
         // Relative  energy threshold prediction