program lane_emden

c**************************************************************
c* program to solve the Lane-Emden equation numerically       *
c*                                                            *
c* the equation to solve is                                   *
c* d^2x/dt^2 = - ((2/t)*(dx/dt) + x**n)                       *
c*                                                            *
c* where x is the dimensionless density function (psi)        *
c* and t is the dimensionless length (xi)                     *
c* initial conditions are      x = 1                          *
c*                         dx/dt = 0                          *
c*                                                            *
c* and the power index is n.                                  *
c* analytical solutions are known for n=0, 1, 5.              *
c*                                                            *
c* CF. Oct 1999.                                              *
c**************************************************************

	real x, dxdt, n, dt, xlast, dxdtlast
	character*50 name

* index n of the polytrope
	
	write(6,*) 'Enter n, polytrope index '
	read(5,*) n
	
	write(6,*) 'Enter name of results file'
	read(5,*) name

* open the file for output of results

	open(20,file=name,status='unknown')

* variables are 
* dimensionless length : t 
* dimensionless density function : x

* specify initial values 

	x    = 1.0   
	dxdt = 0.0
        t    = 1.0E-6

* "radius step" is in length variable, t
* the smaller the step the more accurate the calculation
* and the longer it takes!

	dt = 0.001

* now the MAIN program loop

	do 100 i = 1, 10000000

* calculate the next dxdt  		
	  dxdt = dxdt - ( 2.0*dxdt/t + x**n )*dt
	
* calculate the next x    
	  x = x + dxdt*dt 

* increment the dimensionless radius (move outward in the star)
	  t = t + dt
	
* write out the results

	  write(20,*) t,x,dxdt

* check if the solution has become negative yet  
	  if (x.lt.0.0) goto 999

* save the current values of x and t, for determining the root later
	  xlast = x
	  dxdtlast = dxdt

100	continue

999	continue

* calculate the position of the root, using the current x (which 
* is the first value < 0, and the previous x (stored as xlast) 
* which is the last calculated value which was still > 0

	t1 = t + (x/(xlast-x))*dt
	dphidt1 = dxdt - (t-t1)*(dxdt-dxdtlast)/dt

	write(6,*) 'first root of phi at xi=',t1 
	write(6,*) 'at this point dphi/dt = ',dphidt1

	write(6,*) n,t1,-t1**2*dphidt1

	end