c-------------------------------------------------------------------- 
c                      Vcirc and g(r) for a thin disk
c
      dr=0.025
      x =dr
      do while (x.le.5.)
         vv =sqrt(VC2(x))
         g  =dMass(x)/x**2
         write (*,'(6f9.4)') x,vv,vv**2/x,g,(vv**2/x)/g
         x = x+dr
         if(x.gt.1.)dr=0.1
      enddo

      stop
      end

C--------------------------------------------------------------
      REAL*4 function VC2(x)  !Disk Vcircular**2 
C--------------------------------------------------------------
      External  bessi0,bessk0,bessi1,bessk1 
      y=.5*x  
      VC2 =2.*y*y*( bessi0(y) *bessk0(y)-bessi1(y) *bessk1(y))   
   
      Return  
      end
C---------------------------------------------------------------
      REAL*4 Function bessi0(x) ! Modified Bessel I_o Function  
C-----------------------------------------------------------  
      IMPLICIT REAL*4(A-H,O-Z) 
   
c      Real*4 x                  !bessi0, 
C     Returns the modified Bessel Function  I_o(x) for any real x
c      Real ax
      Real*8 p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y ! Accumulate polinomials in double precision
      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/1.0d0,3.5156229d0,3.0899424d0,
     *     1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
      Data q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,
     *     0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,
     *     0.2635537d-1,-0.1647633d-1,0.392377d-2/ 
      
      if (abs(x).lt.3.75)then
         y=(x/3.75)**2
         bessi0 = p1 + y*(p2+y*(p3 + y*(p4 + y*(p5 + y*(p6 +y*p7)))))
      else
         ax = abs(x)
         y = 3.75/ax
         bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 +y*(q3 + y*(q4 + 
     *        y*(q5 +y*(q6 + y*(q7 + y*(q8 + y*q9)))))))) 
      endif

      Return
      End 


C------------------------------------------------------------
      REAL*4 Function bessk0(x) ! Modified Bessel K_o Function 
C----------------------------------------------------------- 
 
C     Returns the modified Bessel Function  K_o(x) for positive real x  
      IMPLICIT REAL*4(A-H,O-Z)
c      Real*4 x                  
      external bessi0  
      Real*8 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/-0.57721566d0,0.42278420d0,0.23069756d0,
     *     0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/ 
      Data q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1,
     *     -0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/ 
      
      If(x.le.2.0)then  
         y = x*x/4.0
         bessk0 = (-log(x/2.0)* bessi0(x)) + (p1 + y*(p2 +y*(p3 + 
     *        y*(p4 + y*(p5 + y*(p6 + y*p7))))))
      Else 
         y = (2.0/x)
         bessk0 = ( exp(-x)/sqrt(x))*(q1 + y*(q2 + y*(q3 +
     *        y*(q4 + y*(q5 + y*(q6 + y*q7))))))
      Endif
      Return 
      End
C------------------------------------------------------------
      REAL*4 Function bessi1(x) ! Modified Bessel I_1 Function 
C----------------------------------------------------------- 
      IMPLICIT REAL*4(A-H,O-Z)
c      Real*4 x   
C     Returns the modified Bessel Function  I_o(x) for any real x 
c      Real ax
      Real*8 p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y ! Accumulate polinomials in double precision
      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.75)then
         y = (x/3.75)**2  
         bessi1 = x*(p1 + y*(p2 +y*(p3+ y*(p4+ y*(p5+ y*(p6 + y*p7)))))) 
      Else
         ax = abs(x) 
         y = 3.75/ax
         bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 +
     *        y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))))
         if(x.lt.0.)bessi1= -bessi1
      Endif
      Return
      End

C------------------------------------------------------------
      REAL*4 Function bessk1(x) ! Modified Bessel K_1 Function 
C----------------------------------------------------------- 
      IMPLICIT REAL*4(A-H,O-Z)
c      Real*4 x   
c      Real*4 bessi1 
      external bessi1 
      Real*8 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.0)then
         y = x*x/4.0
         bessk1 = (log(x/2.0)*bessi1(x)) + (1.0/x)*(p1 + y*(p2+
     *        y*(p3+y*(p4+y*(p5 + y*(p6 + y*p7)))))) 
      Else
         y = 2.0/x  
         bessk1 = ( exp(-x)/sqrt(x))*(q1 + y*(q2+y*(q3 +
     *        y*(q4+y*(q5 + y*(q6 + y*q7))))))
      Endif
      Return
      End
  
C----------------------------------------------------------
      REAL*4 Function dMass(x) ! This is the dimensionless disk mass. 
C---------------------------------------------------------- 
         dMass = (1.0 - EXP(-X)*(1.0 + X))
      Return  
      End