C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2005
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c     
c-----------------------------------------------------------------------
c     $Log: utils.F,v $
c     Revision 1.1.1.1  2005/09/12 17:28:10  samtaney
c     Original source.
c
c     Revision 1.1  2005/05/11 13:15:34  samtaney
c     Utility functions.
c
c     Original source for utilities 
c-----------------------------------------------------------------------
c=======================================================================
      double precision FUNCTION bessk1(x)
      double precision:: x
CU    USES bessi1
      double precision,external:: bessi1
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
      DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0,
     *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
      DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1,
     *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
      if (x.le.2.D0) then
        y=x*x/4.D0
        bessk1=(dlog(x/2.D0)*bessi1(x))+
     &       (1.D0/x)*(p1+y*(p2+y*(p3+y*(p4+y*
     &       (p5+y*(p6+y*p7))))))
      else
         y=2.D0/x
         bessk1=(dexp(-x)/dsqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     &       q7))))))
      endif
      return
      END
c
c-----------------------------------------------------------------------
      double precision FUNCTION bessi1(x)
      double precision:: x
      double precision:: ax
      DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
      SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
      DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,
     *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
      DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1,
     *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1,
     *0.1787654d-1,-0.420059d-2/
      if (abs(x).lt.3.75D0) then
        y=(x/3.75D0)**2
        bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
      else
        ax=abs(x)
        y=3.75D0/ax
        bessi1=(dexp(ax)/dsqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
     &       (q7+y*(q8+y*q9))))))))
        if(x.lt.0.D0)bessi1=-bessi1
      endif
      return
      END