CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE calc_rez(LAMBDA, L, EPSILON, N_M, P, M, N_RL, B, DK)
      IMPLICIT NONE
      INCLUDE "lun.par"
      REAL*4 LAMBDA, L, EPSILON, N_M, P, M, N_RL, B
      REAL*4 E, GAMMA, BETA
      REAL*4 L_PRIME, DK_RES, PT, DK_MS, DK, K, DP_OVER_P, N_RL_TOTAL
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      L_prime = L*cos(lambda)
      dk_res = epsilon/L_prime**2*sqrt(720/(n_m + 4))
      pt = p*cos(lambda)
      N_RL_TOTAL = N_RL/COS(LAMBDA)
      dk_ms = 0.016/(L*p*beta*cos(lambda)**2)*sqrt(n_rl_TOTAL)
      dk = sqrt(dk_res**2 + dk_ms**2)
      k = 0.3*B/pt
      dp_over_p = dk/k
      WRITE (LUNOUT,*) "p ", p, " dp/p ", dp_over_p
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CALC_AZIMUTH_ANGLE(LAMBDA, L, EPSILON, N_M, P, M, N_RL,
     X     N_RL_FRONT, B, R_M, DK)
      IMPLICIT NONE
      INCLUDE "lun.par"
      REAL*4 LAMBDA, L, EPSILON, N_M, P, M, N_RL, N_RL_FRONT, B, R_M,
     X     DK
      REAL*4 E, GAMMA, BETA
      REAL*4 THETA0, C, DTHETA_MS, DTHETA_RES, DTHETA_K, DTHETA,
     X     N_RL_TOTAL
      DATA C/1.0/
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DTHETA_MS = THETA0/SQRT(3.0)
      DTHETA_RES = EPSILON/L*SQRT(12*(N_M - 1)/(N_M*(N_M + 1)))
      DTHETA_K = R_M*DK
      DTHETA = SQRT(DTHETA_RES**2 + DTHETA_MS**2 + DTHETA_K**2)
      WRITE (LUNOUTA,*) "p ", P, " dtheta ", DTHETA
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CALC_POLAR_ANGLE(LAMBDA, L, EPSILON, STEREO, N_M_POLAR,
     X     P, M, N_RL, N_RL_FRONT, B)
      IMPLICIT NONE
      INCLUDE "lun.par"
      REAL*4 LAMBDA, L, EPSILON, STEREO, N_M_POLAR, P, M, N_RL,
     X     N_RL_FRONT, B
      REAL*4 E, GAMMA, BETA
      REAL*4 C, N_RL_TOTAL, THETA0, DTHETA_MS, DTHETA_RES, DTHETA
      DATA C/1.0/
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DTHETA_MS = THETA0/SQRT(3.0)
      DTHETA_RES = EPSILON*COS(LAMBDA)/(SIN(STEREO)*L)
     X     *SQRT(12*(N_M_POLAR - 1)/(N_M_POLAR*(N_M_POLAR + 1)))
      DTHETA = SQRT(DTHETA_RES**2 + DTHETA_MS**2)
      WRITE (LUNOUTB,*) "p ", P, " dtheta ", DTHETA
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE KINEMATICS(P, M, E, GAMMA, BETA)
      REAL*4 P, M
      REAL*4 E, GAMMA, BETA
      REAL*4 C
      DATA C/1.0/
      E = sqrt(p**2 + m**2)
      gamma = E/m
      beta = sqrt(1.0 - 1/gamma**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC