Index: /trunk/Mars/msim/MAtmosphere.cc =================================================================== --- /trunk/Mars/msim/MAtmosphere.cc (revision 19851) +++ /trunk/Mars/msim/MAtmosphere.cc (revision 19852) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz, 1/2009 ! -! Copyright: CheObs Software Development, 2000-2009 +! Copyright: CheObs Software Development, 2000-2019 ! ! @@ -25,14 +25,7 @@ ////////////////////////////////////////////////////////////////////////////// // -// MSimAtmosphere -// -// Task to calculate wavelength or incident angle dependent absorption -// -// Input Containers: -// MPhotonEvent -// MCorsikaRunHeader -// -// Output Containers: -// MPhotonEvent +// MAtmosphere +// +// Class to calculate atmospheric absorption. // ////////////////////////////////////////////////////////////////////////////// @@ -137,9 +130,11 @@ 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,25 +160,17 @@ // 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,16 +185,18 @@ } -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,5 +209,6 @@ 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,4 +218,7 @@ 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,4 +230,6 @@ 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); } Index: /trunk/Mars/msim/MAtmosphere.h =================================================================== --- /trunk/Mars/msim/MAtmosphere.h (revision 19851) +++ /trunk/Mars/msim/MAtmosphere.h (revision 19852) @@ -16,12 +16,13 @@ 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,5 +31,5 @@ protected: - Double_t fObsLevel; // [cm] observation level a.s.l. + Float_t fObsLevel; // [cm] observation level a.s.l. public: @@ -49,9 +50,31 @@ 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