source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/externals/slast/src/kernel/firstpoint_hadron.F @ 117

      SUBROUTINE FIRSTPOINT_HADRON(Lmax)
      IMPLICIT NONE
#include "event.inc"
#include "detector.inc"
#include "erandom.inc"
#include "constants.inc"
#include "book_histos.inc"
#include "cwntvars.inc"
 
      REAL zmin,z_rand,r_rand,x_in_cone,y_in_cone,z_in_cone
      REAL H_R,THETA_R,LENGTH12,SCALAR11,SCALAR12,SCALAR22
      REAL dir_unit(3),R2(3)
      REAL THETAMAX,Lmax
      REAL Xinteraction,HINT, ATM
      REAL*8 RNDM,aran,bran,cran,r2ran,rinv
      LOGICAL bFirstPoint
*
      zmin = ISS(3)*SIN(FOV)**2 + 
     +     COS(FOV)*SQRT(R_EARTH**2-ISS(3)**2*SIN(FOV)**2)
      bFirstPoint = .FALSE.
      do while (.NOT.bFirstPoint)      
*     generate point on the Earth surface
        aran=0.
        bran=0.
        cran=0.
        r2ran=1. 
        do while(r2ran.gt.0.25)
            aran = ESAFRNDM() -0.5
            bran = ESAFRNDM() -0.5
            cran = ESAFRNDM() -0.5
            r2ran = aran*aran+bran*bran+cran*cran
        enddo
        rinv = R_EARTH/SQRT(r2ran)
        x_in_cone = aran*rinv
        y_in_cone = bran*rinv
        z_in_cone = cran*rinv
        if (z_in_cone.GT.Zmin) bFirstPoint = .TRUE.
      enddo
      z_rand = z_in_cone
      r_rand = x_in_cone**2 + y_in_cone**2
      
      PI = ACOS(-1.)
      RNDM = ESAFRNDM()
      PHI = PHIR(1) + (PHIR(2)-PHIR(1))*RNDM
      H_R = SQRT(r_rand**2+z_rand**2) - R_EARTH
*     ... SIN(2THETA) RANDOM THETA  ...
      RNDM = ESAFRNDM()  
      THETA_R = 0.5*ACOS(COS(2*THETAR(1)/180*PI) - 
     +    RNDM*(COS(2*THETAR(1)/180*PI)-COS(2*THETAR(2)/180*PI)))
      THETA = THETA_R/PI*180
*     ... Find far first point
      dir_unit(1) = SIN(THETA_R)*COS(PHI/180.*PI)
      dir_unit(2) = SIN(THETA_R)*SIN(PHI/180.*PI)
      dir_unit(3) = COS(THETA_R)
      R2(1) = x_in_cone
      R2(2) = y_in_cone
      R2(3) = z_rand
      
      CALL VECTORDOT(3,dir_unit,R2,SCALAR12) 
      CALL VECTORDOT(3,R2,R2,SCALAR22)
      SCALAR22 = SQRT(SCALAR22)
      LENGTH12 = SQRT(80*(SCALAR22+80.) + SCALAR12**2) - SCALAR12
      RINT(1) = R2(1) + LENGTH12*SIN(THETA/180*PI)*COS(PHI/180.*PI)
      RINT(2) = R2(2) + LENGTH12*SIN(THETA/180*PI)*SIN(PHI/180.*PI)
      RINT(3) = R2(3) + LENGTH12*COS(THETA/180*PI) 

      Xinteraction = 0.
*     ... RANDOM FIRST INTERACTION DEPTH ...
      CALL First_Depth_Hadron
      do while (Xinteraction.LT.X1)
          RINT(1) = RINT(1) - 0.1*SIN(THETA/180*PI)*COS(PHI/180.*PI)
          RINT(2) = RINT(2) - 0.1*SIN(THETA/180*PI)*SIN(PHI/180.*PI)
          RINT(3) = RINT(3) - 0.1*COS(THETA/180*PI) 
          CALL VECTORDOT(3,RINT,RINT,SCALAR22)
          HINT = SQRT(SCALAR22) - R_EARTH
          Xinteraction = Xinteraction + 0.1*ATM(HINT)
      enddo    
*     ...  THETA,PHI ARE RETURNED IN DEGREES ...
      CALL ACCEPTANCE(Lmax,ACCEPTANCE_TYPE)
      RINT(3) = RINT(3) - R_EARTH
*
      END