rate_coeff_neutrals_m.f90 Source File

module rate_coeff_neutrals_m
!! Module for calculating fits for rate coefficients k_iz, k_rec
!!   and cooling rate coefficients  w_iz, p_rec
    use precision_grillix_m, only : GP, GP_NAN, GP_EPS
    use error_handling_grillix_m, only: handle_error
    use status_codes_grillix_m, only : GRILLIX_ERR_NAMELIST
    use variable_m, only: variable_t
    use screen_io_m, only :  get_stdout
    
    implicit none
    
    type, public :: rate_coeff_neutrals_t
        !! Class for rate coefficients (ionization & recombination)
        real(GP), private, dimension(0:8,0:8) :: fit
        !! Polynomial coefficients for calculating rate coefficient fit
        real(GP), private :: ne_ref
        !! Reference for electron density
        real(GP), private :: te_ref
        !! Reference for electron temperature
        real(GP), private :: k_ref
        !! Reference for rate coefficient

    contains
        procedure, public :: create => create_rate_coeff
        procedure, public :: eval => eval_rate_coeff
    end type

contains

    subroutine create_rate_coeff(self, file_name, ne_ref, te_ref, k_ref)
        !! Create rate coefficient object with fit coefficients taken from data file
        class(rate_coeff_neutrals_t), intent(inout) :: self
        !! Instance of class
        character(len=*), intent(in) :: file_name
        !! Name of atomic data file to read coefficients from
        real(GP), intent(in) :: ne_ref
        !! Reference for electron density
        real(GP), intent(in) :: te_ref
        !! Reference for electron temperature
        real(GP), intent(in) :: k_ref
        !! Reference for rate coefficient
        
        integer :: l, io_error
        character(len=256) :: io_errmsg

        open(unit = 40, file = file_name , status='old', action='read', iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) then
            call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        endif
        do l = 0, 8
            read(40,*, iostat=io_error, iomsg=io_errmsg) self%fit(l,0:8)
            if (io_error /= 0) then
                call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
            endif
        end do
        close(40)
        self%ne_ref = ne_ref
        self%te_ref = te_ref
        self%k_ref  = k_ref

    end subroutine

    real(GP) function eval_rate_coeff(self, te_in, ne_in)
        !! Evaluate the normalized rate coefficient from a double polynomial fit
        class(rate_coeff_neutrals_t), intent(in) :: self
        !! Instance of class
        real(GP), intent(in) :: te_in
        !! Electron temperature
        real(GP), intent(in) :: ne_in
        !! Plasma density

        real(GP) :: te_fit, ne_fit

        ! Convert inputs to log values, as required for fit (see AMJUEL manual)
        te_fit = log(te_in * self%te_ref)
        ne_fit = log(ne_in * self%ne_ref * 1.0E-8_GP ) ! apply AMJUEL electron density scaling
        eval_rate_coeff = exp( poly8x8_eval(self%fit, te_fit, ne_fit) ) / self%k_ref

    end function

    real(GP) function poly8x8_eval( pc, x, y )
        !! Returns double polynomial of degree 8x8, with coefficients pc(x,y), evaluated at x and y.
        real(GP), dimension(0:8,0:8), intent(in) :: pc
        !! polynomial coefficients
        real(GP), intent(in) :: x, y
        !! input  variables

        real(GP), dimension(0:8) :: pcx
        integer :: m

        do m = 0, 8
            ! collapse to single polynomial
            pcx(m) = poly8_eval( pc(m,0:8), y )
        end do
        ! calculate remaining polynomial
        poly8x8_eval = poly8_eval( pcx, x )

    end function

    real(GP) function poly8_eval( pc, x )
        !! Returns polynomial of degree 8, with coefficients pc(0:8), evaluated at x.
        real(GP), dimension(0:8), intent(in) :: pc
        !! polynomial coefficients
        real(GP), intent(in) :: x
        !! input variable

        integer :: n

        poly8_eval = 0.0_GP
        do n = 8, 1, -1
            poly8_eval = ( poly8_eval + pc(n) ) * x
        end do
        poly8_eval = poly8_eval + pc(0)

    end function

end module