!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
      program black 
!      implicit real(kind=kind(1.d0)) (a-h,o-z)
!      implicit real(kind=selected_real_kind(16)) (a-h,o-z)  ! Works.
      implicit real(kind=selected_real_kind(32)) (a-h,o-z)  ! Works.
!      implicit real(kind=SELECTED_REAL_KIND(64)) (a-h,o-z)  !  Don't Work
! See https://software.intel.com/en-us/forums/intel-fortran-compiler-for-linux-and-mac-os-x/topic/269257
!      implicit real(kind=16 (a-h,o-z)  ! Doesn't work.
!      implicit selected_real_kind(2*precision(1.0_dp)) (a-h,o-z) ! Doesn't work.
!      implicit real(kind=32) (a-h,o-z) ! Does not work.
!      implicit real_kind_precision24 (a-h,o-z) ! Does not work.
!      implicit real(kind=kind(dp)) (a-h,o-z) ! Accepts but what does it do?
!      implicit real(kind=kind(qp)) (a-h,o-z) ! Accepts but what does it do?
!      integer, parameter :: nprecision=kind(1.d0)
!      integer, parameter :: np = selected_real_kind(16)   ! Works.
      integer, parameter :: np = selected_real_kind(32)   ! Works.
!      integer, parameter :: np = selected_real_kind(64)   ! Don't Works.
!      integer, parameter :: nprecision=kind(32)
!      integer, parameter :: nprecision=kind(dp)
!      integer, parameter :: nprecision=kind(qp)
      integer, parameter :: nfa=16
      integer, parameter :: nprime=8
      integer, parameter :: npz=20
      real (kind=np) :: fa(0:nfa)
      integer :: iiprime(nprime)
      real (kind=np) :: pz(npz)
      data pz/1.e-6_np,1.e-5_np,                                                 &
     &        0.0001_np, 0.001_np, 0.01_np, 0.03_np, 1.1_np,                     &
     &        1.15_np,                                                           &
     &        1.2_np,                                                            & 
     &        1.3_np,   1.4_np,  1.5_np,  1.6_np,                                &
     &        1.7_np,    1.8_np,   1.9_np,  2.0_np,   3.0_np,                    &
     &        4.0_np,    5.0_np/
      data iiprime/1,2,3,5, 7,11,13,17/
! From https://en.wikipedia.org/wiki/Factorial
      data fa/1.0_np, 1.0_np, 2.0_np, 6.0_np, 24.0_np,                           &
     &        120.0_np, 720.0_np, 5040.0_np, 40320.0_np,                         &
     &        362880.0_np, 3628800.0_np, 39916800.0_np, 479001600.0_np,          &
     &        6227020800_np, 87178291200_np, 1307674368000_np,                   &
     &           20922789888000_np/
!
      print*
      print*,'Blackbody radiation distribution and integrals'
      print*,'Related to Debye function integrals'
      print*,'which are effectively, Riemann Zeta functions'
      print*,'See Arfken p. 332'
      print*,'Max of integrands exponent p=1.01 to 5'
      write(*,"(1x,'iteration',13x,'iterate',6x,'Newton-Raphson')")
!             https://en.wiktionary.org/wiki/iterate the noun
      con2=1.5936242600400403_np ! Numerical solution for p=2
!
      iexact=0
      if(iexact .eq. 1) then
          pde=1.0e-6_np
          pe=1.0_np+pde
              ! 23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
          xe1=1.99999933333377777745181382492999647E-0006_np  ! what np=32 gives
!             1.99999933333377777745181382501483924E-0006     ! After 1 000 000 iterations converged
!             1.99999933333377777745181382499397862E-0006     ! After converged
!             1.99999933333377777745183802478120206E-06       ! one less that total
!             1.99999933333377777745183802479528595E-0006     ! converged to 29 digits  
!          xe1=1.9933774e-2_np  ! Yes, they do iterate
          xe2=xe1
          xe3=xe1
          ii=1 000 000
          do i=1,ii
            xe1=pe*(1.0_np-exp(-xe1))
            xe2=pe*(exp(xe2)-1.0_np)/exp(xe2)
            sum=0.0_np
            do j=nfa,1,-1
              sum=sum+(-xe3)**j/fa(j)
            end do
            xe3=-pe*sum
            if(i .gt. ii-11) write(*,910) i,xe1,xe2,xe3
          end do 
  910     format(i7,3e45.36) 
          print*,'pe,xe1,xe2,xe3'
          print*,pe,xe1,xe2,xe3
          stop
      end if
!
      print*,'The expansion terms'
         a1=2.0_np           ! determined analytically 1st
         a2=-2.0_np/3.0_np   ! determined analytically 2nd
         a3=4.0_np/9.0_np    ! determined numerically  3rd pd=1.e-5 to relerr 1e-5
         a4=-1.0_np/3.0_np   ! Approximately right, fiducial value
!         a4=-0.325963952847781078141194663752917956_np ! Very small delz fit but not right
         a5=1.0_np/3.0_np    ! Fiducial value at z=1.1
!         a5=0.0_np           ! a5 can't be fit unless, a4 were known exactly, without more cleverness.
!         a5=1.0_np/4.0_np     ! Consistent with a4=-44/135 ?, but it doesn't stablize xcoef6
!         a5=1.0_np           ! determined numerically  5th  this looks right at delz=.01 to 0.5 % 
!         a5=1.0_np/3.0_np    ! determined numerically  5th  this looks right at delz=.1 
         a6=-2.0_np/9.0_np   ! Fiducial value at z=1.1.   
!         a6=0.0_np 
! p = z xlow-x2 
! 1.1  2.43943495082172799567213935641107047E-0009
! 1.2  1.00027505387241362096291602996416415E-0005
! 1.3  8.70623203195004906480843155165127304E-0005
! 1.4  3.31905192729496031072133295792039345E-0004
! 1.5  7.82534201282920938385984169344795317E-0004
! 1.6  1.18189225328027400815451880144153590E-0003
! 1.7  6.07272745384688824701656951475041720E-0004  
! 1.8  -3.04593994617942485727758941976983294E-0003
! 1.9  -1.37308194737776213527942068124790343E-0002
! 2.0  -3.80687044844845367674863203196046635E-0002 probably on the brink of divergence 
! 3.0  -6.82143937212207889340319133029448507  clearly diverging
         a7=2.d0/81.d0   ! Fiducial value at z=1.1 .  I don't trust.  Don't use.
!         a7=0.0_np         ! Can't determine it      
         a8=0.0_np
         a9=0.0_np
         b1=1.0_np   ! Determined by Bell polynomials and confirmed. See p. 19
         b2=1.0_np
         b3=1.5_np
         b4=8.0_np/3.0_np
         b5=125.0_np/24.0_np
         b6=54.0_np/5.0_np  ! Not confirmed
      isatz=2     ! Which interpolation formula to use.  ! https://en.wikipedia.org/wiki/Ansatz 
      icorrect=0  ! 0 for none, 1 for interpolation, 2 for Newton-Raphson
      itest=0     ! Test for the low series coefficients
      if(isatz .eq. 1) then            ! Crude interpolation formula 
          acon=2.0_np
        else if(isatz .eq. 2) then            ! Gem formula 
          acon=2.0_np/(1-exp(-1.0_np))
          bcon=exp(2.0_np)*(1-exp(-1.0_np))/3.0_np
          bcon2=(exp(-1.0_np)                                                &
     &           -bcon*exp(-2.0_np))/(1.0_np-exp(-1.0_np))
          bcon3=bcon2                                                        &
     &     +(-0.5_np*exp(-1.0_np)+bcon*2.0_np*exp(-2.0_np))                  &
     &            /(1.0_np-exp(-1.0_np))
          acon2=(1.0_np/6.0_np)*acon**2
          acon3=0.5_np*acon*(-1.0_np-bcon2)
          ccon=-exp(3.0_np)*(1.0_np-exp(-1.0_np))                            &
     &          *(a3-bcon3-acon2-bcon2-acon3)
          print*,'For gem formula'
          print*,'acon,bcon,ccon,b3'
          print*,acon,bcon,ccon,b3
          x01=0.1937475579949906_np
          pd=0.1_np
          p=1.0_np+pd
          ccon2=-(  (x01/(p*(1.0_np-exp(-acon*pd))))-1.0_np+exp(-p)                &
     &              +bcon*pd*exp(-2.0_np*p) )                                   &
     &          /(pd**2*exp(-3.0_np*p))
          print*,'ccon,ccon2'
          print*,ccon,ccon2
!          stop
        else if(isatz .eq. 3) then    ! Series interpolation formula 
          print*,'isatz=',isatz
          xtrans=1.85_np               ! seems good place to sharp transition
        else if(isatz .eq. 4) then    ! Series-smooth interpolation formula
          pcon=4.0_np                          ! smooth transition power
          efold=.9_np/log(2.0_np)**(1.0_np/pcon)  ! smooth transition point
        else if(isatz .eq. 5) then     ! The newer brilliant flop formula
          print*,'isatz = ',isatz
          acon=2.0_np/(1.0_np-exp(-1.0_np))
          bcon=1.5569247568238684_np
          ccon=1.5442655706398485_np
        else if(isatz .eq. 9) then     ! The new brilliant flop formula
          print*,'isatz = ',isatz
        else if(isatz .eq. 10) then    ! The brilliant flop formula
          print*,'isatz = ',isatz
      end if
      do i=1,npz              ! Run through the chosen p=z's.
        print*
        if(pz(i) .le. 0.11_np) then
            pd=pz(i)
            p=1.0_np+pd
          else 
            p=pz(i)
            pd=p-1.0_np
        end if
        write(*,'(a,i4,a,e10.3,a,e10.3)')                                &
     &      'Iteration cycle ',i,' for z=p=',p,' delz=pd=',pd
        yexp=exp(-p)
        y=p*yexp
        xlow2=pd*(a1+pd*(a2))
        xlow3=pd*(a1+pd*(a2+pd*(a3)))
        xlow4=pd*(a1+pd*(a2+pd*(a3+pd*(a4))))
        xlow5=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5)))))
        xlow6=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6))) )) 
        xlow8=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6                   &
     &              +pd*(a7+pd*a8))) )))     ! High precision form.  
        xlow=xlow6
