##               COPYRIGHT NOTICE

##  Copyright (C) 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



##   LJ6-12.R:   Numerov method 
##          for Lennard Jones 6/12 Potential
##  Ch. 3, Computational Physics with Maxima or R,
##  Boundary Value and Eigenvalue Problems


##       initial global parameters :      
  
    h = 0.01
    h2 = h^2/12    ##   Numerov constants 
    h52 = 5*h2
	gam2 = 2500   ##   square of gam = 50  
	x1decay = 0.3   ##   start yL integration at (x1 - x1decay)   
    x2decay = 1    ##   start yR integration at (x2 + x2decay)   
	y2left = 1e-19     ##  y(x_left + h)  value chosen from E = -0.9 case 
     y2right = 1e-16    ##      y(x_right - h)  value chosen from E = -0.9 case  
	 cat (" gam2 = ", gam2, "\n")
	 cat (" h = ", h, ",   x1decay = ", x1decay, ",  x2decay = ", x2decay, "\n")
	 cat (" y2left = ",y2left, "   y2right = ", y2right, "\n" )
     xm = 2^(1/6)      ##     this is where V(x) =  -1 = minimum value  
     ##    for given energy -1 < E < 0 , these are the turning points  
     xin =  function (E)    xm*(sqrt(E+1)/E-1/E)^(1/6)
     xout  =  function (E)    xm*(sqrt(E+1)+1)^(1/6)/(-E)^(1/6)
      ##    dimensionless potential V(x) for Lennard-Jones 6/12 potential 
      V  =  function (x)     4*(x^(-12) - x^(-6)) 
	  wfdebug = FALSE
	  debug = FALSE

##  wf(E)    creates ** un-normalized **  wave functions
##      for the Lennard-Jones 6/12 potential case.
##      limits of numerical integration, x_left and x_right are
## 	  determined by energy E and xdecay values.
##       The wave functions are stored in global xL, yL, xR, yR.
##	   Program also defines **global**  x2c, nleft, nright, x_left, x_right.
##	   See example run at end.
	   
##	  nleft = the number of steps from x_left to x2c =  grid point nearest to
##	  x1 = xin(E) =  classical turning point < xm = 1.122462
##	   x2 = xout(E) =  classical turning point > xm. 
##	   the global xL grid extends from x_left to x2c + h   and
##	   the global  xR grid extends from  x2c - h to   x_right,
##	   so we can compute y'(x2c) using a 3 pt. symmetric formula.

wf  =  function (E)   {

    if  ((E > 0) | (E < -1))  {
	      cat (" need -1 < E < 0   \n")
		  return(false)   }
	
    
	g  =  function(x)    gam2*( E - V(x) )     ##    coeff. func. in ode:  y''(x) + g(x) y(x) = 0  
	x1 = xin(E)      ## classical turning point for x < xm   
	x2 = xout(E)     ##  classical turning point for x> xm   
	if  (wfdebug)   cat ("  x1 = ", x1, "   x2 = ", x2, "\n")

	x_left = x1 - x1decay
	nleft = round ( (x2 - x_left)/h )    ##   number of steps from x_left to x2c = match point 
	x2c = x_left + h*nleft
	if  (wfdebug)   cat (" x_left = ",x_left," nleft = ",nleft,"  x2c = ", x2c, "\n")
	if  (wfdebug) cat (" y2left = ", y2left, " y2right = ", y2right, "\n")

	nright =  round ( (x2 + x2decay - x2c)/h )   ##   number of steps from x2c to x_right  
	x_right = x2c + h*nright
	if  (wfdebug) cat  ("  nright = ",nright," x_right = ",x_right, "\n" )

	##      find yL for x_left <=  x <= x2c + h using Numerov algorithm  
	xl = vector( mode = "numeric",  length = nleft + 2)
	yl = vector( mode = "numeric",  length = nleft + 2)
	
	xl[1] = x_left
	xl[2] = x_left + h
	yl[1] = 0
	yl[2] = y2left
	 ##  cat ( "yl = ", yl, "\n")
	
	for( j  in  2: (nleft + 1) )   {
	   x = x_left + j*h
	   xl[j+1] = x
	   yl[j+1] = ( 2*(1 - h52*g(x-h)) * yl[ j ] - (1 + h2*g(x-2*h)) * yl[ j-1] ) /  (1 + h2*g(x))   }
	   
	 ##  cat ( "yl = ", yl, "\n")
	   
	##     find yR for  x2c - h  <= x  <= x_right using Numerov method 
	
	xr = vector( mode = "numeric",  length = nright + 2)
	yr = vector( mode = "numeric",  length = nright + 2)
	
	xr[nright + 2 ] = x_right
	xr[nright + 1] =  x_right - h
	yr[nright + 2] = 0
	yr[nright + 1] = y2right
	
	for ( j  in   (nright+1):2 )  {
	   x = x2c + h*(j -3)
	   xr[ j - 1] = x
	   yr[j - 1] = ( 2*(1-h52*g(x+h)) * yr [ j ]  - (1 + h2*g(x+2*h)) * yr [j + 1] ) / (1 + h2*g(x))   }
	   
	fac = yr[2]/yl[nleft + 1]       ##     yR(x2c) / yL(x2c)    		
	yl = fac*yl
	
	##  create globally known stuff	
	x2c <<-   x2c
	nleft <<-   nleft
	nright <<-   nright
	x_left <<-   x_left
	x_right <<-   x_right
	   
	xL <<-  xl
	yL <<-  yl
	xR <<-  xr
	yR <<-  yr }
	
