Changeset 19763

Timestamp:

10/12/19 15:48:08 (5 years ago)

Author:

tbretz

Message:

Moved the atmosphere code to its own files and make it accessible from teh interpreter.

Location:

trunk/Mars/msim

Files:

• MAtmosphere.cc (added)
• MAtmosphere.h (added)
• MSimAtmosphere.cc (modified) (10 diffs)
• Makefile (modified) (1 diff)
• SimLinkDef.h (modified) (1 diff)

trunk/Mars/msim/MSimAtmosphere.cc

@@ -53 +53 @@
 #include "MPhotonData.h"
 
+#include "MAtmosphere.h"
+
 ClassImp(MSimAtmosphere);
 
 using namespace std;
-
+/*
 // ==========================================================================
 //
@@ -253 +255 @@
     }
 
-    /*
-     // The "orginal" code for the vertical thickness
-    Double_t GetVerticalThickness(Double_t obslev, Double_t height) const
-    {
-        // This is a C++-version of the original code from attenu.c
-
-        // Limits (height in cm) of the four layers in which atmosphere is parametrized:
-        const double lahg[5] =
-        {
-            obslev,
-            TMath::Max(obslev, 4e5),
-            1.0e6,
-            4.0e6,
-            1.0e7
-        };
-
-        Double_t Rho_Tot = 0.0;
-        for (int i=0; i<4; i++)
-        {
-            const Double_t b = fAtmB[i];
-            const Double_t c = fAtmC[i];
-
-            const Double_t h0 = TMath::Min(height, lahg[i+1]);
-            const Double_t h1 =                    lahg[i];
-
-            // Calculate rho for the i-th layer from the lower
-            // to the higher layer boundary.
-            // If height is within the layer only calculate up to height.
-            Rho_Tot += b*(exp(-h1/c) - exp(-h0/c));
-
-            if (lahg[i+1] > height)
-                break;
-        }
-
-        return Rho_Tot;
-    }
-    */
     Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const
     {
@@ -334 +299 @@
         // Calculate the traversed "vertical thickness" of air using the
         // US Standard atmosphere:
-        const Double_t Rho_tot = GetVerticalThickness(/*fObsLevel,*/ height);
+        const Double_t Rho_tot = GetVerticalThickness(height);
 
         // We now convert from "vertical thickness" to "slanted thickness"
@@ -387 +352 @@
         // checks are done for each call!
         return g.GetN()==0 ? 0 : g.Eval(wavelength)*1e-5; // from km^-1 to cm^-1
-/*
-        // Linear interpolation
-        // (FIXME: Is it faster to be replaced with a binary search?)
-        // (       This might be faster because we have more photons
-        //         with smaller wavelengths)
-        //int index;
-        //for (index = 1; index <g.GetN()-1; index++)
-        //    if (wavelength < g.GetX()[index])
-        //        break;
-        const Int_t index = TMath::BinarySearch(g.GetN(), g.GetX(), wavelength)+1;
-
-        const Double_t t0 = g.GetY()[index-1];
-        const Double_t t1 = g.GetY()[index];
-
-        const Double_t w0 = g.GetX()[index-1];
-        const Double_t w1 = g.GetX()[index];
-
-        const Double_t beta0 = t0+(t1-t0)*(wavelength-w0)/(w1-w0);
-
-        return beta0 * 1e-5; // from km^-1 to cm^-1
-        */
     }
 
@@ -490 +434 @@
         //const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree
 
-        /* Mie (aerosol): */
+        // Mie (aerosol):
         for (Int_t j = 0; j < 90; j++)
         {
@@ -566 +510 @@
         InitAerosols(name2);
     }
-/*
-    Double_t GetOz(Double_t height, Double_t theta) const
-    {
-        // Distance between two points D = 1km /cos(theta)
-        // Density along y within this km:   f =  (x[i+1]-x[i])/1km * dy
-        // Integral of this density  f =  (x[i+1]-x[i])/1km * (y[i+1]-y[i])
-        // f(h) = int [ (c1-c0)/1km*(h-h0)*dh + c0 ] dh
-        //      = (c-co)*(h-h0)
-
-        Double_t rc = 0;
-        int i;
-        for (i=0; i<49; i++)
-            if (i>=2 && i+1<height/1e5)    // cm -> km
-                rc += oz_conc[i] * 1e5/cos(theta);
-
-        rc -= oz_conc[2]*0.2*1e5/cos(theta);
-        rc += oz_conc[i+1]*fmod(height/1e5,1)*1e5/cos(theta);
-
-        return rc;
-    }
-    */
 
     Double_t CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const
@@ -593 +516 @@
             return 1;
 
-        //******* Ozone absorption *******
+        // ******* Ozone absorption *******
         if (h > 50.e5)
             h = 50.e5;
@@ -617 +540 @@
             return 1;
 
-        //******* Mie (aerosol) *******
+        // ******* Mie (aerosol) *******
         if (h > 30.e5)
             h = 30.e5;
@@ -665 +588 @@
         const Double_t h = TMath::Sqrt(H*H + d*d + 2*H*z) - R();
 
-        //**** Rayleigh scattering: *****
+        // **** Rayleigh scattering: *****
         const Double_t T_Ray = CalcTransmission(h, wavelength, sin2);
         if (T_Ray<0)
             return 0;
 
-        //****** Ozone absorption: ******
+        // ****** Ozone absorption: ******
         const Double_t T_Oz  = CalcOzoneAbsorption(h, wavelength, theta);
 
-        //******** Mie (aerosol) ********
+        // ******** Mie (aerosol) ********
         const Double_t T_Mie = CalcAerosolAbsorption(h, wavelength, theta);
 
@@ -690 +613 @@
 
 const Double_t MAtmosphere::oz_conc[51]={0.3556603E-02, 0.3264150E-02, 0.2933961E-02, 0.2499999E-02, 0.2264150E-02, 0.2207546E-02, 0.2160377E-02, 0.2226414E-02, 0.2283018E-02, 0.2811320E-02, 0.3499999E-02, 0.4603772E-02, 0.6207545E-02, 0.8452828E-02, 0.9528299E-02, 0.9905657E-02, 0.1028302E-01, 0.1113207E-01, 0.1216981E-01, 0.1424528E-01, 0.1641509E-01, 0.1839622E-01, 0.1971697E-01, 0.1981131E-01, 0.1933962E-01, 0.1801886E-01, 0.1632075E-01, 0.1405660E-01, 0.1226415E-01, 0.1066037E-01, 0.9028300E-02, 0.7933960E-02, 0.6830187E-02, 0.5820753E-02, 0.4830188E-02, 0.4311319E-02, 0.3613206E-02, 0.3018867E-02, 0.2528301E-02, 0.2169811E-02, 0.1858490E-02, 0.1518867E-02, 0.1188679E-02, 0.9301884E-03, 0.7443394E-03, 0.5764149E-03, 0.4462263E-03, 0.3528301E-03, 0.2792452E-03, 0.2226415E-03, 0.1858490E-03};
-
+*/
 // ==========================================================================
 

trunk/Mars/msim/Makefile

@@ -29 +29 @@
            MSimAbsorption.cc \
            MSimPointingPos.cc \
-           MClonesArray.cc
+           MClonesArray.cc \
+           MAtmosphere.cc
 
 ############################################################

trunk/Mars/msim/SimLinkDef.h

@@ -20 +20 @@
 #pragma link C++ class MClonesArray-; // - needed as custom streamer is already implemented
 
+#pragma link C++ class MAtmosphere+;
+#pragma link C++ class MAtmRayleigh+;
+
 #endif