Changeset 428 for trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f

Timestamp:

07/25/00 14:48:09 (24 years ago)

Author:

harald

Message:

During the Monte Carlo Software meeting, a new version of attenu.f was
promised by Jose Carlos, Frank Schroeder and Aitor Ibarra. I got this
version for the repository. A small change in the Makefile was also
neccessary. Thanks to Jose Carlos for his support to run it on alpha!

File:

• trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f (modified) (6 diffs)

trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f

@@ -37 +37 @@
      + RHO_TOT, RHO_FI, RHOFP
       
+      DOUBLE PRECISION Rcos2, Rsin2
+      DOUBLE PRECISION x1, x2, x3, x4
+      DOUBLE PRECISION e1, e2, e3, e4
+
       REAL wavelength, height, theta, tr_atmos
       real trr, trm, tro
@@ -278 +282 @@
       DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0,0.D0  /
       DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0,1.D-9/
-c fs 190400: lowest Layer is 2200.d0 (observation level)
-c      DATA LAHG / 2000.0D2,4.0D5,1.0D6,4.0D6,1.0D7 /  
+
       DATA LAHG / 2200.0D2,4.0D5,1.0D6,4.0D6,1.0D7 /
       DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
-     + 5.99D5,6.33D5, 6.67D5,7.04D5 /
+     +         5.99D5,6.33D5, 6.67D5,7.04D5 /
       DATA LONG /
      +    280.0, 
@@ -297 +300 @@
      +    650.0 /
       DATA PI / 3.141592654D0 /
-c fs 191999: take same Earth radius as in CORSIKA
-c      PARAMETER (RT=6348.0D5)
+
+c--   Take same Earth radius as in CORSIKA
       parameter (rt = 6371315.D2)
+
+c--   Scale-height for Exponential density profile
+      parameter (hscale = 7.4d5)
 
 ***********************************************************************
@@ -309 +315 @@
 ***********************************************************************
         
-c fs observation level at La Palma
+c-- Observation level at La Palma
       parameter (obslev = 2200.d2)
       
@@ -316 +322 @@
       T_Oz = 1.0
       
-c fs : error: height with obslev = 0. !!!!
+c-- Height calculated using an obslev = H(LaPalma) <> 0
 c      H = -RT + SQRT(RT**2 + (height/COS(theta))**2 +
 c     +     (2.0D0*RT*height))
@@ -333 +339 @@
 c      => pure geometric correction on 
 c         absorption lenght measured at vertical height [km^-1]   
-c      m= (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
+c--
+c      ma = (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
 c     +     (RT*COS(theta)))
-      m = (1.d0/(h-obslev)) * 
-     &   ( -(rt+obslev)*cos(theta) 
-     &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
+c--
+c      mb = (1.d0/(h-obslev)) * 
+c     &     ( -(rt+obslev)*cos(theta) 
+c     &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
+c--
+
+c-- Air mass "m" calcualted using a one-exponential density
+c-- profile for the atmosphere, rho = rho_0 exp(-h/Hs)
+c-- with Hs = hscale = 7.4 km
+
+      Rcos2 = rt * cos(theta)**2
+      Rsin2 = rt * sin(theta)**2
+      
+      x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
+      x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
+      x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
+      x4 = sqrt((2. * h      + rt)    / (2. * hscale))
+
+      e1 = derfc(x1)
+      e2 = derfc(x2)
+      e3 = derfc(x3)
+      e4 = derfc(x4)
+
+      m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))
+
 **********************************************************************   
-          
-
+*          
 *     RAYLEIGH SCATTERING