!/ module fsig
!.
!.               %%%%%%%%%%%%%%%%%%%%%%%%%
!.               %                       %
!.               %       fsig            %
!.               %                       %
!.               %%%%%%%%%%%%%%%%%%%%%%%%%
!.
function fsig(zeb,ef,af,ksig)
  !  25-mar-81 added ksig=5 option
  !  27-may-82    corrected gryzinski proton ionization formula
  !============
  use iso_c_binding, only: fp => c_double
  implicit none
  integer :: ksig
  !============
  real(fp) :: zeb,ef,af,fsig,za0,za1,za2,za3,za4,za5,za6,zp0,zp1,zp2
  real(fp) :: zp3,zp4,zp5,zp6,fg2,z,galf,ze,zlgsig,zx,zz,blimt,ulimt
  real(fp) :: alph
  !============
  data za0,za1,za2,za3,za4,za5,za6/-3.17385e1_fp,1.143818e1_fp, &
       -3.833998_fp,0.7046692_fp,-7.431486e-2_fp,4.153749e-3_fp, &
       -9.486967e-5_fp/
  data zp0,zp1,zp2,zp3,zp4,zp5,zp6/-42.03309_fp,3.557321_fp, &
       -1.045134_fp,0.3139238_fp,-7.454475e-2_fp,8.459113e-3_fp, &
       -3.495444e-4_fp/
  !  corrected fg2: -1/galf before alog
  fg2(z)=(z/(z+1._fp))**1.5_fp*(1._fp/z)*(z/(z+1._fp)+ &
       (2._fp/3._fp)*(1._fp-1._fp/galf)*log(2.7_fp+sqrt(z)))* &
       (1._fp-1._fp/galf)*(1._fp-(1._fp/galf)**(1+z))
  !
  !-----------------------------------------------------------------------
  !
  fsig=0.0_fp
  !
  if(ksig .eq. 0) then
     !        proton ionization from freeman and jones.
     ze= log(ef/(af*1000._fp))
     zlgsig= zp0 +ze*(zp1+ze*(zp2+ze*(zp3+ze*( &
          zp4+ze*(zp5+ze*zp6)))))
     fsig= exp(zlgsig)
     return
  end if
  if(ksig .eq. 1) then
     !        proton ionization from gryzinski
     zx=0.0233_fp*sqrt(ef/(af*zeb))
     zz=zx**2
     galf=4._fp*zz*(1._fp+1._fp/zx)
     fsig=8.79e-17_fp*746._fp*fg2(zz)/zeb**2
     if(galf.le.1._fp) fsig=0._fp
     return
  end if
  if(ksig .eq. 2) then
     !        proton ionization from garcia,gerjouy, and welker
     blimt= (sqrt(2._fp)-1._fp)/2._fp
     ulimt= blimt +1._fp
     alph= .0233_fp*sqrt(ef/(af*zeb))
     if(alph .lt. blimt) fsig= 0._fp
     if(alph .lt. blimt) return
     if(alph .gt. ulimt) then
        fsig= (5._fp+3._fp/(8._fp*(alph+1._fp)) + &
             alph/(8._fp*(alph-1._fp)**2) &
             -(4._fp*alph-1._fp)/(8._fp*alph**2) -1._fp/(8._fp*alph* &
             alph*(alph-1._fp)**2) -3._fp/(4._fp*alph*(alph-1._fp)))/ &
             (3._fp*(alph**2)*(zeb**2))
     else
        fsig=(-0.328_fp+2._fp*alph**3 +3._fp/(8._fp*(alph+1._fp)))/ &
             (3._fp*(alph**2)*(zeb**2))
     endif
     fsig=fsig*740._fp*8.79e-17_fp
     return
  end if
  if(ksig .eq. 3) then
     !        electron ionization (sigma-v-bar) from freeman and jones
     ze=log(ef*(13.6_fp/zeb))
     zlgsig=za0+ze*(za1+ze*(za2+ze*(za3+ze* &
          (za4+ze*(za5+ze*za6)))))
     fsig=exp(zlgsig)*(13.6_fp/zeb)**1.5_fp
     return
  end if
  if(ksig .eq. 4) then
     !        hydrogen charge-transfer from freeman and jones
     ze= log(ef/af)/2.30258_fp
     fsig= 6.937e-15_fp*(1._fp-0.155_fp*ze)**2/ &
          (1._fp+1.112e-15_fp*(ef/af)**3.3_fp)
     return
  end if
  if(ksig .eq. 5) then
     !        electron loss by impurity ions
     !        from olson. et al., phys. rev. letters, 41 (1978) 163
     ze=ef*13.6_fp/(1000._fp*af*zeb*32._fp)
     fsig=4.6e-16_fp*(13.6_fp/zeb)**2*(1._fp-exp(-ze))/ze
     return
  end if
end function fsig
