PRO test
COMMON LANE, n
COMMON LANE_DATA, data, ndata, h, zmax
COMMON LANE_INTERPOL, a, b, c, d

IF N_ELEMENTS(zmax) EQ 0 THEN lane_emden, 3.D0

msun = 1.989D33         ;mass of the sun in grams
mx = 1000.D0            ;mass of new star, in solar masses
mstar = mx*msun         ;calculates mass of star for polytrope
rho_c = 1.D0            ;central density in g/cm^3
radius = ((-mstar*zmax)/(4.D0*!dpi*rho_c*(theta(zmax))[1]))^(1.D0/3.D0)         ;calculates radius from polytrope, mass, and central density

ngrid = 500
x = dindgen(ngrid+1)/ngrid*zmax
z = theta(x)

rfac=(-mstar/(4.D0*!dpi*rho_c*zmax^2*theta(zmax))[1])^(1.D0/3.D0)
mfac=-4.D0*!dpi*((radius/zmax)^3)*rho_c

max_change = 0.1D0              ;maximum fractional change in density between adjacent zones
cutoff = 1.D-20                 ;defines density of outmost zone of the star (approaches zero)

iter=0

REPEAT BEGIN
        m = x^2*z[*,1]*mfac             ;creates physical mass profile from polytrope
        r = x*rfac              ;creates radial profile

;test prints
;print, 'x is',n_elements(x)
;print, 'theta(zmax)',n_elements(theta(zmax))

        array_end = N_ELEMENTS(r)-1
        r1 = r[1:array_end]
        m1 = m[1:array_end]
        r0 = r[0:array_end-1]
        m0 = m[0:array_end-1]
        dV = 4.D0*!dpi*(r1-r0)*(r1^2+r1*r0+r0^2)/3.D0
        dM = m1 - m0

        mrat=2.D0*abs(dM)/(m1+m0)

        ii=WHERE(mrat LT 1.d-7)

        rho = dM/dV             ;density profile
        
; the following is very crude and should be refined using a 3rd order
; polynomial fit to density and do an analytic integral over the zone
        rho[ii]=rho_c*((z[ii,0]+z[ii+1,0])*0.5D0)^n

        rho_end = N_ELEMENTS(rho)-1
        rho1 = rho[1:rho_end]
        rho0 = rho[0:rho_end-1]
        delta_rho = (rho1 - rho0)/rho0          ;determines change in density in adjacent zones
        out_range = WHERE((ABS(delta_rho) GT max_change) and (rho[0:rho_end-1] GT cutoff))      ;determines where to create new zones

        IF out_range[0] NE -1 THEN BEGIN                ;creates new zones
           ii=[out_range,out_range+1]
           ii=ii[UNIQ(ii,SORT(ii))]
           new_x = [x[ii] + x[ii+1]]*0.5D0

;           new_x = [x[out_range +1] + x[out_range +2],x[out_range]+x[out_range+1]]*0.5D0
           x = [x, new_x]
           k = SORT(x)
           z = [z, theta(new_x)]
           z = z[k, *]
           x = x[k]
           PRINT,'iteration=',iter+1,', zones=',N_ELEMENTS(x)-1
        ENDIF
        iter+=1
ENDREP UNTIL (out_range[0] EQ -1) OR N_ELEMENTS(r) GT 1000 
plot, r, rho, /ylog
oplot, r, rho, psym=1

; finalize mass coordinate
nr=N_ELEMENTS(r)

xm=rho*4.D0*!dpi*(r1-r0)*(r1^2+r1*r0+r0^2)/3.D0
ym=[REVERSE(TOTAL(REVERSE(xm),/CUMULATIVE,/DOUBLE)),0]
zm=ym[0]-ym

plot,alog10(ym/msun), rho, /yl

END
FUNCTION differential, x, y
COMMON LANE, n


;function to calculate differential used in RK4
diff = [y[1], 0.D0]
IF x NE 0 THEN diff[1] = -2.D0 / x * y[1] - y[0]^n
RETURN, diff
END

;-------------------------------------------------

FUNCTION theta, z, n=n0

COMMON LANE, n
COMMON LANE_DATA, data, ndata, h, zmax
COMMON LANE_INTERPOL, a, b, c, d

IF N_ELEMENTS(n0) EQ 0 THEN n0=0

IF (N_ELEMENTS(n0) EQ 0) THEN LANE_EMDEN
IF (n0 NE 0) THEN BEGIN
        IF n0 GT 1 THEN BEGIN
                LANE_EMDEN, n0
                n0 = 0 
        ENDIF
ENDIF

zh = (z<zmax)/h
ii = FLOOR(zh)
f = zh-ii
;f = [[f],[f]]
;RETURN, data[ii+1,0:1]*f+data[ii,0:1]*(1.D0-f)
RETURN, [[a[ii]+b[ii]*f+c[ii]*f^2+d[ii]*f^3],[(b[ii]+2.D0*c[ii]*f+3.D0*d[ii]*f^2)/h]]

END

;------------------------------------------------

PRO LANE_EMDEN, n0
COMMON LANE, n
COMMON LANE_DATA, data, ndata, h, zmax
COMMON LANE_INTERPOL, a, b, c, d 

IF N_ELEMENTS(n0) EQ 0 THEN n0 = 3.D0
n = n0

h = 0.5D0

acc = 1.D-9
 
yy=0

;fh = (sqrt(5.D0) - 1.D0) * 0.5D0 
fh = 0.5D0

REPEAT BEGIN
        h*=fh
        data = DBLARR(ROUND(10.D0/h),2)
        ndata = 0L
        x = 0.D0
        y = [1.D0, 0.D0]
        data[ndata, 0:1]=y[0:1]
        WHILE y[0] GT 0.D0 DO BEGIN
                dydx = differential(X, Y)
                result = RK4(y, dydx, x, h, 'differential',/DOUBLE)
                x+= h
                y0 = y
                y = result
                ndata+=1
                data[ndata, 0:1] = y[0:1]
        ENDWHILE
        yy0 = yy
        yy=x+h*y[0]/(y0[0]-y[0])
        PRINT, x, y[1]
ENDREP UNTIL (ABS(yy-yy0) LT acc)

;calculate the coefficients for a third order interpolation
data1 = data[1:ndata, *]
data0 = data[0:ndata-1, *]

a = data0[*, 0]
b = data0[*, 1]*h
c = 3.D0*data1[*, 0] - data1[*, 1]*h - 3.D0*data0[*, 0] - 2.D0*data0[*, 1]*h
d = data1[*, 1]*h - 2.D0*data1[*, 0] + 2.D0*data0[*, 0] + data0[*, 1]*h

data=data[0:ndata,*]
zmax = yy

END

;------------------------------------------------