!        xhigh=p-(y+y**2+b3*y**3+b4*y**4+b5*y**5+b6*y**6)
        xhigh=p-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) ))   ! High precision form
!
        if(isatz .eq. 1) then                              ! Crude interpolation formula
            x0=p*(1.0_np-exp(-acon*pd))
! prec max = 4, double prec max =7, relerr max= 9.0014819399899054E-002  at 1.7  good
! prec max = 3, double prec max =6, relerr max= 4.5495036393012151E-003 at 1.6   good with NR correction
          else if(isatz .eq. 2) then                       ! Gem formula
            x0=p*(1.0_np-exp(-p)-bcon*pd*exp(-2.0_np*p)                                       &
! prec max = 3, double prec max =6, relerr max= 6.5047149472154484E-003 at 1.7  good
!     &                -b3*pd**2*exp(-3.0_np*p)                                              &
! prec max = 3, double prec max =6, relerr max= 7.3877224167747660E-004 at 1.8  much better
! prec max = 2, double prec max =5, relerr max= 3.4485709062780802E-007 at 1.7  excellent with NR correction
     &                -ccon*pd**2*exp(-3.0_np*p)                                            &
! prec max = 3, double prec max =6, relerr max= 5.7797367485568822E-004  at 1.8  slightly significantly better 
!     &                -0.1_np*pd**2*exp(-3.0_np*p)                                           &
! prec max = 3, double prec max =6, relerr max= -4.8571899523800627E-004 at 3.0 marginal improve with fit
!     &                -b4*pd**3*exp(-4.0_np*p)                                              &
! prec max = 3, double prec max =6, relerr max= 6.8422439795701086E-004 at 1.4  marginally better
!                                   relerr    = 5.8364429522015239E-004 at 1.7  probably negligible.
!     &                -b4*pd**2*p*exp(-4.0_np*p)                                            &
! prec max = 3, double prec max =6, relerr max= -2.9347868449807234E-003 at 1.5  definitely worse 
!     &                -b4*pd**3*exp(-4.0_np*p)-b5*pd**4*exp(-5.0_np*p)                       &
! prec max = 3, double prec max =6, relerr max= -9.3720229008774933E-004 at 1.6  marginally worse
! 
     &                                           )                                     &
     &         *(1.0_np-exp(-acon*pd))
           else if(isatz .eq. 3) then                        ! Series interpolation formula
              if(p .le. xtrans) then  ! Numerically clear best just to use a sharp transition for accuracy.
                  x0=xlow
                else
                  x0=xhigh
              end if
