fRE_H_3D Function

  FUNCTION fRE_H_3D(F,eta,g,R,Z,R0,Z0)
    TYPE(FIELDS), INTENT(IN)    :: F
    REAL(rp), INTENT(IN) 	:: eta ! pitch angle in degrees
    REAL(rp), INTENT(IN) 	:: g ! Relativistic gamma factor
    REAL(rp), INTENT(IN)  :: R,Z,R0,Z0
    REAL(rp) 				:: fRE_H_3D
    REAL(rp) 				:: D
    REAL(rp) 				:: g0
    REAL(rp) 				:: g1
    REAL(rp) 				:: f0
    REAL(rp) 				:: f1
    REAL(rp) 				:: m
    REAL(rp) 				:: feta
    REAL(rp) 				:: A
    REAL(rp) 				:: rm,E_G,Z_G
    INTEGER 				:: index

    index = MINLOC(ABS(h_params%g - g),1)
    ! index of gamma supplied to function in Hollmann input gamma range
    D = h_params%g(index) - g

    ! linear interpolation of Hollmann input gamma range to gamma supplied
    ! to function
    if (D.GT.0) then
       f0 = h_params%fRE_E(index-1)
       g0 = h_params%g(index-1)

       f1 = h_params%fRE_E(index)
       g1 = h_params%g(index)
    else
       f0 = h_params%fRE_E(index)
       g0 = h_params%g(index)

       f1 = h_params%fRE_E(index+1)
       g1 = h_params%g(index+1)
    end if

    m = (f1-f0)/(g1-g0)

    fRE_H_3D = f0 + m*(g - g0)
    ! end of linear interpolation, fRE_H is evaluation of input Hollman energy
    ! distribution PDF at gamma supplied to function

    rm=sqrt((R-R0)**2+(Z-Z0)**2)

    E_G=F%Ro*h_params%Eo/R

!    E_G=h_params%E*exp(-(rm/h_params%sigma_E)**2/2)
    Z_G=h_params%Zeff*exp(-(rm/h_params%sigma_Z)**2/2)

!    write(6,'("rm: ",E17.10)') rm
    
    A = (2.0_rp*E_G/(Z_G + 1.0_rp))*(g**2 - 1.0_rp)/g
    A = A*h_params%A_fact

    feta = A*EXP(-A*(1.0_rp - COS(deg2rad(eta))))/(1.0_rp - EXP(-2.0_rp*A))     ! MRC
    !	feta = 0.5_rp*A*EXP(A*COS(deg2rad(eta)))/SINH(A)                            ! MRC

    fRE_H_3D = fRE_H_3D*feta
  END FUNCTION fRE_H_3D