******************************************************************************
      program pol4
c ----- integrates a polytrope with parameters to be typed from the terminal
c ----- displays the results on the screen and stores on the disk in a file
c ----- "polytrop.res"
c CONSTANTS:
c   fn  =  n  = polytropic index
c VARIABLES:
c   x(1)  = dimensionless polytropic radius
c   x(2)  = dimensionless polytropic mass
c   x(3)  = dimensionless polytropic "temperature"
      implicit double precision (a-h,o-z)
      common/const/fn
      dimension x(10)

      write(*,*)' program "pol4.f", results ==> "pol4.res" '
      open(1,file='pol4.res')
      write(1,105)
      write(*,105)
  105 format('   is  -  number of the current integration step',/,
     *       '   x1  -  dimensionless polytropic radius',/,
     *       '   x2  -  dimensionless polytropic mass',/,
     *       '   x3  -  dimensionless polytropic "temperature"')
      x3min=1.0d-3
c       x3min = the minimum value of the variable "x3", at which the
c               outward integrations are stopped; the region beyound
c               is integrated with analytical expansion

    3 continue
c ----- input the parameters for the polytrope from the keyboard:
c   fn  =  n  = polytropic index
c   h         = integration step size
c   ni        = print control, when ni<0 no printout

      write(*,101)
  101 format('         polytropic index:   n = ',$)
      read(*,*)fn
      if(fn.le.0)close(1)
      if(fn.le.0)stop
      write(*,102)
  102 format('    integration step size:   h = ',$)
      read(*,*)h
      write(*,103)
  103 format('  output every "ni" steps:  ni = ',$)
      read(*,*)ni

c ----- expand variables in a power series near the stellar center, with
c ----- the expansion radius three times larger than the integration step:
      x(1)=h*3
      x(2)=x(1)**3/3*(1-fn/10*x(1)**2)
      x(3)=1-x(1)**2/6*(1-fn/20*x(1)**2)

c ----- prepare the counters:
      is=2
      ip=-ni+2
      write(1,106)
      write(*,106)
  106 format('     is    x1      x2      x3')

c ----- integrate out untill the surface is close :
    1 continue
      is=is+1
      if(is.gt.100000)go to 3
      ip=ip+1
      if(ip.lt.0)go to 2

c ----- print the results:
      ip=-ni
      write(1,100)is,x(1),x(2),x(3)
      write(*,100)is,x(1),x(2),x(3)
  100 format(1x,i6,3f8.4)
    2 continue

c ----- make one integration step:
      x2p=x(2)
      x3p=x(3)
      call step(x,h,3)
      if(x(3).gt.x3min)go to 1

c ----- polytopic "temperature" x(3) < x3min
c ----- interpolate the dimensionless radius "r" to x(3)=x3min,
      r=x(1)+h*(x3min-x(3))/(x3p-x(3))
      fm=x(2)+(x2p-x(2))*(x3min-x(3))/(x3p-x(3))
      t=x3min

      isn=-is
      write(1,104)isn,r,fm,t,fn
      write(*,104)isn,r,fm,t,fn
c       r = dimensionless radius at the end of numerical integrations
c      fm = dimensionless mass at the end of numerical integrations
c       t = x3max = dimensionless "temperature" at the end of numerical integ.

c expand the solution all the way to the surface using analytic approximation
      rs=1/(1/r-t/fm)
      fms=fm*(1+r**4*t**(fn+1)/(fn+1))
      r=rs
      fm=fms
      t=0.0d0
c       r = dimensionless radius at the surface, i.e. total radius
c      fm = dimensionless mass at the surface, i.e. total mass
c       t = 0 = dimensionless "temperature" at the surface
      
      write(1,104)isn,r,fm,t,fn
      write(*,104)isn,r,fm,t,fn
  104 format(1x,i6,3f8.4,'    polytropic index: n =',f8.4)

c ----- go for a new set of input parameters:
      go to 3

      end
******************************************************************************
******************************************************************************
      subroutine step(x,h,k)
c ------- subroutine step makes one integration step with a fourth order Runge
c ------- Kutta method in a double precision. It assumes that the right hand
c ------- sides are calculated with a subroutine "rhs".

c INPUT:
c           x(1), x(2), ... x(k) = variables at the beginning of the step
c           h = integration step in variable  x(1)
c           k = the number of variables, required to be no more than 20

c OUTPUT:
c           x(1), x(2), ... x(k) = variables at the end of the step

c NOTICE: the subroutine "rhs" has "x" array as input and "y" array as output,
c           y(i) = d x(i) / d x(1) ,  i = 1, 2, ... k

      implicit double precision (a-h,o-z)
      dimension x(20),y(20),xp(20),xx(20)

      call rhs(x,y)

      do 1 i=1,k 
      xp(i)=x(i)+h*y(i)/2
      xx(i)=x(i)+h*y(i)/6
    1 continue

      call rhs(xp,y)

      do 2 i=1,k
      xp(i)=x(i)+h*y(i)/2
      xx(i)=xx(i)+h*y(i)/3
    2 continue

      call rhs(xp,y)

      do 3 i=1,k 
      xp(i)=x(i)+h*y(i)
      xx(i)=xx(i)+h*y(i)/3
    3 continue

      call rhs(xp,y)

      do 4 i=1,k
      x(i)=xx(i)+h*y(i)/6
    4 continue

      return
      end
******************************************************************************
******************************************************************************
      subroutine rhs(x,y)
c       calculates derivatives, i.e. right hand sides of the differential
c       equations
c INPUT:
c   x(1)  = dimensionless polytropic radius
c   x(2)  = dimensionless polytropic mass
c   x(3)  = polytropic "temperature"
c   fn    = polytropic index   (in common/const/)
c OUTPUT:
c   y(i)  = dx(i)/dx(1),   i = 1, 2, 3
      implicit double precision (a-h,o-z)
      dimension x(10),y(10)
      common/const/fn

      y(1)=1.0
      y(2)=x(1)**2*dabs(x(3))**fn
      y(3)=-x(2)/x(1)**2

      return
      end
******************************************************************************