function cion(n,j,e,t)
c
c     collisional ionization rate from Raymond
c     rate is returned in units cm^3/s
c     inputs are: n=nuclear charge,
c                 j=ion stage
c                 e=valence shell ionization threshold (in eV)
c                 t=temperature in K.
c
c Tim Kallman's routine for calculating collisional ionization rates.
c    Note that this routine only accounts for valence shell ionization.
c    It should be called only once for each ion stage and the valence
c    shell ionization threshold should be the lowest one (e.g. the 2p
c    ionization potential for CI).
c   
c
      INCLUDE 'IMPLIC.FOR'
c
c  sm younger jqsrt 26, 329; 27, 541; 29, 61   with moores for undone
c  a0 for b-like ion has twice 2s plus one 2p  as in summers et al
c  chi = kt / i
c
      dimension a0(30),a1(30),a2(30),a3(30),b0(30),b1(30),
     &          b2(30),b3(30),c0(30),c1(30),c2(30),c3(30),
     &          d0(30),d1(30),d2(30),d3(30)
c
      data a0/13.5,27.0,9.07,11.80,20.2,28.6,37.0,45.4,
     &     53.8,62.2,11.7,38.8,37.27,46.7,57.4,67.0,
     &     77.8,90.1,106.,120.8,135.6,150.4,165.2,180.0,
     &     194.8,209.6,224.4,239.2,154.0,268.8/
      data a1/ - 14.2,-60.1,4.30,27*0./
      data a2/40.6,140.,7.69,27*0./
      data a3/ - 17.1,-89.8,-7.53,27*0./
c
      data b0/ - 4.81,-9.62,-2.47,-3.28,-5.96,-8.64,-11.32,
     &     -14.00,-16.68,-19.36,-4.29,-16.7,-14.58,-16.95,
     &     -19.93,-23.05,-26.00,-29.45,-34.25,-38.92,
     &     -43.59,-48.26,-52.93,-57.60,-62.27,-66.94,
     &     -71.62,-76.29,-80.96,-85.63/
      data b1/9.77,33.1,-3.78,27*0./
      data b2/ - 28.3,-82.5,-3.59,27*0./
      data b3/11.4,54.6,3.34,27*0./
c
      data c0/1.85,3.69,1.34,1.64,2.31,2.984,3.656,4.328,
     &     5.00,5.672,1.061,1.87,3.26,5.07,6.67,8.10,
     &     9.92,11.79,7.953,8.408,8.863,9.318,9.773,
     &     10.228,10.683,11.138,11.593,12.048,12.505,12.96/
      data c1/0.,4.32,.343,27*0./
      data c2/0.,-2.527,-2.46,27*0./
      data c3/0.,.262,1.38,27*0./
c
      data d0/ - 10.9,-21.7,-5.37,-7.58,-12.66,-17.74,
     &     -22.82,-27.9,-32.98,-38.06,-7.34,-28.8,-24.87,
     &     -30.5,-37.9,-45.3,-53.8,-64.6,-54.54,-61.70,
     &     -68.86,-76.02,-83.18,-90.34,-97.50,-104.66,
     &     -111.82,-118.98,-126.14,-133.32/
      data d1/8.90,42.5,-12.4,27*0./
      data d2/ - 35.7,-131.,-8.09,27*0./
      data d3/16.5,87.4,1.23,27*0./
c
      cion = 0.
      chir = t/(11590.*e)
      if ( chir.le..0115 ) return
      chi=chir
      if(chi.lt.0.1) ch=0.1
      ch2 = chi*chi
      ch3 = ch2*chi
      alpha = (.001193+.9764*chi+.6604*ch2+.02590*ch3)
     &        /(1.0+1.488*chi+.2972*ch2+.004925*ch3)
      beta = (-.0005725+.01345*chi+.8691*ch2+.03404*ch3)
     &       /(1.0+2.197*chi+.2457*ch2+.002503*ch3)
      j2 = j*j
      j3 = j2*j
      iso = n - j + 1
c
      a = a0(iso) + a1(iso)/j + a2(iso)/j2 + a3(iso)/j3
      b = b0(iso) + b1(iso)/j + b2(iso)/j2 + b3(iso)/j3
      c = c0(iso) + c1(iso)/j + c2(iso)/j2 + c3(iso)/j3
      d = d0(iso) + d1(iso)/j + d2(iso)/j2 + d3(iso)/j3
c
c  fe ii experimental ionization montague et al: d. neufeld fit
      if ( n.eq.26 .and. j.eq.2 ) then
         a = -13.825
         b = -11.0395
         c = 21.07262
         d = 0.
      endif
c
      ch = 1./chi
      fchi = 0.3*ch*(a+b*(1.+ch)+(c-(a+b*(2.+ch))*ch)*alpha+d*beta*ch)
      cion = 2.2e-6*sqrt(chir)*fchi*expo(-1./chir)/(e*sqrt(e))
      return
      end