##  > wf(-0.9)
##  > x2c
##  [1] 1.19222
##  > nleft
##  [1] 42
##  > nright
##  [1] 100
##  > fll(xL)
##    0.772218   1.20222   44 
##  > fll(yL)
##    0   12.4282   44 
##  > fll(xR)
##    1.18222   2.19222   102 
##  > fll(yR)
##    15.5682   0   102 

##    dy_diff() uses global yL, yR,nleft, h 
##     computes numerical y'(x2c)  using
##	 symmetric three point method for
##	 both yL and yR, and returns the difference
##	 divided by y(x2c)  


dy_diff  =  function()   {
     ypL =   ( last(yL) - yL[nleft] ) / (2*h)
	ypR =   (yR[3] - yR[1]) / (2*h)
	(ypL - ypR)/ abs (yR[2])  }		
	
	
##   num_nodes()
##    count the number of nodes in yL
##     ignore region where elements of yL are
##	 tiny in magnitude.
##    takes advantage of the fact that xL elements steadily increase  

num_nodes  =  function ()    {
	x11 = x_left + x1decay
	j0 = which(xL > x11) [1]   ## position in xL where x > x11	
	n = 0
	for (j in  j0: (length(yL) - 3) )  { if (yL[ j ] * yL[ j + 1 ] < 0)   n = n + 1 }
	n  }
	
	

##  F(E)
##    energy eigenvalue if global function F(E) = 0 .
##     F(E) calls wf(E) then returns dy_diff(), but
##      returns false if E > 0.
  
 
	
F  =  function (E)    {
      if (E > 0)   {
           cat (" in F(E), E = ",E,"  should be negative   \n")
  	       return(FALSE)   }
    wf(E)
	dy_diff()  }
	
##  > F(-0.91)
##  [1] 7.19967
##  > F(-0.88)
##  [1] -24.5015
##  > uniroot(F, c(-0.91,-0.88),tol = 1e-16)$root
##  [1] -0.896404	
	
	

##   wf_plot (E) calls wf(E) and plot2d 	   
##	       creates ** un-normalized **  wave functions
##          stored in global xL, yL, xR, yR.
##           prints out number of nodes in yL
##		    prints out dy_diff .
	
wf_plot  =  function(E)  {
      wf(E)	  
	  cat ( " E = ", E, ",      ymax = ", max(yL), "\n" )
	  cat (" number of nodes = ",num_nodes(), ",      dy_diff = ",dy_diff(), "\n" )
      plot(0,type="n",xlim=c(min(xL),max(xR)),ylim=c(min(yL),max(yL)),
                                  xlab = "x",ylab = "y")
       lines(xL[1:(nleft+1)],yL[1:(nleft+1)],lwd=2,col="blue")
       lines(xR[2:(nright+2)],yR[2:(nright+2)],lwd=2,col="red")
       mygrid()	}
	   
	   
	   
##   normalize() uses the current global xL,yL, xR, yR and
##    the  utility function simp (Simpson's 1/3 rule)  to define global
##    xn and yn, with the latter being normalized.    

  
  

normalize  =  function()    {
     xn = c ( xL[1:(length(xL) - 1)], xR[3:length(xR)] )
	 yn = c ( yL[1:(length(yL) - 1)], yR[3:length(yR)] )	 
     
	 ##   we need xn to have odd number of elements to use simp 
	 if (is.even ( length (xn) ) )  {
	     xn = xn[2 : length(xn)]
		 yn = yn[2 : length(yn)]    }
		 
	 AA = simp(xn,yn^2)
	 cat ( " AA = ",AA, "\n")
	 yn = yn/sqrt(AA)	 
	 
	 x_mean = simp(xn, xn * yn^2)
     cat (" x_mean = ", x_mean, "\n")

     x2_mean = simp(xn, xn^2 * yn^2)

     delx2 = x2_mean - x_mean^2     ##   this should be positive!  

     delx = sqrt(delx2)
	 
	 cat (" delx = ", delx, "\n")

      xn  <<-   xn
	  yn  <<-   yn	 }
	 
##  for ground state:	 
##  > normalize()
##   AA =  85.0069 
##   x_mean =  1.14069 
##   delx =  0.0455039 
##  > fll(xn)
##    0.781455   2.20145   143 
##  > fll(yn)
##    -3.14183e-20   0   143 
##  > max(yn)
##  [1] 2.97425


##   yn_plot_current() uses the currently defined normalized set (xn,yn) 

	 
yn_plot_current  =  function ()   {
	 ymn = floor( min(yn) )
	 ymx =  1 + floor( max (yn) ) 
	 cat (" ymax = ", max(yn), "\n" )
	 plot(xn,yn,type = "l",ylim = c(ymn,ymx),lwd=3,col="blue",xlab="x",ylab="y")
	 mygrid() }
	 
