Changeset 19852 for trunk/Mars
- Timestamp:
- 11/01/19 15:11:07 (5 years ago)
- Location:
- trunk/Mars/msim
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/msim/MAtmosphere.cc
r19765 r19852 18 18 ! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de> 19 19 ! 20 ! Copyright: CheObs Software Development, 2000-20 0920 ! Copyright: CheObs Software Development, 2000-2019 21 21 ! 22 22 ! … … 25 25 ////////////////////////////////////////////////////////////////////////////// 26 26 // 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. 37 30 // 38 31 ////////////////////////////////////////////////////////////////////////////// … … 137 130 138 131 MAtmRayleigh::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}, 143 137 fObsLevel(-1) 138 144 139 { 145 140 } … … 165 160 // the observer 166 161 167 // FIXME: Could be replaced by 0, AtmLay[0]-fAtmLay[3]168 169 const Double_t h[5] =170 {171 fObsLevel, // fObsLevel, // 0km172 TMath::Max(fObsLevel, 7.75e5), // TMath::Max(fObsLevel, 4e5), // 4km173 16.5e5, // 10e5, // 10km174 50.0e5, // 40e5, // 40km175 105.0e5, // 100e5 // 100km176 };177 178 memcpy(fHeight, h, sizeof(Double_t)*5);179 180 162 fRho[0] = 0; 181 163 for (int i=0; i<4; i++) 182 164 { 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]; 188 175 189 176 const Double_t B = b*TMath::Exp(-h0/c); … … 198 185 } 199 186 200 void MAtmRayleigh::Init( Double_t obs, const Float_t *atmb, const Float_t *atmc)187 void MAtmRayleigh::Init(Float_t obs) 201 188 { 202 189 // Observation level above earth radius 203 190 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 211 191 PreCalcRho(); 192 } 193 194 void 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); 212 201 } 213 202 … … 220 209 fR = h.EarthRadius(); 221 210 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()); 223 213 } 224 214 … … 228 218 Double_t MAtmRayleigh::GetVerticalThickness(Double_t height) const 229 219 { 220 if (height<=fObsLevel) 221 return 0; 222 230 223 // FIXME: We could store the start point above obs-level 231 224 // (Does this really gain anything?) … … 237 230 const Double_t c = fAtmC[i]; 238 231 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 239 234 return fRho[i] - b*TMath::Exp(-height/c); 240 235 } -
trunk/Mars/msim/MAtmosphere.h
r19766 r19852 16 16 Double_t fR; // [cm] Earth radius to be used 17 17 18 Double_t fHeight[5]; // Layer boundaries (atmospheric layer)18 Float_t fHeight[5]; // Layer boundaries (atmospheric layer) 19 19 20 20 // Parameters of the different atmospheres. We use the same parametrization 21 21 // shape as in Corsika atmospheric models (see Corsika manual, appendix D). 22 22 // 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) 26 27 27 28 Double_t fRho[5]; // Precalculated integrals for rayleigh scatterning … … 30 31 31 32 protected: 32 Double_t fObsLevel; // [cm] observation level a.s.l.33 Float_t fObsLevel; // [cm] observation level a.s.l. 33 34 34 35 public: … … 49 50 Double_t R() const { return fR; } 50 51 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); 52 54 void Init(const MCorsikaRunHeader &h); 53 55 54 56 Double_t GetVerticalThickness(Double_t height) const; 55 57 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 } 56 79 }; 57 80 … … 98 121 Init(h);//, "ozone.txt", "aerosols.txt"); 99 122 } 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) 101 124 { 102 125 MAtmRayleigh::Init(obs);
Note:
See TracChangeset
for help on using the changeset viewer.