calculate_k_matrix Subroutine

public subroutine calculate_k_matrix(number_of_channels, log_der_matrix, number_of_open_channels, channel_indices, channel_l_values, r_, k_matrix)

calculates the K-matrix from log-derivative matrix using Eq. (4) in "Solution of coupled equations"

Arguments

Type IntentOptional Attributes Name
integer(kind=int32), intent(in) :: number_of_channels

y-matrix is of number_of_channels x number_of_channels size
real(kind=dp), intent(in) :: log_der_matrix(number_of_channels,number_of_channels)

asymptotic log-derivative matrix
integer(kind=int32), intent(in) :: number_of_open_channels

number of open channels
integer(kind=int32), intent(in) :: channel_indices(number_of_channels)

holds the indices pointing to the basis arrays
integer(kind=int32), intent(in) :: channel_l_values(number_of_channels)

holds all values of l
real(kind=dp), intent(in) :: r_

r_max
real(kind=dp), intent(inout) :: k_matrix(number_of_open_channels,number_of_open_channels)

K-matrix

Calls

proc~~calculate_k_matrix~~CallsGraph proc~calculate_k_matrix calculate_k_matrix proc~is_open is_open proc~calculate_k_matrix->proc~is_open dgesv dgesv proc~calculate_k_matrix->dgesv dgemm dgemm proc~calculate_k_matrix->dgemm proc~riccati_bessel_j riccati_bessel_j proc~calculate_k_matrix->proc~riccati_bessel_j proc~modified_bessel_k_ratio modified_bessel_k_ratio proc~calculate_k_matrix->proc~modified_bessel_k_ratio proc~riccati_bessel_y riccati_bessel_y proc~calculate_k_matrix->proc~riccati_bessel_y proc~wavevector_squared_from_energy wavevector_squared_from_energy proc~calculate_k_matrix->proc~wavevector_squared_from_energy proc~total_energy total_energy proc~is_open->proc~total_energy proc~rctj rctj proc~riccati_bessel_j->proc~rctj proc~write_warning write_warning proc~riccati_bessel_j->proc~write_warning proc~integer_to_character integer_to_character proc~riccati_bessel_j->proc~integer_to_character proc~modified_bessel_temme_algorithm modified_bessel_temme_algorithm proc~modified_bessel_k_ratio->proc~modified_bessel_temme_algorithm proc~riccati_bessel_y->proc~write_warning proc~rcty rcty proc~riccati_bessel_y->proc~rcty proc~riccati_bessel_y->proc~integer_to_character proc~wavevector_squared_from_energy->proc~total_energy proc~write_error write_error proc~wavevector_squared_from_energy->proc~write_error proc~modified_bessel_temme_algorithm->proc~write_error proc~rgamma rgamma proc~modified_bessel_temme_algorithm->proc~rgamma proc~msta2 msta2 proc~rctj->proc~msta2 proc~msta1 msta1 proc~rctj->proc~msta1 proc~write_message write_message proc~write_warning->proc~write_message proc~write_error->proc~write_message proc~envj envj proc~msta2->proc~envj proc~msta1->proc~envj

Called by

proc~~calculate_k_matrix~~CalledByGraph proc~calculate_k_matrix calculate_k_matrix program~scattering SCATTERING program~scattering->proc~calculate_k_matrix

Contents

Variables

Source Code

calculate_k_matrix

Variables

Type Visibility Attributes Name Initial
integer(kind=int32), private :: channel_index_
integer(kind=int32), private :: closed_channel_index_
integer(kind=int32), private :: closed_channels_indices(number_of_channels-number_of_open_channels)
real(kind=dp), private :: diag_j_matrix(number_of_channels,number_of_open_channels)
real(kind=dp), private :: diag_jp_matrix(number_of_channels,number_of_open_channels)
real(kind=dp), private :: diag_n_matrix(number_of_channels,number_of_channels)
real(kind=dp), private :: diag_np_matrix(number_of_channels,number_of_channels)
real(kind=dp), private :: j_element_
real(kind=dp), private :: jp_element_
integer(kind=int32), private :: l_
real(kind=dp), private :: n_element_
real(kind=dp), private :: np_element_
integer(kind=int32), private :: open_channel_index_
integer(kind=int32), private :: open_channels_indices(number_of_open_channels)
real(kind=dp), private :: ratio
integer(kind=int32), private :: status_
real(kind=dp), private :: wavevector
real(kind=dp), private :: x