! prec max = 3, double prec max =6, relerr max=  6.9558650904086279E-003 at 1.9 good
!   So the maximum is worse by a factor 10 than the magic formula.
!   Maybe the magic formula is just semi-accidentally a good fit.
          else if(isatz .eq. 4) then                       ! Series smooth-transition interpolation formula
              w=exp(-(pd/efold)**pcon)
              x0=xlow*w+xhigh*(1.0_np-w)
          else if(isatz .eq. 5) then                         ! Newer briiliant flop formula. isatz = 5
            y=pd*exp(-p)
            wlow=exp(-acon*pd)
            w=1.0_np-exp(-acon*pd)
!            x0=1.0_np-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) ))/pd   ! High precision form
            x0=1.0_np-y*(b1+y*(bcon+y*(ccon+y*(b4+y*(b5+y*b6))) ))/pd   ! High precision form hybrid 
            x0=p*x0*w 
! prec max = 4, double prec max =9, relerr max= 1.8141849762850569E-002 at 1.4 not competitive 
! prec max = 4, double prec max =6, relerr max= -1.1862074441575882E-003 at 1.6  hybrid, worse 
!            x0=x0*w+xlow*wlow
! prec max = 4, double prec max =6, relerr max= 1.4647559416957010E-002 at 1.6 marginal about the same
          else if(isatz .eq. 9) then                         ! New briiliant flop formula.   isatz = 9 
             pd2=2.0_np*pd                                   ! Never bad, but up to 10 % error.
             y2=pd2*exp(-pd2)                                ! Exactly right only at z=p=2
             xpd2=pd2-pd*y2*(b1+y2*(b2+y2*(b3+y2*(b4+y2*(b5+y2*b6))) ))  ! The gem formula was better even at
             x0=p*(1.0_np-exp(-xpd2))                        ! z=p=2 where this is supposed to be
