;   This file checked against the Univac version, July 1987.  MA
;   converted to IDL, Feb 1999 lmw
      FUNCTION AI0, X
;
;   CALCULATES MODIFIED BESSEL FUNCTION OF FIRST KIND, I0(X).
;   USES FORMULAE 9.8.1 AND 9.8.2 FROM ABRAMOWITZ & STEGUN
;
      T = X/3.75
      T2 = T*T
      IF (X lt 3.75) then begin
       A = 1. + T2*(3.5156229 + T2*(3.0899424 + T2*(1.2067492 + $
        T2*(.2659732 + T2*(.0360768 + .0045813*T2)))))

      endif else begin
       TEMP = .39894228 + (.01328592 + (.00225319 + (-.00157565 + $
           (.00916281 + (-.02057706 + (.02635537 + (-.01647633 + $
            .00392377/T)/T)/T)/T)/T)/T)/T)/T
       if (X GT 85.) then A =  temp else A = TEMP*EXP(X)/SQRT(X)
      endelse
     RETURN, A
      END