!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! D.J. Jeffery, 2020jan06
! An acceleration subroutine for solving the Saha electron density equation.
!      include 'xne_acceleration.f'
! xne = number density of electrons in my jargon.
! x   = reduced number density of electrons in my jargon.
! Since this code fragment works for either, I use "x" inside, but use "xne"
!   as part of the code fragment name.  I've compromised. 
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      include 'xne_quadratic.f'
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
      subroutine xne_acceleration(iacceleration,iter,                    &
     &                 pair,x,x_out,                                     &
     &                 x_min,x_max,x_mid,x_L,x_q,x_accel)
      use precision_set
      use contants
      include 'module_implicit.f'
!      implicit none
      integer, intent(in) :: iacceleration ! 0 for none, 1 for midpoint, 2 for linear, 3 for quadratic.
      integer, intent(in) :: iter          ! The current iteration count.
      real (kind=np), intent(in) :: pair(2,3)  ! The last three (input,output) pairs for the 
!                                              !   the Saha electron density equation evaluation.
      real (kind=np) :: slope ! The slope for the linear acceleration.
      real (kind=np), intent(in) :: x     ! The input for last the Saha electron density equation evaluation.
      real (kind=np), intent(in) :: x_out ! The output for the last Saha electron density equation evaluation.
      real (kind=np), intent(out) :: x_accel ! The accelerated electron density input for the next
!                                            !  Saha electron density equation evaluation. 
      real (kind=np), intent(inout) :: x_min,x_max ! Range limits on the true x 
!                                                  !   whatever that is.
!                                                  !   These can always be found
!                                                  !   since the Saha electron density equation
!                                                  !   is strictly decreasing (I think),
!                                                  !   except perhaps for x input = infinity. 
      real (kind=np), intent(out) :: x_mid ! midpoint or average acceleration value
!                                            !   which may not be any acceleration at all
!                                            !   but is a safe fall back choice.
      real (kind=np), intent(out) :: x_L   ! Linear fit acceleration value. 
      real (kind=np), intent(out) :: x_q   ! Quadratic fit acceleration value
!
!      print*,'In  xne_acceleration'
      x_min=min(x,x_out) 
      x_max=max(x,x_out)
      x_mid=0.5_np*(x_min+x_max)      ! The allowed range mipoint.
!      print*,'Before no-acceleration and midpoint acceleration.'
!      print*,'iacceleration =',iacceleration
!      print*,'x_min,x_max,x_mid,x_accel'
!      print*,x_min,x_max,x_mid,x_accel
      if(iacceleration .eq. 0) then   ! The straight iteration case which
          x_accel=x_out               !   I believe will always converge
!                                     !   but very slowly.
!                                     !   But it's actually only guaranteed
!                                     !   not to diverege, I think.
! From the analytic H case:
!   Results case iacceleration= 0:  Avg iter to prct= 6.053:  Avg iter to mpre=41.632
!   So 6.053 iterations to converge to 1 %, and 41.632 iterations to converge to machine precision.
!
         else
!           print*,'Before the x_accel assignment.'
           x_accel=x_mid  ! Which is probably no acceleration, but it
!                         !   is guaranteed to converge, and so
!                         !   but its a fail-safe if the linear or 
!                         !   quadratic acceleration gives a NaN or is out
!                         !   of the known range.
! From the analytic H case:
!   Results case iacceleration= 1:  Avg iter to prct= 6.842:  Avg iter to mpre=36.158
!   So 6.842 iterations to converge to 1 %, and 36.158 iterations to converge to machine precision.
!   Except that it's guaranteed to coverge, it is probably no better than straight iteration.
!   
      end if
!
!      print*,'Before linear and quadratic acceleration.'
      if(iter .ge. 2) then
          slope=(pair(2,2)-pair(2,1))/(pair(1,2)-pair(1,1))
          x_L=(pair(2,1)-slope*pair(1,1))/(1.0_np-slope)
          if(iacceleration .ge. 2  .and.                                          &
     &       (x_L .ge. x_min .and. x_L .le. x_max)) x_accel=x_L
! From the analytic H case:
!   Results case iacceleration= 2:  Avg iter to prct= 6.579:  Avg iter to mpre=12.737
!   So 6.579 iterations to converge to 1 %, and 12.737 iterations to converge to machine precision.
!   For convergence to 1 %, it's no better than straight iteration and midpoint acceleration. 
!   But it does to converge to machine precision much faster.
!   The x_mid is the fail-safe value if x_L has become a NaN which can happen
!   as machine precision is approached.
!   My compiler does not allow direct tests for NaN's, but they cannot be true for
!   if statements.
!
        else
          x_L=0.0_np
          slope=0.0_np
      end if
      if(iter .ge. 3) then
!          print*,'Before calling xne_quadratic'
          call xne_quadratic(pair,x_q)
          if(iacceleration .ge. 3  .and.                                          &
     &       (x_q .ge. x_min .and. x_q .le. x_max)) x_accel=x_q
! From the analytic H case:
!   Results case iacceleration= 3:  Avg iter to prct= 4.105:  Avg iter to mpre= 8.211
!   So 4.105 iterations to converge to 1 %, and 8.211 iterations to converge to machine precision.
!   This is a significant improvement over the other cases.
!   I think that for large numbers of elements and ionization stages, the extra time overhead
!   of quadratic acceleration is worthwhile.  Also, it's elegant.
!   If the iteration starts with a good initial electron density, covergence to of order 1 % 
!   might be done by the iteration 2 and the quadratic acceleration will never be called.
!   Or it might be called just once allowing convergence on the iteration 4. 
!   I think it's worth the cpu time overhead.  My guess is that quadratic iteration
!   costs about the cpu time of including 2 or 3 elements.  So if one is including 20
!   or more, it's probably worth it.
!   The x_mid or x_L values are fail-safe values if x_q has become a NaN which can happen
!   as machine precision is approached.
!   My compiler does not allow direct tests for NaN's, but they cannot be true for
!   if statements.
!
        else
          x_q=0.0_np
      end if
!
      end subroutine xne_acceleration
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