!             print*,'pd2,y2,xpd2,x0'                        ! the exact large series value.
!             print*,'pd2,y2,xpd2,x0,1.0_np'                 ! Well the large series isn't that good at 2.
!             print*,pd2,y2,xpd2,x0,1.0_np
!             stop
          else if(isatz .eq. 10) then                   ! Briiliant flop formula.
            if(p .lt. 2.0_np) then                       ! Spectacular error. 
                w=pd**2
              else
                w=1.0_np
            end if
            wlow=max(1.0_np-w,0.0_np)  ! Should never be numerically negative.
            c0=p*w  ! The low series coeficient is 0 and the high is z=1+delz=p.
            c1=-a1*wlow+b1*w   ! The minus are because the it's natural to put a minus in front of the 
            c2=-a2*wlow+b2*w   ! whole nonzero order set of terms.
            c3=-a3*wlow+b3*w
            c3=-a4*wlow+b4*w
            c4=-a5*wlow+b5*w
            c6=-a6*wlow+b6*w
            s=pd*wlow+y*w     ! The interpolation expansion parameter.
            x0=c0-s*(c1+s*(c2+s*(c3+s*(c4+s*(c5+s*c6)) )))
              xlow=pd*(a1+pd*(a2+pd*(a3+pd*(a4+pd*(a5+pd*a6
     &                       +pd*a7))) ))     ! High precision form.  ! a7=0 except for tests
              xhigh=p-y*(b1+y*(b2+y*(b3+y*(b4+y*(b5+y*b6))) ))   ! High precision form. 
              print*,'pd,y,s'
              print*,pd,y,s
        end if
        if(icorrect .eq. 0) then
            xcor=0.0_np
          else if(icorrect .eq. 1) then
            xcor=p*(1.0_np-exp(-x0))                           ! Interation formula correction 
            x0=xcor
          else if(icorrect .eq. 2) then
!            xcor=p*(1.0_np-2.0_np*exp(-x0))/(1.0_np-p*exp(-x0))       ! Newton-Raphson correction, numerically rotten
            xcor=x0-( (x0-p*(1.0_np-exp(-x0)))/(1.0_np-p*exp(-x0)) ) ! Newton-Raphson equation 
            x0=xcor
        end if 
        print*,'xlow,xhigh,(xlow-xhigh)/xhigh,xcor,x0'
        print*,xlow,xhigh,(xlow-xhigh)/xhigh,xcor,x0
!
        x1=x0  ! x1 for the iteration equation 
        x2=x0  ! x2 for the Newton-Raphson 
        jcon2a=0   ! Iterations to single precision convergence
        jcon2b=0   ! Iterations to machine precision convergence 
        dcon1=1.d6  ! Not used at present, but for variable symmetry.
        dcon2=1.d6  ! Tests to see if Newton-Raphson has reached jitter stage