##   yn_plot(E,xmin,xmax) first calls wf(E) to create 
##     un-normalized wave functions corresponding to the
##	 given energy E. Then normalizes those wave functions
##	 to produce the vectors xn and yn. Finally makes a plot
##	 of yn over only the region (xmn, xmx)   */
					  
yn_plot =  function (E,xmn,xmx)   {
     wf(E)
	 cat (" E = ",E, "\n" )
	 cat (" number of nodes = ",num_nodes(),",       dy_diff = ",dy_diff(), "\n" )
	 normalize()
	 cat (" normalized ymax = ", max(yn), "\n" )
	 ymn = floor( min(yn) )
	 ymx =  1 + floor( max (yn) ) 
	 plot(xn,yn,type = "l",ylim = c(ymn,ymx), xlim = c(xmn, xmx),
	                lwd=3,col="blue",xlab="x",ylab="y")
	 mygrid() }
	 
				  
				  

##  bracket is a modified version of bracket_basic, designed to work with
##    the function F(E)  which can return FALSE.
	
##    bracket 	looks for a sign change in func,
##	starting with xx,  and increasing xx by dxx each step.
##	If sign change is found,  then we back up to the previous xx 
##	and search with new  dxx value one half of the previous value.
##   normally returns [ea,eb], but if can't find change in sign,
##   then returns [0,0], and if func returns FALSE, then 
##    bracket returns FALSE.

bracket =  function (func,xx,dxx,xacc)  {
     x = xx
	 dx = dxx
	 it = 0 
	 itmax = 1000
     anerror = FALSE
	 anerror2 = FALSE

    repeat  {
    	  it = it + 1  
	      if (it > itmax)  {
              cat (" can't find change in sign  \n")
		      anerror = TRUE
		      break}
	      
	       x1 = x
	      x2 = x + dx
	      if ( debug )  cat (" it = ",it," x1 = ",x1," x2 = ",x2," dx = ", dx, "\n")
	  
	      f1 = func(x1)
	      if ( f1 == FALSE)   {
	          cat (" in bracket, f1 = FALSE , x1 = ",x1, "  dx = ", dx, " \n ")
		      anerror2 = TRUE
		      break }
		  
	      f2 = func(x2)
	      if ( f2 == FALSE)   {
	          cat (" in bracket, f2 = FALSE , x2 = ",x2, "  dx = ", dx, " \n ")
		      anerror2 = TRUE
		      break }
		  
	  
	      if  ( f1 * f2  < 0 )    {
	           if  ( abs(dx) < xacc )   break
		       x = x - dx
		       dx = dx/2 }  else x = x2 }
	  
    if (anerror) c(0,0) else if (anerror2) FALSE  else c(x1,x2) 	  }
	
	
	
##  levels(Emin,Emax,dE, Eacc ) returns a vector  of energy levels with
##     increasingly larger number of nodes in energy range (Emin, Emax).
##	 uses F(E) to find  roots, and calls wf, num_nodes() and dy_diff() for each root found,
##    Uses bracket and uniroot.
##	The arguments (dE, Eacc) are used to call bracket, and do not describe
##	  the accuracy of the energy levels found.
##	Once a good energy e.v. is found we look for the region of
##	energies with one more node and search there.
##	interactive continue or stop.

levels  =  function (Emin,Emax,dE, Eacc )   {     
     rmax = 20
	 eL = rep(NA, rmax)	 ##  vector  eL will hold energy eigenvalues found  
     e = Emin
	 nlast =  -1
	 j = 1
	 
	 repeat {
         	 if  (e > Emax  |  j > rmax)  break    ## exit do loop  
	         cat ("---------------- levels -------------------\n")
	        cat (" nlast = ", nlast,"\n")
	        cat (" Estart = ", e," dE = ", dE, "\n" )
            out =  bracket(F,e, dE, Eacc)
			cat (" ea = ",out[1]," eb = ",out[2],"\n")
			if (out[1] == 0)    {
			    cat (" can't find bracket interval \n")
				cat (" e = ",e, "\n")				
				break  }
			eroot = uniroot(F, out, tol = 1e-16)$root
            cat (" eroot = ", eroot, "\n")			
			wf(eroot)
			nn = num_nodes()
			cat (" number of nodes = ", nn, "\n")
			cat (" dy_diff at x = x2c is    ", dy_diff(), "\n" )
			eL[ j ]  = eroot
            nlast = nn	
			j = j + 1
			r = readline (" input c or s \n ")
			if (r == "s")  break	    ## exit do loop  

			##    search for an  e value greater than eb which produces
			##      a wave function with nn + 1 nodes 
			enext  =  out[2] + dE
			repeat  {
			        wf(enext)
					if  (num_nodes() > nlast)   {
					     e = enext
						 break }  else enext =  enext + dE }  }   ## end of outer repeat loop
		##   remove na's at end of eL				 
	 eL[!is.na(eL)]  }