!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! D.J. Jeffery, 2020jan06
! Functions used by saha_solver.f
!      include 'saha_solver_functions.f'  ! Includes some saha_solver.f functions.
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      function func_phi(u,u1,xne,xion_en,tem_use)    ! xne*phi-tilde function (Mihalas-1978-113)
      use precision_set
      use contants
      include 'module_implicit.f'
      include 'function_type_declaration.f'
      real (kind=np), intent(in) :: xne      ! electron density.
      real (kind=np), intent(in) :: tem_use  ! Temperature. 
      real (kind=np), intent(in) :: u        ! j partition function. 
      real (kind=np), intent(in) :: u1       ! j+1 partition.
      real (kind=np), intent(in) :: xion_en  ! Ionization energy from j to j+1.
!
      func_phi=(u/u1)*func_phi_r(xne,xion_en,tem_use) 
      end function func_phi 
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      function func_phi_r(xne,xion_en,tem_use)    ! Same as above without partition function. 
      use precision_set
      use contants
      include 'module_implicit.f'
      include 'function_type_declaration.f'
      real (kind=np), intent(in) :: xne      ! electron density.
      real (kind=np), intent(in) :: tem_use  ! Temperature.
      real (kind=np), intent(in) :: xion_en  ! Ionization energy from j to j+1.
!
      func_phi_r=xne*( (saha_con*tem_use**(-1.5_np))                          &
     &                *exp(xion_en/(xk_boltzmann_ev*tem_use))  )
      end function func_phi_r
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      function func_phi_r_inverse(xne,xion_en,tem_use)    ! Inverse of the above. 
      use precision_set
      use contants
      include 'module_implicit.f'
      include 'function_type_declaration.f'
      real (kind=np), intent(in) :: xne      ! electron density.
      real (kind=np), intent(in) :: tem_use  ! Temperature.
      real (kind=np), intent(in) :: xion_en  ! Ionization energy from j to j+1.
!
      func_phi_r_inverse=(1.0_np/xne)                                          &
     &            *( (saha_con_inverse*tem_use**(1.5_np) )                     &
     &               *exp(-xion_en/(xk_boltzmann_ev*tem_use))  )
      end function func_phi_r_inverse
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12