!                   ! implying machine precision convergence.
        j=0
        write(*,990) j,p,x1,x2  ! Initial values
        x1old=x1
        x2old=x2
        kcon=10               ! Maximum number of iterations.  If more needed, something is wrong.
        do k=1,kcon
          j=j+1
          x1=p*(1.0_np-exp(-x1old))    ! Iteration equation
          x2=x2old-( (x2old-p*(1.0_np-exp(-x2old)))                        &
     &               /(1.0_np-p*exp(-x2old)) ) ! Newton-Raphson equation
!          x2=p*(1.0_np-2.0_np*exp(-x2old))/(1.0_np-p*exp(-x2old)) ! NR equation compact, but numerically rotten
          write(*,990) j,p,x1,x2
          if(abs(x2-x2old) .le. 1.0e-8_np*x2                              &
     &         .and. jcon2a .eq. 0) jcon2a=j ! Single precision test
          if(abs(x2-x2old) .ge. dcon2  .or.  jcon2b .gt. 0) then        ! Machine-precsion/jitter tests
              if(jcon2b .eq. 0) jcon2b=j
              if(k .ge. 4) exit                                       ! Always do at least 4 iterations
          end if
          dcon1=abs(x1-x1old)
          dcon2=abs(x2-x2old)   ! For machine-precision/jitter test
          x1old=x1
          x2old=x2
        end do
!        write(*,990) j,p,x1,x2
        print*,'NR iteration, relerr in fit',
     &            ' relerr in it-function=',kcon
        print*,jcon2a,jcon2b,(x0-x2)/x2,abs(x1-x2)/x2
        if(isatz .ge. 2) then
           print*,'xlow-x2,xhigh-x2,x0-x2,x2'
           print*,xlow-x2,xhigh-x2,x0-x2,x2
        end if
        if(itest .eq. 1  .and.                                             &
     &     (
     &      abs(pd-1.0e-6_np) .lt. 1.0e-10_np .or.                               &  precision 32 looses marbles
     &      abs(pd-1.0e-5_np) .lt. 1.0e-10_np .or.                               &  precision 32 looses marbles
     &      abs(pd-0.0001_np) .lt. 1.0e-10_np .or.                               &  precision 32 looses marbles
!     &      abs(pd-0.003_np) .lt. 0.0001_np .or.                               &  precision 32 looses marbles
     &      abs(pd-0.01_np) .lt. 0.001_np .or.                               &  
!     &      abs(pd-0.03_np) .lt. 0.001_np .or.                               &  
     &      abs(pd-0.1_np) .lt. 0.001_np .or.                                &
     &      abs(pd-0.13_np) .lt. 0.001_np .or.                                &
     &      abs(pd-0.15_np) .lt. 0.001_np .or.                                &
     &      abs(pd-0.2_np) .lt. 0.001_np .or.                                &
     &      abs(pd-0.3_np) .lt. 0.001_np .or.                                &
     &      abs(pd-0.5_np) .lt. 0.001_np .or.                                & 
     &      abs(pd-0.8_np) .lt. 0.001_np .or.                                & 
     &      abs(pd-1.0_np) .lt. 0.001_np                                     & 
     &                                        )     ) then   
            if(abs(pd-1.0e-6_np) .lt. 1.0e-10_np) then
                x2_in=1.99999933333377777745183802479528595E-0006_np
                print*,'Replace yes/no?'
                print*,x2,x2_in,x2-x2_in
!                x2=x2_in
            end if
            print*,'For Testing for the low series coefficients.'
            xcoef2=(x2-xlow)/pd**2 + a2
            xcoef3=(x2-xlow2)/pd**3 
            xcoef4=(x2-xlow3)/pd**4 
            xcoef5=(x2-xlow4)/pd**5
            xcoef6=(x2-xlow5)/pd**6
            xcoef7=(x2-xlow)/pd**7 + a7 
            xcoef8=(x2-xlow)/pd**8 + a8
            xcoef9=(x2-xlow)/pd**9 + a9
            print*,'xcoef2,a2,(a2-xcoef2)/xcoef2'
            print*,xcoef2,a2,(a2-xcoef2)/xcoef2
            print*,'xcoef3,a3,(a3-xcoef3)/xcoef3'
            print*,xcoef3,a3,(a3-xcoef3)/xcoef3
            print*,'xcoef4,a4,(a4-xcoef4)/xcoef4'
