!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! D.J. Jeffery, 2020jan06
!
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
! Includes standard modules:  precision,numerical,others eventually
      include 'module_standard.f'
      include 'saha_solver.f'  ! The Saha solver.
!      include 'xne_acceleration.f' ! Electron density solution iteration acceleration.
!        This subroutine is included in saha_solver.f.
      include 'xne_h_analytic.f' ! The analytic H subprograms.  
!                                ! H- is not included in this analytic solution.
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
!
      program test
      use precision_set
      use contants
!      use ieee_exceptions  ! My compiler does not support it.  MRC-221.
!      use ieee_arithmetic  ! My compiler does not support it.  MRC-221.
      include 'test_declarations.f'
!
! Execution lines begin 
!
!  1) Tests of precision and range of types and some variable values.
      print*
      print*,'1) Tests of precision and range of types'
      print*,'     and some variable values.'
      print*,'iprecision_np'
      print*,iprecision_np
      print*,'iminexponent_np,imaxexponent_np'
      print*,iminexponent_np,imaxexponent_np
      print*,'iprecision_ns'
      print*,iprecision_ns
      print*,'iminexponent_ns,imaxexponent_ns'
      print*,iminexponent_ns,imaxexponent_ns
      print*,'phi_0_golden_ratio,phi_1_golden_ratio'
      print*,phi_0_golden_ratio,phi_1_golden_ratio
      print*,'e_rydberg_ev,e_rydberg_ev_reduced'
      print*,e_rydberg_ev,e_rydberg_ev_reduced
!
      print*,'saha_con,saha_con_inverse'
      print*,saha_con,saha_con_inverse
      print*,'   saha_con = constant used in '
      print*,'     Saha ionization stage ratios'
      print*,'    (Mihalas-1978-113:  his value is 2.07e-16 cm**3.)'
      tmp=2.0_np*(pi_sqroot*alpha_fine_structure*xm_e                        &
     &            *clight/h_planck)**3
      print*,'tmp,saha_con_inverse_natural,'
      print*,'tmp,saha_con_inverse_natural,temperature_rydberg'
      print*,tmp,saha_con_inverse_natural,temperature_rydberg
!
!      flag=ieee_is_nan(1.0_np) ! Compiler does not support it.  MRC-221.
!      if(flag .eqv. .true.) stop
!      if(flag .eqv. .false.) stop 'flag stop' ! This works.
      if(flag .eqv. .false.) print*,'Logical variables work.'
!
! 2) Tests for the analytic H ionization or Saha function.
      print*
      print*,'2) Tests for the analytic H ionization'
      print*,'     or Saha function.'
      do410 : do i=-10,10
        xni_reduced=10.0_np**(i)
        xne_e=func_xne_h_analytic_exact(xni_reduced)  ! Exact formula numerically good.
        xne_e2=func_xne_h_analytic_exact2(xni_reduced)! Exact formula numerical bad.
        xne_small=func_xne_h_analytic_series_small(xni_reduced) 
        xne_large=func_xne_h_analytic_series_large(xni_reduced) 
                                                      ! Series result to 3rd order
!                                                     !   for an altnerating series
!                                                     !   where the error is less than
!                                                     !   last term neglected (Arf-294).
        relerr=(xne_e-xne_small)/xne_small
        tmp=(xne_e-xne_large)/xne_large
        tmp2=(xne_e2-xne_small)/xne_small
        write(*,'(i5,8e15.7)') i,xni_reduced,xne_e,relerr,tmp,              &
     &                                       xne_e2,tmp2
      end do do410
!      stop
!
! 3) Tests for convergence: iteration,mipoint,linear,quadratic'
      print*
      print*,'3) Tests for convergence:'
      print*,'     iteration,mipoint,linear,quadratic'
      steps=3.0_np
      steps_factor=10.0_np**(1.0_np/steps)
      xne_initial_min=10.0_np**(-3)
      xne_initial_max=10.0_np**3
      isteps=log10(xne_initial_max/xne_initial_min)                          &
     &       /log10(steps_factor) + 1.0_np
! The math for isteps comes from B=Af**isteps
      xni_reduced=1.0_np  ! Exact hydrogenic Saha equation solution for this is (golden ratio)-1
!                         !   and this our converged solution.
      do412 : do iacceleration_index=0,3     ! Go through 0--3 acceleration cases.
        print*
        write(*,'(a,i2)'),'Acceleration case ',iacceleration_index
        xne_initial=xne_initial_min
        xne_iter_prct=0.0_np    !  Average number of iterations to xne_prct for each acceleration case.
        xne_iter_mpre=0.0_np    !  Average number of iterations to xne_mpre for each acceleration case.
        do420 : do i=1,isteps
          write(*,'(a,i2,a,i2)')                                             &
     &     '  Iteration ',i,                                                 &
     &     ' for iacceleration=',iacceleration_index
          if(iacceleration_index .ge. 3) write(*,fmt=910)
  910     format(' Iter',12x,'xne',8x,'xne_out',8x,'xne_mid',                &
     &              10x,'xne_L',10x,'xne_q',6x,'xne_accel',9x,'relerr')
          xne=xne_initial 
          pair(:,:)=0.0_np   ! Initiate all pairs to zero.
          iter_prct=0
          iter_mpre=0
          iter=0
          do430 : do
            iter=iter+1
            xne_out=func_xne_h_analytic_equation(xne,xni_reduced)  ! The H Saha equation, not its solution.
