Ignore:
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!
Location:
trunk/MagicSoft/Simulation/Detector/Reflector
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Simulation/Detector/Reflector/Makefile

    r348 r428  
    1919#
    2020# $RCSfile: Makefile,v $
    21 # $Revision: 1.4 $
     21# $Revision: 1.5 $
    2222# $Author: harald $
    23 # $Date: 2000-01-28 09:19:54 $
     23# $Date: 2000-07-25 13:48:09 $
    2424#
    2525##################################################################
     
    4040           -I${INCLUDE_MC}   \
    4141           -I${INCLUDE_REFL} \
    42            -I/usr/include  -I/usr/include/g++
     42           -I/usr/include  -I/usr/include/cxx
    4343
    4444#CERNLIBDIR = ${CERNDIR}/pro/lib/
     
    4949# special flags
    5050
    51 osf_FORLIBS = -lUfor -lfor -lutil -lots -lm
     51osf_FORLIBS = -L/usr/lib -lUfor -lFutil -lfor -lutil -lots -lm
    5252linux_FORLIBS = -lm
    5353generic_FORLIBS = -lm
     
    6161FFLAGS    = ${CXXFLAGS}
    6262LIBS      = ${RANLIB} ${FORLIBS}
    63 #LIBS      = ${CERNLIB} ${RANLIB} ${FORLIBS}
     63LIBS      = ${CERNLIB} ${RANLIB} ${FORLIBS}
    6464
    6565#------------------------------------------------------------------------------
     
    100100        attenu.o \
    101101        readparam.o \
    102         reflector.o     
     102        reflector.o
    103103
    104104PROGRAM=reflector
  • trunk/MagicSoft/Simulation/Detector/Reflector/attenu.f

    r407 r428  
    3737     + RHO_TOT, RHO_FI, RHOFP
    3838     
     39      DOUBLE PRECISION Rcos2, Rsin2
     40      DOUBLE PRECISION x1, x2, x3, x4
     41      DOUBLE PRECISION e1, e2, e3, e4
     42
    3943      REAL wavelength, height, theta, tr_atmos
    4044      real trr, trm, tro
     
    278282      DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0,0.D0  /
    279283      DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0,1.D-9/
    280 c fs 190400: lowest Layer is 2200.d0 (observation level)
    281 c      DATA LAHG / 2000.0D2,4.0D5,1.0D6,4.0D6,1.0D7 / 
     284
    282285      DATA LAHG / 2200.0D2,4.0D5,1.0D6,4.0D6,1.0D7 /
    283286      DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
    284      + 5.99D5,6.33D5, 6.67D5,7.04D5 /
     287     +         5.99D5,6.33D5, 6.67D5,7.04D5 /
    285288      DATA LONG /
    286289     +    280.0,
     
    297300     +    650.0 /
    298301      DATA PI / 3.141592654D0 /
    299 c fs 191999: take same Earth radius as in CORSIKA
    300 c      PARAMETER (RT=6348.0D5)
     302
     303c--   Take same Earth radius as in CORSIKA
    301304      parameter (rt = 6371315.D2)
     305
     306c--   Scale-height for Exponential density profile
     307      parameter (hscale = 7.4d5)
    302308
    303309***********************************************************************
     
    309315***********************************************************************
    310316       
    311 c fs observation level at La Palma
     317c-- Observation level at La Palma
    312318      parameter (obslev = 2200.d2)
    313319     
     
    316322      T_Oz = 1.0
    317323     
    318 c fs : error: height with obslev = 0. !!!!
     324c-- Height calculated using an obslev = H(LaPalma) <> 0
    319325c      H = -RT + SQRT(RT**2 + (height/COS(theta))**2 +
    320326c     +     (2.0D0*RT*height))
     
    333339c      => pure geometric correction on
    334340c         absorption lenght measured at vertical height [km^-1]   
    335 c      m= (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
     341c--
     342c      ma = (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
    336343c     +     (RT*COS(theta)))
    337       m = (1.d0/(h-obslev)) *
    338      &   ( -(rt+obslev)*cos(theta)
    339      &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
     344c--
     345c      mb = (1.d0/(h-obslev)) *
     346c     &     ( -(rt+obslev)*cos(theta)
     347c     &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
     348c--
     349
     350c-- Air mass "m" calcualted using a one-exponential density
     351c-- profile for the atmosphere, rho = rho_0 exp(-h/Hs)
     352c-- with Hs = hscale = 7.4 km
     353
     354      Rcos2 = rt * cos(theta)**2
     355      Rsin2 = rt * sin(theta)**2
     356     
     357      x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
     358      x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
     359      x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
     360      x4 = sqrt((2. * h      + rt)    / (2. * hscale))
     361
     362      e1 = derfc(x1)
     363      e2 = derfc(x2)
     364      e3 = derfc(x3)
     365      e4 = derfc(x4)
     366
     367      m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))
     368
    340369**********************************************************************   
    341          
    342 
     370*         
    343371*     RAYLEIGH SCATTERING
    344372
Note: See TracChangeset for help on using the changeset viewer.