Changeset 19852 for trunk/Mars


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/msim/MAtmosphere.cc

    r19765 r19852  
    1818!   Author(s): Thomas Bretz,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: CheObs Software Development, 2000-2009
     20!   Copyright: CheObs Software Development, 2000-2019
    2121!
    2222!
     
    2525//////////////////////////////////////////////////////////////////////////////
    2626//
    27 //  MSimAtmosphere
    28 //
    29 //  Task to calculate wavelength or incident angle dependent absorption
    30 //
    31 //  Input Containers:
    32 //   MPhotonEvent
    33 //   MCorsikaRunHeader
    34 //
    35 //  Output Containers:
    36 //   MPhotonEvent
     27//  MAtmosphere
     28//
     29//  Class to calculate atmospheric absorption.
    3730//
    3831//////////////////////////////////////////////////////////////////////////////
     
    137130
    138131MAtmRayleigh::MAtmRayleigh() : fR(MCorsikaRunHeader::fgEarthRadius),
    139     //fHeight{0, 300000, 1e+06, 5e+06, 1e+07},
    140     //fAtmA{-144.838, -124.071, 0.360027, -0.000824761, 0.00221589},
    141     fAtmB{1192.34, 1173.98, 1412.08, 810.682, 1},
    142     fAtmC{994186, 967530, 636143, 735640, 4.92961e9},
     132    // U.S. Standard Atmosphere (after Keilhauer)
     133    fHeight{0, 7.0e5, 11.4e5, 37.0e5, 100e5},
     134    //fAtmA{-149.801663, -57.932486, 0.63631894, 4.35453690e-4},
     135    fAtmB{1183.6071, 1143.0425, 1322.9748, 655.67307},
     136    fAtmC{954248.34, 800005.34, 629568.93, 737521.77},
    143137    fObsLevel(-1)
     138
    144139{
    145140}
     
    165160    // the observer
    166161
    167     // FIXME: Could be replaced by 0, AtmLay[0]-fAtmLay[3]
    168 
    169     const Double_t h[5] =
    170     {
    171         fObsLevel,                     // fObsLevel,                   //   0km
    172         TMath::Max(fObsLevel, 7.75e5), // TMath::Max(fObsLevel, 4e5),  //   4km
    173         16.5e5,                        //  10e5,                       //  10km
    174         50.0e5,                        //  40e5,                       //  40km
    175         105.0e5,                       // 100e5                        // 100km
    176     };
    177 
    178     memcpy(fHeight, h, sizeof(Double_t)*5);
    179 
    180162    fRho[0] = 0;
    181163    for (int i=0; i<4; i++)
    182164    {
    183         const Double_t b = fAtmB[i];
    184         const Double_t c = fAtmC[i];
    185 
    186         const Double_t h1 = h[i+1];
    187         const Double_t h0 = h[i];
     165        // Below the observation level, the atmospheric overburden is 0
     166        const Float_t &h1 = fHeight[i+1];
     167        if (h1<=fObsLevel)
     168            fRho[i] = 0;
     169
     170        // Start integration of the overburden at fObsLevel
     171        const Float_t &h0 = TMath::Max(fHeight[i], fObsLevel);
     172
     173        const Float_t &b = fAtmB[i];
     174        const Float_t &c = fAtmC[i];
    188175
    189176        const Double_t B = b*TMath::Exp(-h0/c);
     
    198185}
    199186
    200 void MAtmRayleigh::Init(Double_t obs, const Float_t *atmb, const Float_t *atmc)
     187void MAtmRayleigh::Init(Float_t obs)
    201188{
    202189    // Observation level above earth radius
    203190    fObsLevel = obs;
    204 
    205     //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
    206     if (atmb)
    207         memcpy(fAtmB, atmb, sizeof(Float_t)*5);
    208     if (atmc)
    209         memcpy(fAtmC, atmc, sizeof(Float_t)*5);
    210 
    211191    PreCalcRho();
     192}
     193
     194void MAtmRayleigh::Init(Float_t obs, const Float_t *atmb, const Float_t *atmc, const Float_t *height)
     195{
     196    memcpy(fAtmB,   atmb,   sizeof(Float_t)*4);
     197    memcpy(fAtmC,   atmc,   sizeof(Float_t)*4);
     198    memcpy(fHeight, height, sizeof(Float_t)*5);
     199
     200    Init(obs);
    212201}
    213202
     
    220209    fR = h.EarthRadius();
    221210
    222     Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC());
     211    //memcpy(fAtmA, h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
     212    Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC(), h.GetAtmosphericLayers());
    223213}
    224214
     
    228218Double_t MAtmRayleigh::GetVerticalThickness(Double_t height) const
    229219{
     220    if (height<=fObsLevel)
     221        return 0;
     222
    230223    // FIXME: We could store the start point above obs-level
    231224    //        (Does this really gain anything?)
     
    237230    const Double_t c = fAtmC[i];
    238231
     232    // fRho is the integral up to the lower bound of the layer or the
     233    // observation level is the obervation level is within the bin
    239234    return fRho[i] - b*TMath::Exp(-height/c);
    240235}
  • trunk/Mars/msim/MAtmosphere.h

    r19766 r19852  
    1616    Double_t fR;         // [cm] Earth radius to be used
    1717
    18     Double_t fHeight[5]; // Layer boundaries (atmospheric layer)
     18    Float_t fHeight[5]; // Layer boundaries (atmospheric layer)
    1919
    2020    // Parameters of the different atmospheres. We use the same parametrization
    2121    // shape as in Corsika atmospheric models (see Corsika manual, appendix D).
    2222    // The values here can be/are obtained from the Corsika output
    23     //Float_t  fAtmA[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
    24     Float_t  fAtmB[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
    25     Float_t  fAtmC[5];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
     23    // fAtmA is not required as it is an offset which cancels out
     24    //Float_t  fAtmA[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
     25    Float_t  fAtmB[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
     26    Float_t  fAtmC[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
    2627
    2728    Double_t fRho[5];    // Precalculated integrals for rayleigh scatterning
     
    3031
    3132protected:
    32     Double_t fObsLevel; // [cm] observation level a.s.l.
     33    Float_t fObsLevel; // [cm] observation level a.s.l.
    3334
    3435public:
     
    4950    Double_t R() const { return fR; }
    5051
    51     void Init(Double_t obs, const Float_t *atmb=0, const Float_t *atmc=0);
     52    void Init(Float_t obs);
     53    void Init(Float_t obs, const Float_t *atmb, const Float_t *atmc, const Float_t *height);
    5254    void Init(const MCorsikaRunHeader &h);
    5355
    5456    Double_t GetVerticalThickness(Double_t height) const;
    5557    Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const;
     58
     59    static Double_t Fcn(const Double_t *xx, const Double_t *par)
     60    {
     61        const Double_t *A = par;
     62        const Double_t *B = par+5;
     63        const Double_t *C = par+10;
     64
     65        double H[6] = { 0, 0, 0, 0, 0, 1e100 };
     66
     67        const Double_t x = xx[0]*1e5; // km -> cm
     68
     69        for (int i=0; i<4; i++)
     70        {
     71            H[i+1] = log(C[i+1]/C[i] * B[i]/B[i+1]) / (1./C[i] - 1./C[i+1]);
     72
     73            if (x>=H[i] && x<H[i+1])
     74                return log10(A[i] + B[i] * exp(-x/C[i]));
     75        }
     76
     77        return -1;
     78    }
    5679};
    5780
     
    98121        Init(h);//, "ozone.txt", "aerosols.txt");
    99122    }
    100     MAtmosphere(Double_t obs=0, const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
     123    MAtmosphere(Float_t obs=0, const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
    101124    {
    102125        MAtmRayleigh::Init(obs);
Note: See TracChangeset for help on using the changeset viewer.