!                                                                  ! It is used here as iteration equation.
            relerr=(xne_out-xne)/xne
            if(abs(relerr) .le. xne_mpre) then
                xne=xne_out         ! Convergence achieved
                exit do430
            end if
            pair(:,1)=pair(:,2)
            pair(:,2)=pair(:,3)
            pair(1,3)=xne
            pair(2,3)=xne_out
!            print*,'Before calling xne_acceleration'
            call xne_acceleration(iacceleration_index,
     &                  iter,pair,xne,xne_out,                               &
     &                  xne_min,xne_max,xne_mid,xne_L,xne_q,xne_accel)
            if(iter_prct .eq. 0  .and.                                       &
     &       abs(xne_accel-phi_1_golden_ratio)/phi_1_golden_ratio            &
     &       .le. xne_prct) iter_prct=iter
            if(iter_mpre .eq. 0  .and.                                       &
     &       abs(xne_accel-phi_1_golden_ratio)/phi_1_golden_ratio            &
     &       .le. xne_mpre) iter_mpre=iter
            if(iacceleration_index .ge. 3)                                   &
     &        write(*,'(1p,i5,8e15.7)') iter,xne,xne_out,xne_mid,            &
     &                                  xne_L,xne_q,xne_accel,relerr
            xne=xne_accel  ! For the next evaluation of the electron density equation.
            if(iter .eq. 50) then
              print*,'iacceleration case ',iacceleration_index
              stop 'Exceeded allowed iterations for test cases.'
            end if
          end do do430
          write(*,920) iacceleration_index,iter,iter_prct,iter_mpre
  920     format(2x,'iacceleration=',i2,' iteration count=',i3,               &
     &           ' iter_prct=',i3,' iter_mpre=',i3)
          xne_iter_prct=xne_iter_prct+real(iter_prct,kind=np)
          xne_iter_mpre=xne_iter_mpre+real(iter_mpre,kind=np)
          xne_initial=xne_initial*steps_factor
        end do do420
        write(*,930) iacceleration_index,                                      &
     &               xne_iter_prct/real(isteps,kind=np),                     &
     &               xne_iter_mpre/real(isteps,kind=np)
  930   format('  Results case iacceleration=',i2,                           &
     &         ':  Avg iter to prct=',f6.3,                                  &
     &         ':  Avg iter to mpre=',f6.3)
      end do do412
      stop
!
!  4) Use saha_solver just to read-in and print out the partition functions,
!       both from partition_function_read.f and saha_solver_f to check agreement.
      print*,
      print*,'4) Use saha_solver just to read-in'
      print*,'     and print out the partition functions.' 
      idiagnostic=3
      include 'saha_solver_call.f'  ! The Saha solver.
!
!  5) Use saha_solver for pure H (without H-).
      print*
      print*,'5) Use saha_solver for pure H (without H-).'
      tem=7.e+3_np
      kkel=1
      kel(1)=1 
      xel(1)=1.0_np
      idiagnostic=2
      ih_minus_not=1  ! Means do not use H-.
      xne_converged=xne_mpre
      ixne_initial=1
      idefault=0
      include 'saha_solver_call.f'  ! The Saha solver.
! xne_h,xnum_total,phi(0,1)
!  4.9524433E+10  6.0221408E+10  4.3613598E-12
! The Saha solution for:
! iacceleration=    3 rho=  1.0000000E-13 T=  7.0000000E+03 xne_initial=  4.9524433E+10
! kel(:)           1
! xel(:)   1.0000000000000000000
! iter            xne        xne_out      xne_accel         relerr     Ionization
!    1  4.9524433E+10  4.9524433E+10  0.0000000E+00 -3.0088504E-19  8.2237256E-01
!
      xne_solution=xne  ! Use the output xne of the last test and the huge
!                               !   range for analytic hydrogen test above.
      write(*,"('The solution xne=',1p,e15.7)") xne_solution
      xne_tmp=xne_solution*xne_initial_min  ! Use the output xne of the last test and the huge
!                                           !   range for analytic hydrogen test above.
      do510 : do i=1,isteps
        idiagnostic=1
        ih_minus_not=1
        ixne_initial=2
        print*
        xne=xne_tmp
        write(*,'(a,i2,a,i2,1p,2(a,e15.7))') 'Calculation=',i,              &
     &                         ' with iacceleration=',iacceleration,        &
     &                         ' and solution xne=',xne_solution,           &
     &                         ' for initial xne=',xne
        idefault=0
        include 'saha_solver_call.f'  ! The Saha solver.
        xne_tmp=xne_tmp*steps_factor
      end do do510
!
! 6) Solve for Carbon ionization fractions.
!    http://cococubed.asu.edu/code_pages/eos_ionize.shtml for calculated results.
!
      print*
      print*,'6) Solve for Carbon ionization fractions.'
      kkel=1
      kel(1)=6   ! Pure carbon.
      xel(1)=1.0_np
      steps=3.0_np
      steps_factor=10.0_np**(1.0_np/steps)
      rho=1.0e-4_np
      tem_min=1.0e+3_np
      tem_max=1.0e+7_np
      isteps=log10(tem_max/tem_min)                                           &
     &       /log10(steps_factor) + 1.0_np
      tem=tem_min
!
      do610 : do i=1,isteps
        idefault=0
        idiagnostic=2
        ixne_initial=2
        include 'saha_solver_call.f'  ! The Saha solver.
        write(*,942) i,tem,fraction(0:6,1)
        tem=tem*steps_factor
      end do do610 
  942 format(i5,1p,8e15.7) 
!
      end program test
!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
