##               COPYRIGHT NOTICE

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

##  example1.R is a utility file needed by
##  Example 1 (semiclassical quantization of molecular vibrations)
##  associated with Ch. 1 (Numerical Differentiation, Quadrature, and Roots)
##   in the series Computational Physics with Maxima or R, 
##   See  example1.pdf for more info.

##  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







##   example1.R 
##   R code for semiclassical quantization of molecular vibrations

##   An example for Ch.1 of Computational Physics with R and Maxima
##   Ted Woollett, Oct. 2013
##   http://www.csulb.edu/~woollett/


wkb = function(gam,e,n) {
	if (e <= -1.0) stop(" e must be greater than -1 ")
	if (e >= 0.0)  stop(" e must be less than zero ")	
	xmin = 2^(1/6)	
	vfun = function(x) 4*(1/x^12 - 1/x^6)	
	xin = xmin*(sqrt(e+1)/e-1/e)^(1/6)
	xout = xmin*(sqrt(e+1)+1)^(1/6)/(-e)^(1/6)
	integrand = function(x) sqrt(e - vfun(x))
	nint = integrate(integrand,xin,xout,abs.tol=1e-14,subdivisions=500L)$value
	nint - (n + 1/2)*pi/gam}
	
find.e = function(gam,ebott,etop,n) uniroot(function(e) wkb(gam,e,n),
                 c(ebott,etop),tol=1e-14)$root
                 

levels = function(gam,num) {
    bott = -0.99
    top = -5e-4
    for (nn in 0:(num-1)){
        en = find.e(gam,bott,top,nn)
        cat("  ",nn,"  ",en,"\n")
        bott = en}}
        
        
levels.info = function(gam,num) {
    bott = -0.99
    top = -5e-4
    xmin = 2^(1/6)
    x1 = function(e) xmin*(sqrt(e+1)/e-1/e)^(1/6)
    x2 = function(e) xmin*(sqrt(e+1)+1)^(1/6)/(-e)^(1/6)
    for (nn in 0:(num-1)) {
        en = find.e(gam,bott,top,nn)
        cat("  ",nn,"  ",en," ",abs(wkb(gam,en,nn))," ",x1(en)," ",x2(en),"\n")
        bott = en}}
        
        
levels.plot = function(gam,num) {
    bott = -0.99
    top = -5e-4
    plot(0:1, -1:0, type="n",xlab="",ylab="energy")
    for (nn in 0:(num-1)){
        en = find.e(gam,bott,top,nn)
        cat("  ",nn,"  ",en,"\n")
        abline(h = en,col = "blue",lwd=2)
        bott = en}}
        
        
find.xin = function(e,xstart) {
    dx = 0.001
    fdiff = function(x)  e - 4*(1/x^12 - 1/x^6)
    fstart = fdiff(xstart)
    xnew = xstart - dx
    fnew = fdiff(xnew)
    repeat {
        if (fnew*fstart < 0)  break
        xnew = xnew - dx
        fnew = fdiff(xnew)}
    uniroot(fdiff,c(xnew,xstart),tol=1e-14)$root}
    
find.xout = function(e,xstart) {
    dx = 0.001
    fdiff = function(x)  e - 4*(1/x^12 - 1/x^6)
    fstart = fdiff(xstart)
    xnew = xstart + dx
    fnew = fdiff(xnew)
    repeat {
        if (fnew*fstart < 0)  break
        xnew = xnew + dx
        fnew = fdiff(xnew)}
    uniroot(fdiff,c(xstart,xnew),tol=1e-14)$root}
    
    
wkb1 = function(gam,e,n) {	
	xmin = 2^(1/6)
	xin = find.xin(e,xmin)
	xout = find.xout(e,xmin)		
	vfun = function(x) 4*(1/x^12 - 1/x^6)	
	integrand = function(x) sqrt(e - vfun(x))
	nint = integrate(integrand,xin,xout,abs.tol=1e-14,subdivisions=500L)$value
	nint - (n + 1/2)*pi/gam}
	
find1.e = function(gam,ebott,etop,n) uniroot(function(e) wkb1(gam,e,n),
                 c(ebott,etop),tol=1e-14)$root
                 
levels1 = function(gam,num) {
    bott = -0.99
    top = -5e-4
    for (nn in 0:(num-1)){
        en = find1.e(gam,bott,top,nn)
        cat("  ",nn,"  ",en,"\n")
        bott = en}}