!
            print*,xcoef4,a4,(a4-xcoef4)/xcoef4
            print*
            print*,'xcoef4*9.0_np,xcoef4*18.0_np'
            print*,xcoef4*9.0_np,xcoef4*18.0_np
            print*,'xcoef4*27.0_np,xcoef4*81.0_np'
            print*,xcoef4*27.0_np,xcoef4*81.0_np
            print*,'xcoef4*11.0_np,xcoef4*14.0_np'
            print*,xcoef4*11.0_np,xcoef4*14.0_np
            print*,'xcoef4*21.0_np,xcoef4*24.0_np'
            print*,xcoef4*21.0_np,xcoef4*24.0_np
            xc=abs(xcoef4)
            xsign=sign(1.0_np,xcoef4)
            in=0
            if(in .eq. 1) then
                do k=1,2000
                xd=real(k)
                xn=xd*xc
                xdiff=xn-real(int(xn+0.5_np))
                relerr=xsign*(xdiff/xd)/pd
                itell=0
                iprime=nprime
                do kk=2,iprime   ! Must have one prime factor in this range 
                   jtell=mod(k,iiprime(kk)) 
                   if(jtell .eq. 0) itell=itell+1
                end do
                if(abs(relerr) .lt. 20.0_np                                         &
     &            .and. itell .gt. 0)
     &            print*,itell,k,xn,xdiff,relerr
                end do
                stop
             end if
!
            print*
            print*,'xcoef5,a5,(a5-xcoef5)/xcoef5'
            print*,xcoef5,a5,(a5-xcoef5)/xcoef5
            xc=abs(xcoef5)
            xsign=sign(1.0_np,xcoef5)
            in=0
            if(in .eq. 1) then
                do k=1,10 000
                xd=real(k)
                xn=xd*xc
                xdiff=xn-real(int(xn+0.5_np))
                relerr=xsign*(xdiff/xd)/pd
                if(abs(relerr) .lt. 20.0_np) print*,k,xn,xdiff,relerr
                end do
                stop
             end if
!
            print*,'xcoef6,a6,(a6-xcoef6)/xcoef6'
            print*,xcoef6,a6,(a6-xcoef6)/xcoef6
!
            print*
            print*,'Remainder ',(x2-xlow)/pd**6
            print*
            print*,'xcoef7,a7,(a7-xcoef7)/xcoef7'
            print*,xcoef7,a7,(a7-xcoef7)/xcoef7
            print*,'xcoef7*9.0_np,xcoef7*18.0_np'
            print*,xcoef7*9.0_np,xcoef7*18.0_np
            print*,'xcoef7*27.0_np,xcoef7*81.0_np'
            print*,xcoef7*27.0_np,xcoef7*81.0_np
!
            print*
            print*,'xcoef8,a8,(a8-xcoef8)/xcoef8'
            print*,xcoef8,a8,(a8-xcoef8)/xcoef8
            print*,'xcoef8*9.0_np,xcoef8*18.0_np'
            print*,xcoef8*9.0_np,xcoef8*18.0_np
            print*,'xcoef8*27.0_np,xcoef8*81.0_np'
            print*,xcoef8*27.0_np,xcoef8*81.0_np
!
            print*
            print*,'xcoef9,a9,(a9-xcoef9)/xcoef9'
            print*,xcoef9,a9,(a9-xcoef9)/xcoef9
            print*,'xcoef9*9.0_np,xcoef9*18.0_np'
            print*,xcoef9*9.0_np,xcoef9*18.0_np
            print*,'xcoef9*27.0_np,xcoef9*81.0_np'
            print*,xcoef9*27.0_np,xcoef9*81.0_np
        end if
      end do
  990 format(i10,f10.5,2f20.16)
!
                             ! http://physics.nist.gov/cuu/Constants/Table/allascii.txt
      boltev=8.6173303e-5_np  ! 8 digits boltzmann's constant NIST
      ryd=13.6056930090_np   ! 11 digits rydberg NIST
!       ryd=13.605692530_np   ! http://en.wikipedia.org/wiki/Rydberg_constant
      t=ryd/(boltev*x)
      print*,'Comparison of Liddle reionization temperature and mine.'
      print*,5700.0_np,t
!
      end

