c integrates the equations for polytropic equilibrium from theta=pi/2
c to theta=0 
c last modified 26.8.96
      external derivs, rkqc
      parameter (npt=200,nvar=4)
      dimension ystart(nvar)
      common/par/ gm, cc, pi
      pi=2.*asin(1.)
c----------------------------------------------------------
      gm=1.       ! gamma
      cc=0.5      ! H_0
c----------------------------------------------------------
C Parameters for the integration routine
      eps=1.0e-4
      h1=0.005
      hmin=0.0
c ------------------------------------------------------------------------
c Initial conditions
c ------------------------------------------------------------------------
      write (6,*) 'enter y1(pi/2)'
      read (5,*) a0
      write (6,*) 'enter y3(pi/2)'
      read (5,*) b0
      tin=asin(1.)  ! initial angle, in rad
      ystart(1)=a0  
      ystart(2)=.0
      ystart(3)=b0  
      ystart(4)=.0
c ------------------------------------------------------------------------
c Output parameters
      tfin = .0   !final angle, in rad
c------------------------------------------------------------------------
      t=tin
      bc=sin(t)*(2.*cc*ystart(1)**((2.-gm)/(3.*gm-4.))*ystart(2)
     *   -ystart(3)**(gm-2.)*ystart(4)/(2.*pi)**(gm-1.))
      write(6,20) t/tin,ystart(1),ystart(2),ystart(3),ystart(4),bc
      write(67,20) t/tin,ystart(1),ystart(2),ystart(3),ystart(4),bc
      do 40 iout = 1,npt-1
c ------------------------------------------------------------------------
       tout=tin+(tfin-tin)*iout/npt
       call odeint(ystart,nvar,tin,tout,eps,h1,hmin,nok,
     *             nbad,derivs,rkqc)
       bc=sin(t)*(2.*cc*ystart(1)**((2.-gm)/(3.*gm-4.))*ystart(2)
     *   -ystart(3)**(gm-2.)*ystart(4)/(2.*pi)**(gm-1.))
       write(6,20) tout/asin(1.),ystart(1),ystart(2),
     *             ystart(3),ystart(4),bc
       write(67,20) tout/asin(1.),ystart(1),ystart(2),
     *              ystart(3),ystart(4),bc
c three-point extrapolation to theta=0
      if (iout.eq.npt-3) then
       y1c=ystart(1)
       y2c=ystart(2)
       y3c=ystart(3)
       y4c=ystart(4)
       bcc=bc
      endif
      if (iout.eq.npt-2) then
       y1b=ystart(1)
       y2b=ystart(2)
       y3b=ystart(3)
       y4b=ystart(4)
       bcb=bc
      endif
      if (iout.eq.npt-1) then
       y1a=ystart(1)
       y2a=ystart(2)
       y3a=ystart(3)
       y4a=ystart(4)
       bca=bc
      endif
  40  continue
      to=0.0
      y1o=3.*y1a-3.*y1b+y1c
      y2o=3.*y2a-3.*y2b+y2c
      y3o=3.*y3a-3.*y3b+y3c
      y4o=3.*y4a-3.*y4b+y4c
      bco=3.*bca-3.*bcb+bcc
      write(6,20) to,y1o,y2o,y3o,y4o,bco
      write(67,20) to,y1o,y2o,y3o,y4o,bco
c ------------------------------------------------------------
  20  format(1x,1pe12.4,3x,1p5e12.4)
c ------------------------------------------------------------------------
      stop
      end
C 
      subroutine derivs (t, y, dydt)
      parameter (nvar=4)
      dimension y(nvar), dydt(nvar)
      common/par/ gm, cc, pi
c ------------------------------------------------------------------------
      dydt(1)=y(2)
      dydt(2)=y(2)*cos(t)/sin(t)-2.*(3.*gm-4.)*(gm-1.)/(gm-2.)**2*y(1)
     * -cc*y(3)*y(1)**((2.-gm)/(3.*gm-4.))*sin(t)**2
      dydt(3)=y(4)
      dydt(4)=-(2.*pi)**(gm-1.)*(
     *  2.*y(3)**(3.-gm)  
     * +2.*(3.*gm-4.)/(gm-2.)**2*y(3)/(2.*pi)**(gm-1.) 
     * -2.*((3.*gm-4.)/(gm-2.))**2*cc*y(1)**(2.*(gm-1.)/(3.*gm-4.)) 
     *        *y(3)**(2.-gm)
     * -2.*cc*cos(t)/sin(t)*y(1)**((2.-gm)/(3.*gm-4.))*y(2)**2
     *        *y(3)**(2.-gm)
     * +1./(2.*pi)**(gm-1.)*cos(t)/sin(t)*y(4)
     * -2.*cc*(2.-gm)/(3.*gm-4.)*y(1)**((6.-4.*gm)/(3.*gm-4.))*y(2)**2
     *        *y(3)**(2.-gm)
     * -2.*cc*y(1)**((2.-gm)/(3.*gm-4.))*dydt(2)*y(3)**(2.-gm)
     * +(gm-2.)/(2.*pi)**(gm-1.)*y(4)**2/y(3) )
c ------------------------------------------------------------------------
      return
      end