##               COPYRIGHT NOTICE

##  Copyright (C) 2014, 2015  Edwin L. (Ted) Woollett
##  http://www.csulb.edu/~woollett

##  This program is free software; you can redistribute
##  it and/or modify it under the terms of the
##  GNU General Public License as published by
##  the Free Software Foundation; either version 2 
##  of the License, or (at your option) any later version. 

##  This program is distributed in the hope that it
##  will be useful, but WITHOUT ANY WARRANTY;
##  without even the implied warranty of MERCHANTABILITY
##  or FITNESS FOR A PARTICULAR PURPOSE. See the 
##  GNU General Public License for more details at
##  http://www.gnu.org/copyleft/gpl.html




##  project2.R  July, 2014
##  Ted Woollett:     http://www.csulb.edu/~woollett
## project2: Structure of a White Dwarf Star
##  Ch. 2, Computational Physics with Maxima or R,
##   Initial Value Problems

##   rk4_step(init,t,dt,func) returns y(t+dt) if one dependent variable
##    or returns the vector c(x(t+dt), vx(t+dt)) if simple harmonic oscillator (two dependent
##     variables), etc.

rk4_step = function(init, t, dt, func)  {
   num.var = length(init)   
   h = dt   # step size 
   yL = init   
   k1 = func(t, yL)  #  vector of derivatives
   k2 = func(t + h/2, yL + h*k1/2)  #  vector of derivatives
   k3 = func(t + h/2, yL + h*k2/2)  #  vector of derivatives
   k4 = func(t + h, yL + h*k3)  #  vector of derivatives
   for (k in 1:num.var) yL[k] = yL[k] + h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 
   if (num.var==1) yL[1] else yL}
   
##  example 1: one dependent variable, integrate dy/dt = -y, one step dt = 0.1
##  from t=0, y(0)=1, using same derivatives function deriv as used with euler1 above.
##  > rkstep = rk4_step(1, 0, 0.1, deriv)
##  > rkstep
##  [1] 0.9048375

## example2: two dependent variables, simple harmonic oscillator with unit period
##  same derivatives function sho as used above with x(0)=1, vx(0) = 0.
##  > rkstep = rk4_step(c(1,0), 0, 0.1, sho)
##  > rkstep
##  [1]  0.8091019 -3.6880842

##  derivatives function 'derivs' for dm/dr = rho*r^2, drho/dr = -m*rho/gamma/r^2
##  y[1] = m, y[2] = rho
##  gamma = rho^(2/3)/3/sqrt(1+rho^(2/3))

derivs = function(r,y) {
    with( as.list(y), { 
         y23 = y[2]^(2/3)
         gamma = y23/3/sqrt(1+y23)
         dm = r^2*y[2]
         drho = -y[1]*y[2]/gamma/r^2
         c(dm, drho)})}
		 
##  R returns rho_new = NaN from rkstep if rho goes negative
##  > rkstep = rk4_step(c(0.70582, 0.0036087), 2.5, 0.1, derivs)
##  > rkstep
##  [1] 0.7067797       NaN
##  > is.nan(rkstep[2])
##  [1] TRUE


##    dwarf1  --> list(rL, mL, rhoL)
##     With rbar -> r, rhobar -> rho,  mbar -> m, rhobar_c -> rhoc, drbar -> dr
##     dwarf1(rhoc,dr,rtol,nmax) returns a R list with the vector elements rL,mL,and rhoL,
##     given the central density rhoc at r = 0, the step size dr, and the rtol value and the maximum
##       number of elements per vector, nmax,
##     (using the externally defined derivs function when  calling rk4_step)
##     after taylor expanding the dependent variables out to r1=small,  and then watching
##     for rho to pass from positive to negative value.
##            If rho is found to be negative, we back up and let rtol = smaller be the
##     step size, integrating forward again until rho is found to be negative, then
##     taking the last good values as the final values of r,m and rho   
##     and using that last good value of r to define the radius R and that 
##     last good value of m(r) to define M.
##    Inside dwarf1,  rk4_step returns the vector c(m,rho).

