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