FUNCTION fRE_H_3D(params,F,eta,g,R,Z,R0,Z0,EPHI,rho_ind)
TYPE(KORC_PARAMS), INTENT(IN) :: params
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), INTENT(IN),optional :: EPHI
INTEGER, INTENT(IN),optional :: rho_ind
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
REAL(rp) :: CLog0,CLogee,E_CH,k=5._rp,VTe
INTEGER :: index
index = MINLOC(ABS(h_params%g - g),1)
!write(6,*) index
! 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(present(rho_ind)) then
if (D.GT.0) then
f0 = h_params%fRE_E_2D(rho_ind,index-1)
g0 = h_params%g(index-1)
f1 = h_params%fRE_E_2D(rho_ind,index)
g1 = h_params%g(index)
else
f0 = h_params%fRE_E_2D(rho_ind,index)
g0 = h_params%g(index)
f1 = h_params%fRE_E_2D(rho_ind,index+1)
g1 = h_params%g(index+1)
end if
else
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
endif
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)
if(present(EPHI)) then
CLog0=14.9_rp - LOG(1E-20_rp*n_ne)/2._rp + &
LOG(1E-3_rp*Teo)
VTe=sqrt(2._rp*Teo*C_E/C_ME)
CLogee=CLog0+log(1+(2*(g-1)/(VTe/C_C)**2)**(k/2._rp))/k
E_CH=n_ne*C_E**3*CLogee/(4*C_PI*C_E0**2*C_ME*C_C**2)
!write(output_unit_write,*) 'gamma',g
!write(output_unit_write,*) 'ne',n_ne, &
! 'Te',Teo,'VTe',VTe
!write(output_unit_write,*) 'EPHI',EPHI,'CLog',CLogee,'E_CH',E_CH
!flush(output_unit_write)
E_G=abs(EPHI)/E_CH
!write(output_unit_write,*) 'E_G',E_G
!flush(output_unit_write)
else
E_G=F%Ro*h_params%Eo/R
endif
! 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(output_unit_write,'("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
!write(output_unit_write,*) 'A',A
!flush(output_unit_write)
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
!if (fRE_H_3D.eq.0._rp) then
! write(6,*) 'f_RE_H',f0 + m*(g - g0)
! write(6,*) 'feta',feta
! write(6,*) 'A',A
! write(6,*) 'g',g
! write(6,*) 'E_G',E_G
! write(6,*) 'E_CH',E_CH
! write(6,*) 'Z_G',Z_G
!end if
END FUNCTION fRE_H_3D