Changeset 19852 for trunk/Mars/msim/MAtmosphere.cc

Timestamp:

11/01/19 15:11:07 (5 years ago)

Author:

tbretz

Message:

Changed the default atmosphere without user initailization to Keilhauer US Std. Make sure that the overburden below the obs level is 0. Split the Init functions such that the observation level can be changed indendently from the atmospshere. Copy the height levels from MRunHeader.

File:

• trunk/Mars/msim/MAtmosphere.cc (modified) (8 diffs)

trunk/Mars/msim/MAtmosphere.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: CheObs Software Development, 2000-2009
+!   Copyright: CheObs Software Development, 2000-2019
 !
 !
@@ -25 +25 @@
 //////////////////////////////////////////////////////////////////////////////
 //
-//  MSimAtmosphere
-//
-//  Task to calculate wavelength or incident angle dependent absorption
-//
-//  Input Containers:
-//   MPhotonEvent
-//   MCorsikaRunHeader
-//
-//  Output Containers:
-//   MPhotonEvent
+//  MAtmosphere
+//
+//  Class to calculate atmospheric absorption.
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -137 +130 @@
 
 MAtmRayleigh::MAtmRayleigh() : fR(MCorsikaRunHeader::fgEarthRadius),
-    //fHeight{0, 300000, 1e+06, 5e+06, 1e+07},
-    //fAtmA{-144.838, -124.071, 0.360027, -0.000824761, 0.00221589},
-    fAtmB{1192.34, 1173.98, 1412.08, 810.682, 1},
-    fAtmC{994186, 967530, 636143, 735640, 4.92961e9},
+    // U.S. Standard Atmosphere (after Keilhauer)
+    fHeight{0, 7.0e5, 11.4e5, 37.0e5, 100e5},
+    //fAtmA{-149.801663, -57.932486, 0.63631894, 4.35453690e-4},
+    fAtmB{1183.6071, 1143.0425, 1322.9748, 655.67307},
+    fAtmC{954248.34, 800005.34, 629568.93, 737521.77},
     fObsLevel(-1)
+
 {
 }
@@ -165 +160 @@
     // the observer
 
-    // FIXME: Could be replaced by 0, AtmLay[0]-fAtmLay[3]
-
-    const Double_t h[5] =
-    {
-        fObsLevel,                     // fObsLevel,                   //   0km
-        TMath::Max(fObsLevel, 7.75e5), // TMath::Max(fObsLevel, 4e5),  //   4km
-        16.5e5,                        //  10e5,                       //  10km
-        50.0e5,                        //  40e5,                       //  40km
-        105.0e5,                       // 100e5                        // 100km
-    };
-
-    memcpy(fHeight, h, sizeof(Double_t)*5);
-
     fRho[0] = 0;
     for (int i=0; i<4; i++)
     {
-        const Double_t b = fAtmB[i];
-        const Double_t c = fAtmC[i];
-
-        const Double_t h1 = h[i+1];
-        const Double_t h0 = h[i];
+        // Below the observation level, the atmospheric overburden is 0
+        const Float_t &h1 = fHeight[i+1];
+        if (h1<=fObsLevel)
+            fRho[i] = 0;
+
+        // Start integration of the overburden at fObsLevel
+        const Float_t &h0 = TMath::Max(fHeight[i], fObsLevel);
+
+        const Float_t &b = fAtmB[i];
+        const Float_t &c = fAtmC[i];
 
         const Double_t B = b*TMath::Exp(-h0/c);
@@ -198 +185 @@
 }
 
-void MAtmRayleigh::Init(Double_t obs, const Float_t *atmb, const Float_t *atmc)
+void MAtmRayleigh::Init(Float_t obs)
 {
     // Observation level above earth radius
     fObsLevel = obs;
-
-    //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
-    if (atmb)
-        memcpy(fAtmB, atmb, sizeof(Float_t)*5);
-    if (atmc)
-        memcpy(fAtmC, atmc, sizeof(Float_t)*5);
-
     PreCalcRho();
+}
+
+void MAtmRayleigh::Init(Float_t obs, const Float_t *atmb, const Float_t *atmc, const Float_t *height)
+{
+    memcpy(fAtmB,   atmb,   sizeof(Float_t)*4);
+    memcpy(fAtmC,   atmc,   sizeof(Float_t)*4);
+    memcpy(fHeight, height, sizeof(Float_t)*5);
+
+    Init(obs);
 }
 
@@ -220 +209 @@
     fR = h.EarthRadius();
 
-    Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC());
+    //memcpy(fAtmA, h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
+    Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC(), h.GetAtmosphericLayers());
 }
 
@@ -228 +218 @@
 Double_t MAtmRayleigh::GetVerticalThickness(Double_t height) const
 {
+    if (height<=fObsLevel)
+        return 0;
+
     // FIXME: We could store the start point above obs-level
     //        (Does this really gain anything?)
@@ -237 +230 @@
     const Double_t c = fAtmC[i];
 
+    // fRho is the integral up to the lower bound of the layer or the
+    // observation level is the obervation level is within the bin
     return fRho[i] - b*TMath::Exp(-height/c);
 }