dwarf1 = function(rhoc, dr, rtol, nmax)  {

##  we don't know how long these vectors need to be		

		rL = vector(length = nmax)   # radius values
		mL = vector(length = nmax)     # mass values
		rhoL = vector(length = nmax)     #  mass density values
		
		rL[1] = 0
		mL[1] = 0
		rhoL[1] = rhoc

##  	    taylor expand away from r = 0 to avoid division by zero problems
		rh23 = rhoc^(2/3)
		gamc = rh23/3/sqrt(1+rh23) 
		r1 = 1e-12
		m1 =  rhoc*r1^3/3
		rho1 =  rhoc - r1^2*rhoc^2/gamc/6
		 
		rL[2] = r1
		mL[2] = m1
		rhoL[2] = rho1		
		   
##         first loop: search for rho = 0 position using step dr 
        num  = 2     # number of vector elements defined so far
##		cat(" nmax = ",nmax," num = ",num," \n")
##		cat(" loop 1 \n")
		
		repeat {
		    rkstep = rk4_step( c(m1, rho1), r1, dr, derivs)
##            watch for rho becoming negative 
			if (is.nan(rkstep[2])) {
##			    cat(" with r1 = ",r1," m1 = ",m1," rho1 = ",rho1," \n")
##			    cat(" new rho is negative-- go back and try again with smaller step \n")				
				break }
			
			if (num == nmax) {
##			    cat(" num = ",num," = nmax; abort integration  \n")
				return() }
				
			r1 = r1 + dr
			m1 = rkstep[1]
			rho1 = rkstep[2]
			num = num + 1
##			cat(" num = ",num, " \n")
			rL[num] = r1
			mL[num] = m1
			rhoL[num] = rho1} 
				
##      second  loop: search for rho=0 location using smaller step size  rtol
##         cat(" start loop 2 with num = ",num," \n")
##		 cat(" start loop 2 with with r1 = ",r1," m1 = ",m1," rho1 = ",rho1," \n")
		 repeat {
		    rkstep = rk4_step( c(m1, rho1), r1, rtol, derivs)
##            watch for rho becoming negative 
			if (is.nan(rkstep[2])) {
##			    cat(" rho is negative -- use last good value as surface \n")
				break }
			
			if (num == nmax) {
##			    cat(" num = nmax; abort integration  \n")
				return() }
				
			r1 = r1 + rtol
##			cat(" r1 = ",r1, " \n")
			m1 = rkstep[1]
			rho1 = rkstep[2]
			num = num + 1
##			cat(" num = ",num," \n")
			rL[num] = r1
			mL[num] = m1
			rhoL[num] = rho1} 
		 
##        construct vectors of length num for return  
##		 cat(" final num = ",num," \n")
		 
		 rrL = vector (length = num)
		 rmL = vector (length = num)
		 rrhoL = vector (length = num)
		 
		 for (kk in 1:num) {
		      rrL[kk]  = rL[kk]
			  rmL[kk] = mL[kk]
			  rrhoL[kk] = rhoL[kk] }
			  
		 list(rrL, rmL,  rrhoL) }
		 
		 
##   dwarf3(rhoc, dr, rtol) calls dwarf1 and returns c(Mbar, Rbar) 

dwarf3  =  function(rhoc, dr, rtol)  {
     nmax = 500
     soln = dwarf1(rhoc, dr, rtol, nmax)   # list(rL,mL,rhoL)
	 surfL = get_last(soln)
	 c(surfL[2], surfL[1]) }
	 
##  > dwarf3(1,0.1,0.01)
##  [1] 0.7066441 2.4900000
		 
##   dwarf4() calls dwarf3 to accumulate  vectors of values of M and R
##     returns list(ML, RL)  

dwarf4 = function(rho_vals)  {
       drv = 0.1
	   rtolv = 0.01
	   num = length(rho_vals)
	   mv = vector(length=num)
	   rv = vector(length=num)
	   k = 1
	   for (rhoc in rho_vals) {
	        rd3 = dwarf3(rhoc,drv,rtolv)
			mv[k] = rd3[1]
			rv[k] = rd3[2]
			k = k + 1}
	   list(mv, rv) }

	   
##  Trapezoidal Rule: allows non-uniform grid, from Ch. 1 file cp1code.R

trapz = function(xv,yv) {
   idx = 2:length(xv)
   as.numeric( (xv[idx] - xv[idx - 1]) %*% (yv[idx - 1] + yv[idx])/2)}




##  fll is a home-made vector utility: prints  first, last, and length of a vector
   
fll = function(xL) {
    xlen = length(xL)
    cat(" ",xL[1]," ",xL[xlen]," ",xlen,"\n") }
	
##  > fll(xL)
##    1   0.9959199   11 

mygrid = function() grid(lty = "solid", col = "darkgray")

last = function(xL) { xL[ length(xL)] }

##  get_last(list.of.vecs) returns a vector whose elements
##  are the last components of the vectors    
    
get_last = function(rkvecs) sapply(rkvecs, last)

##  > myL = list(c(1,2,3,4),c(5,6,7,8),c(9,10,11,12))
##  > get_last(myL)
##  [1]  4  8 12