Changeset 19852 for trunk

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.

Location:

trunk/Mars/msim

Files:

• MAtmosphere.cc (modified) (8 diffs)
• MAtmosphere.h (modified) (4 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);
 }

trunk/Mars/msim/MAtmosphere.h

@@ -16 +16 @@
     Double_t fR;         // [cm] Earth radius to be used
 
-    Double_t fHeight[5]; // Layer boundaries (atmospheric layer)
+    Float_t fHeight[5]; // Layer boundaries (atmospheric layer)
 
     // Parameters of the different atmospheres. We use the same parametrization
     // shape as in Corsika atmospheric models (see Corsika manual, appendix D).
     // The values here can be/are obtained from the Corsika output
-    //Float_t  fAtmA[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
-    Float_t  fAtmB[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
-    Float_t  fAtmC[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
+    // fAtmA is not required as it is an offset which cancels out
+    //Float_t  fAtmA[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
+    Float_t  fAtmB[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
+    Float_t  fAtmC[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
 
     Double_t fRho[5];    // Precalculated integrals for rayleigh scatterning
@@ -30 +31 @@
 
 protected:
-    Double_t fObsLevel; // [cm] observation level a.s.l.
+    Float_t fObsLevel; // [cm] observation level a.s.l.
 
 public:
@@ -49 +50 @@
     Double_t R() const { return fR; }
 
-    void Init(Double_t obs, const Float_t *atmb=0, const Float_t *atmc=0);
+    void Init(Float_t obs);
+    void Init(Float_t obs, const Float_t *atmb, const Float_t *atmc, const Float_t *height);
     void Init(const MCorsikaRunHeader &h);
 
     Double_t GetVerticalThickness(Double_t height) const;
     Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const;
+
+    static Double_t Fcn(const Double_t *xx, const Double_t *par)
+    {
+        const Double_t *A = par;
+        const Double_t *B = par+5;
+        const Double_t *C = par+10;
+
+        double H[6] = { 0, 0, 0, 0, 0, 1e100 };
+
+        const Double_t x = xx[0]*1e5; // km -> cm
+
+        for (int i=0; i<4; i++)
+        {
+            H[i+1] = log(C[i+1]/C[i] * B[i]/B[i+1]) / (1./C[i] - 1./C[i+1]);
+
+            if (x>=H[i] && x<H[i+1])
+                return log10(A[i] + B[i] * exp(-x/C[i]));
+        }
+
+        return -1;
+    }
 };
 
@@ -98 +121 @@
         Init(h);//, "ozone.txt", "aerosols.txt");
     }
-    MAtmosphere(Double_t obs=0, const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
+    MAtmosphere(Float_t obs=0, const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
     {
         MAtmRayleigh::Init(obs);