program solardensity
C
C  this program makes a crude estimate of the density profile 
C  of a star
C  with polytrope of gamma 
C  with both explicit Euler and taylor integration
C
C  uses scaled quantities
C 
C  integrates out rho(r), m(r)
C
C  drho/dr = - (1/Gamma) rho**(2-gamma) *m/ r^2
C  d m/dr  = 3 r^2 rho
C
C  CWJ SDSU April 2005 revised Oct 2006
C

      implicit none

      real rho,mass,r0
      real rho0,mass0
      real rho00   ! starting density
      real dr
      integer i,n 

      real gamma

      print*,' enter index Gamma '
      read*,gamma
      print*,' Enter starting density (of order 1) '
      read*,rho00

      print*,' Enter dr '
      read*,dr
      n = int(20./dr)+1
 
      r0 = dr


C.............Explicit Euler.......

      print*,' '
      print*,' Euler '
      rho0 = rho00
      mass0 = rho0*4*3.14159*dr**3/3.
      do i = 1, n+1
        print*,r0,rho0,mass0
        write(16,*)r0,rho0
        mass = mass0+dr*3*r0**2*rho0
        rho = rho0+dr*(-(1./gamma)*rho0**(2-gamma)*mass0/r0**2)
        if(rho .le. 0.0)then
          print*,' radius = ',r0+dr,', mass = ',mass
          write(44,*)mass,r0+dr
          goto 99
        endif
         
        mass0 = mass
 
        rho0 = rho
        
        r0 = r0+dr
      enddo

99    continue

C.............Explicit taylor.......

      print*,' '
      print*,' Taylor '
      rho0 = rho00
      mass0 = rho0*4*3.14159*dr**3/3. 
      r0 = dr      
      do i = 1, n+1
        print*,r0,rho0,mass0
        write(17,*)r0,rho0
        mass = mass0+dr*3*r0**2*rho0 + 0.5*dr**2*
     & ( 6*r0*rho0 -3./gamma*(rho0)**(2-gamma)*mass0)
        rho = rho0+dr*(-(1./gamma)*rho0**(2-gamma)*mass0/r0**2)
     & +0.5*dr**2*(2/gamma*mass0*(rho0)**(2-gamma)/r0**3
     & - 3/gamma*rho0**(3-gamma)+
     & (2-gamma)/gamma**2*rho0**(3-2*gamma)*mass0**2/r0**4)
        if(rho .le. 1.0e-5*rho00)then
          print*,' radius = ',r0+dr,', mass = ',mass
          write(45,*)mass,r0+dr
          goto 199
        endif
         
        mass0 = mass
 
        rho0 = rho
        
        r0 = r0+dr
      enddo

199    continue


      end