Source Code

      subroutine calculate_k_matrix(number_of_channels, log_der_matrix,        &
         number_of_open_channels, channel_indices, channel_l_values, r_,       &
         k_matrix)
         !! calculates the K-matrix from log-derivative matrix using Eq. (4) in
         !! "Solution of coupled equations" 
         !---------------------------------------------------------------------!
         integer(int32), intent(in) :: number_of_channels
            !! y-matrix is of number_of_channels x number_of_channels size
         real(dp), intent(in) :: log_der_matrix(number_of_channels,            &
            number_of_channels)
            !! asymptotic log-derivative matrix
         integer(int32), intent(in) :: number_of_open_channels
            !! number of open channels
         integer(int32), intent(in) :: channel_indices(number_of_channels)
            !! holds the indices pointing to the basis arrays
         integer(int32), intent(in) :: channel_l_values(number_of_channels)
            !! holds all values of l
         real(dp), intent(in) :: r_
            !! r_max
         real(dp), intent(inout) :: k_matrix(number_of_open_channels,          &
            number_of_open_channels)
            !! K-matrix
         !---------------------------------------------------------------------!
         integer(int32) :: open_channel_index_, closed_channel_index_,         &
            channel_index_, status_, l_
         real(dp) :: wavevector, x, j_element_, jp_element_, n_element_,       &
            np_element_, ratio
         integer(int32) :: open_channels_indices(number_of_open_channels)
         integer(int32) :: closed_channels_indices(                            &
            number_of_channels-number_of_open_channels)
         real(dp) ::  diag_n_matrix(number_of_channels, number_of_channels),   &
            diag_np_matrix(number_of_channels, number_of_channels),            &
            diag_j_matrix(number_of_channels, number_of_open_channels),        &
            diag_jp_matrix(number_of_channels,number_of_open_channels)
         !---------------------------------------------------------------------!
         ! diag_j_matrix   -  diagonal J matrix (Eqs. 5, 7)
         ! diag_jp_matrix  -  diagonal J`matrix (derivative of J)
         ! diag_n_matrix   -  diagonal N matrix (Eqs. 6, 8)
         ! diag_np_matrix  -  diagonal N`matrix (derivative of N)
         !---------------------------------------------------------------------!
         diag_j_matrix  = 0
         diag_jp_matrix = 0
         diag_n_matrix  = 0
         diag_np_matrix = 0
         !---------------------------------------------------------------------!
         open_channel_index_   = 0
         closed_channel_index_ = 0
         !---------------------------------------------------------------------!
         ! save indices to open and closed channels
         ! this is because channels might not be sorted eneregetically
         !---------------------------------------------------------------------!
         do channel_index_ = 1, number_of_channels
            if (is_open(internal_energies(channel_indices(channel_index_)))) then
               open_channel_index_ = open_channel_index_+1
               open_channels_indices(open_channel_index_) = channel_index_
            else
               closed_channel_index_ = closed_channel_index_+1
               closed_channels_indices(closed_channel_index_) = channel_index_
            endif
         enddo
         !---------------------------------------------------------------------!
         ! Prepare J, J', N and N' matrices (Eqs. 5-6)
         ! open channels:
         !---------------------------------------------------------------------!
         do open_channel_index_ = 1, number_of_open_channels
            wavevector = sqrt( wavevector_squared_from_energy(                 &
               internal_energies(channel_indices(open_channels_indices(open_channel_index_)))))
            x  = wavevector*r_
            l_ = channel_l_values(open_channels_indices(open_channel_index_))
            call riccati_bessel_j(                                             &
               l_, x, j_element_, jp_element_)
            diag_j_matrix(open_channel_index_,open_channel_index_)             &
               = (wavevector)**(-0.5_dp)*j_element_
            diag_jp_matrix(open_channel_index_,open_channel_index_)            &
               = (wavevector)**(0.5_dp)*jp_element_

            call riccati_bessel_y(l_, x, n_element_, np_element_)
            diag_n_matrix(open_channel_index_,open_channel_index_)             &
               = (wavevector)**(-0.5_dp)*n_element_
            diag_np_matrix(open_channel_index_,open_channel_index_)            &
               = (wavevector)**(0.5_dp)*np_element_
         enddo
         !---------------------------------------------------------------------!
         ! Prepare J, J', N and N' matrices (Eqs. 7-8)
         ! closed channels:
         !---------------------------------------------------------------------!
         do closed_channel_index_ = 1,number_of_channels-number_of_open_channels  
            wavevector = sqrt( abs( wavevector_squared_from_energy(            &
               internal_energies(channel_indices(closed_channels_indices(closed_channel_index_))))))
            x  = wavevector*r_
            l_ = channel_l_values(closed_channels_indices(closed_channel_index_))
            call modified_bessel_k_ratio(l_,x,ratio)
            !------------------------------------------------------------------!
            ! substitution for closed channels, (Eqs. 10 - 11)             
            !------------------------------------------------------------------!
            diag_n_matrix(number_of_open_channels+closed_channel_index_,       &
               number_of_open_channels+closed_channel_index_) = 1.0_dp
            diag_np_matrix(number_of_open_channels+closed_channel_index_,      &
               number_of_open_channels+closed_channel_index_) = wavevector*ratio
         enddo
         !---------------------------------------------------------------------!
         call DGEMM('N','N',number_of_channels,number_of_channels,             &
            number_of_channels,1.0_dp,log_der_matrix,number_of_channels,       &
            diag_n_matrix,number_of_channels,-1.0_dp,diag_np_matrix,           &
            number_of_channels)
         call DGEMM('N','N',number_of_channels,number_of_open_channels,        &
            number_of_channels,-1.0_dp,log_der_matrix,number_of_channels,      &
            diag_j_matrix,number_of_channels,1.0_dp,diag_jp_matrix,            &
            number_of_channels)
         !---------------------------------------------------------------------!
         call DGESV(number_of_channels,number_of_open_channels,diag_np_matrix, &
            number_of_channels,diag_j_matrix,diag_jp_matrix,number_of_channels,status_)
         !---------------------------------------------------------------------!
         k_matrix = diag_jp_matrix(1:number_of_open_channels, 1:number_of_open_channels)
         !---------------------------------------------------------------------!
      end subroutine calculate_k_matrix