/*               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
*/





/* FW.mac uses Runge-Kutta
   for finite well.
   Ch.3 Computational Physics with Maxima or R:
       Boundary Value and Eigenvalue Problems
   
   dimensionless units
   V = 1 for x < 0 and x > 1
   V = 0 for 0 < x < 1
   y''(x) + gam2*(E - V(x))*y(x) = 0
   gam2 = gam^2 = 2500
   gam = 50 = sqrt(2*m*L^2*V0/ hbar^2)
   
*/
  
 /*     initial global parameters:      */
  
 (  N : 100,
    h : 0.01,
    gam : 50,
	gam2 : gam^2, 
	xdecay : 0.5,  /* start yL1 integration at x = -xdecay  */
	   /* start yR integration at x = 1 + xdecay */
	ypleft : 1e-8,   
	ypright : 1e-8,
	 print("  gam = ",gam, "    gam2 = ", gam2),
	 print(" h = ", h, ",   xdecay = ", xdecay, ",  ypleft = ", ypleft,"  ypright = ", ypright)   )$
	 
/*  analytic method functions  */ 		

Fa(k) := (float(k - (2*k^2 - gam2)*tan(k)/2/sqrt(gam2 - k^2)) )$

Frhs(k) := (float( (2*k^2 - gam2)*tan(k)/2/sqrt(gam2 - k^2)) )$
				
kroot_plot(kkmin,kkmax, ymn, ymx) := 
	 ( plot2d([kk, Frhs(kk)], [kk, kkmin, kkmax],
	          [y, ymn, ymx],  [style, [lines,2]], [xlabel, "k"],
               [legend, false], [ylabel, ""],
                [gnuplot_preamble," set grid"]))$

kroot(kmin,kmax) := (find_root(Fa, kmin, kmax))$

ktoE (kv) := (kv^2/gam2)$

/* bracket_basic looks for a sign change in func, starting with x,
   and increasing x by dx each step. If sign change is found,
   then we back up to the previous x and search with new
   dx value one half of the previous value.
   Used by levels_analytic(kmax).
*/

bracket_basic(func,xx,dxx,xacc) :=
block([ x:xx, dx:dxx,x1,x2,it:0,itmax:1000],        	
    do (
	  it : it + 1,	  
	  if it > itmax then (
          print(" can't find change in sign "),    
	      return([0, 0 ])),
	  x1 : x,
	  x2 : x + dx,
	  if debug then print(" it = ",it," x1 = ",x1," x2 = ",x2," dx = ", dx),
	  
	  if func(x1)*func(x2) < 0 then (
	       if abs(dx) < xacc then return([x1,x2]),
		   x : x - dx,
		   dx : dx/2)
	  else x : x2))$
	  
/*
example using func = sin:

(%i44) [xa,xb] : bracket_basic(sin,3,0.01,0.001);
(%o44) [3.14125,3.141875]
(%i45) xv : find_root(sin,xa,xb);
(%o45) 3.1415927


zero node ground state using func = Fa:

(%i4) [ka,kb] : bracket_basic(Fa,1.6,0.1,0.05);
(%o4) [3.0,3.025]
(%i5) kv : find_root(Fa,ka,kb);
(%o5) 3.0206914

*/


	 
/* analytic_wf(E)   creates global functions  yn1(x), yn2(x), yn0(x), yna(x) ,
       and computes and prints analytic values of delx, delp.
     yn1(x) is the wave function for x < 0
	 yn2(x) is the wave function for x > 1
	 yn0(x) is the wave function for 0 <= x <= 1
	 yna(x) is the wave function for all x, usable by plot2d,
	     but not by integrate.	 

	   */

analytic_wf(E) :=
block([D,sd,skd,i1,i2,iR,x_mean,x2_mean,
             delx,delx2,ydy_sum,delp,delp2,
			k,  k1,k2,  delta,fac1,A1,A2,numer],numer:true,
     print("  E  =  ", E),
	 k : gam*sqrt(E),
     k1 : gam*sqrt( 1 - E ),	 
	 k2 : k1,
	 delta : atan(k/k1),	 
	 sd : sin(delta),	 
	 skd : sin(k + delta),	 
	 
	 D : sd^2/2/k1 + skd^2/2/k2 - sin(2*k+2*delta)/4/k + sin(2*delta)/4/k + 1/2, 	
	fac1 : 1/sqrt(D),
	A1 : fac1*sd,
	A2 : fac1*skd,
	
		/* calc < x > */
		   
	i1 : - A1^2/4/k1^2,
	
	iR : A2^2 * (1 + 2*k2)/4/k2^2,
	
	i2 : fac1^2*( -sin(2*k+2*delta)/(4*k)-cos(2*k+2*delta)/(8*k^2)
	                +cos(2*delta)/(8*k^2)+1/4 ),	
					
	x_mean : i1 + i2 + iR,
	print(" x_mean = ", x_mean),			
	
	/*  calc < x^2 >  */
		
	i1 :  A1^2/4/k1^3,
	
	iR :  A2^2 * (1 + 2*k2 + 2*k2^2)/4/k2^3,
	
	i2 : fac1^2 * (  -sin(2*k+2*delta)/(4*k)+sin(2*k+2*delta)/(8*k^3)-cos(2*k+2*delta)/(4*k^2)
                          -sin(2*delta)/(8*k^3)+1/6  ),
						  
	x2_mean : i1 + i2 + iR,	
		
	delx2 : x2_mean - x_mean^2,   /* this should be positive!  */

     delx : sqrt(delx2),
	 
	 print(" delx = ", delx),	

     /*  calc dyd_sum  */	 
	 
	 i1 : A1^2 / 2,
	 
	 iR : -A2^2 / 2,
	 
	 i2 : fac1^2 * ( cos(delta)^2 - cos(k + delta)^2 ) / 2,
	 
	 ydy_sum : i1 + i2 + iR,
	 
	 print(" ydy_sum = ", ydy_sum),	 
	 
	 /* calc <p^2>  */
	 
	 i1 : -A1^2 * k1/2,
	 
	 iR : - A2^2 * k2/2,
	 
	 i2 : (fac1*k)^2 * (  - sin(2*k+2*delta)/4/k + sin(2*delta)/4/k + 1/2  ),
	 
	 delp2 : i1 + i2 + iR,
			   
	delp : sqrt(delp2),
	
	print(" delp = ",delp,"  delx*delp = ", delx*delp),					

	/* use 'define' to create globally usable functions
	    yn1(x), yn2(x), yn0(x), yna(x)  */
		
	define(yn1(x), ev(A1*exp(k1*x))),			

	define(yn2(x), ev(A2*exp(k2)*exp(-k2*x))),	

	define(yn0(x), ev(fac1*sin(k*x + delta))),
	
	yna(x) := 
          ( if x < 0 then yn1(x)
             else if x > 1 then yn2(x)
             else  yn0(x) ),	

	done)$
	
	 
/*
11 node solution:

(%i6) kroot_plot(36,36.1);
(%o6) ""
(%i7) k0 : kroot(36,36.1);
(%o7) 36.086514
(%i8) ktoE (kv) := (kv^2/gam2)$
(%i9) E0 : ktoE(k0);
(%o9) 0.520895
(%i18) plot2d(yn0(x),[x,0,1],[ylabel,"y0(x)"],
           [style,[lines,2]],[gnuplot_preamble,"set grid"])$

*/


/* nodes_analytic(dx) counts the number of nodes implied by
   the function yn0(x) created by analytic_wf(E)  */
   
   
nodes_analytic ( dx ) :=
block( [num:0, xv:0, xnew, f1, f2, numer], numer:true,
    do (
         f1 : yn0(xv),
         xnew : xv + dx,
	     if xnew > 1  then return(),
         f2 :  yn0(xv + dx),
         if f1*f2 < 0 then num : num + 1,
         xv : xnew),
	num)$
	
/*  
  analytic
  zero nodes  ground state  solution  
  
(%i23) kroot_plot(2.8,3.2);
(%o23) ""
(%i24) kval : kroot(2.9,3.1);
(%o24) 3.0206914
(%i25) Eval : ktoE(kval);
(%o25) 0.00364983
(%i26) analytic_wf(Eval);
 x_mean =  0.5 
 delx =  0.18802 
 ydy_sum =  -1.78676518E-16 
 delp =  2.9619274   delx*delp =  0.556902 
(%o26) done

(%i29) nodes_analytic(0.01);
(%o29) 0

(%i30) plot2d(yn0(x),[x,0,1],[ylabel,"y0(x)"],
           [style,[lines,2]],[gnuplot_preamble,"set grid"])$
		   
		   shows zero node ground state solution over the
		   interval [x,0,1]
		   
plot_analytic (Eval)  calls analytic_wf(E)
		   
(%i40) plot_analytic(Eval);
 x_mean =  0.5 
 delx =  0.18802 
 ydy_sum =  -1.78676518E-16 
 delp =  2.9619274   delx*delp =  0.556902 
 number of nodes =  0 
(%o40) ""
		   
*/	

/* plot_analytic(E) calls analytic_wf(E), prints out 
    the number of nodes, and makes a plot of yna(x)
	*/
	
plot_analytic(E) :=
block ( [ddx : 0.001, xvL, ynL, ymn, ymx,
              xmn: -0.25, xmx : 1.25,  numer], numer:true,
	 analytic_wf(E),
     xvL : makelist(x,x,0,1,ddx),
     ynL : map(yn0, xvL),     
	 ymn : floor( lmin(ynL) ),
	 ymx : 1 + floor( lmax (ynL) ),     
	 print(" number of nodes = ", nodes_analytic(ddx) ),
	plot2d(yna(x), [x,xmn,xmx], [ylabel, "yna(x)"],[y,ymn,ymx],
           [style,[lines,2]],[gnuplot_preamble,"set grid"]))$
		   
/*  plot_ps_help() creates global ymn, ymx assoc
     with current **analytic** yn0  */
		   
plot_ps_help() :=
block([xvL,ynL,numer],numer:true,
     xvL : makelist(x,x,0,1,0.001),
     ynL : map(yn0, xvL),     
	 ymn : floor( lmin(ynL) ),
	 ymx : 1 + floor( lmax (ynL) ))$
	 
/*
(%i19) plot_ps_help()$
(%i20) ymn;
(%o20) 0
(%i21) ymx;
(%o21) 2
(%i22) plot2d(yna(x), [x,-0.5,1.5], [ylabel, "yna(x)"],[y,ymn,ymx],
           [style,[lines,5]],[gnuplot_preamble,"set grid"],
           [psfile,"c:/k3/fw3.eps"]);
(%o22) "c:/k3/fw3.eps"
*/	 
	
	
/*
  using ground state energy Eval from above, 
  if we call plot_analytic(Eval + 0.01) we get
  a one node solution which is not smooth:
  
(%i44) plot_analytic(Eval+0.01);
 x_mean =  0.516877 
 delx =  0.285471 
 ydy_sum =  1.66533454E-16 
 delp =  5.3591208   delx*delp =  1.5298736 
 number of nodes =  1 
(%o44) ""
*/

/* levels_analytic(kmax) returns a list of analytic eigen energies
   related to k via:  k = gam*sqrt(E). The corresponding eigenvalues
   of k are separated by roughly dk = 3 and lie in the middle of the intervals
   [n1*%pi/2, n2*%pi/2], where n1 and n2 are adjacent odd integers, n2 > n1.
   The ground state (zero node soln) corresponds to k = 3.02  (n1=1,n2=3).
*/


levels_analytic(kmax) :=
block( [ka,kb,kv,level:0, nmax, Elist : [],  numer], numer:true,
     nmax : ceiling(2*kmax/%pi),
	 print("    nodes      E    "),
	 print( "  "),
	 /* make nmax an odd integer */
	 if evenp (nmax) then nmax : nmax + 1,
	 /* analytic kv using n*pi/2 + 0.1 with n odd  */
	 for j:1 step 2  thru nmax do (
	    [ka, kb] : bracket_basic( Fa,   j*%pi/2 + 0.1,  0.01,0.005),
		kv : find_root( Fa, ka, kb),
		Ev : ktoE(kv),
		Elist : cons(Ev, Elist),
		print( "   ", level,  "   ", Ev ),
		level : level + 1),
	 reverse(Elist) )$
	 
/* 
(%i1) load(cp3);
(%o1) "c:/k3/cp3.mac"
(%i2) load(FW);
  gam =  50     gam2 =  2500 
 h =  0.01 ,   xdecay =  0.5 ,  ypleft =  1.0E-8   ypright =  1.0E-8 
(%o2) "c:/k3/FW.mac"
(%i3) levels_analytic(20);
    nodes      E     
   
    0     0.00364983 
    1     0.0145973 
    2     0.032836 
    3     0.0583552 
    4     0.0911389 
    5     0.131165 
    6     0.178405 
(%o3) [0.00364983,0.0145973,0.032836,0.0583552,0.0911389,0.131165,0.178405]
*/





	
	 
/***************  functions for numerical integration  **************************/ 
	 
/*  wf(E)    creates ** un-normalized **  numerical wave functions
         using Runge-Kutta routine rk4.       
       The wave functions are stored in global xL1, yL1,ypL1, xL2, yL2, ypL2, xR, yR, ypR .
	   Program also defines **global**   nleft, nright, ncenter	  	   	  
	   the global xL1 grid extends from -xdecay to 0   and
	   the global xL2 grid extends from 0 to 1 and
	   the global  xR grid extends from  1 to 1 + xdecay
	   
*/



wf(E)  :=
block( 	[  glr,gc, outL,  fac,  numer], numer : true,
        
     if (E < 0) or (E > 1) then (
	      print(" need 0 < E < 1 "),
	      return(false)),
	 ncenter : N,
	 if not integerp(ncenter) then (
	     print (" ncenter = ",ncenter,"  is not an integer"),
		 return(false)),
	 nleft : round(xdecay/h),
	 nright : nleft,
	 if wfdebug then print(" nleft = ",nleft," ncenter = ",ncenter," nright = ",nright),

	 glr : gam2*(E - 1),          /* g(x) for x < 0 and x > 1   */
	 gc : gam2*E,                 /* g(x) for 0 < x < 1  */
	 if wfdebug then  print(" glr = ",glr,"   gc = ", gc),
	 
	 
	 
	 /* construct xL1, yL1, and ypL1 for -xdecay < x < 0  */
	 
     outL : rk4(['y2, - glr* 'y1], ['y1, 'y2], [ 0, ypleft], ['x, -xdecay, 0, h] ),
	 xL1 : take(outL,1),
	 yL1 : take(outL,2),
	 ypL1 : take(outL,3),	 
	 
	 /*  construct xL2, yL2, and ypL2 for 0 < x < 1  */

	 outL : rk4(['y2, - gc* 'y1], ['y1, 'y2], [ last(yL1), last(ypL1)], ['x, 0, 1, h] ),
	 xL2 : take(outL,1),
	 yL2 : take(outL,2),
	 ypL2 : take(outL,3),	 
	 
	 
	 /*  construct xR, yR, and ypR for 1 < x < 1 + xdecay   */

	 outL : rk4(['y2, - glr* 'y1], ['y1, 'y2], [ 0, -ypright], ['x, 1 + xdecay, 1, -h] ),
	 xR : take(outL,1),
	 yR : take(outL,2),
	 ypR : take(outL,3),	 
	 
	 xR : reverse(xR),
	 yR : reverse(yR),
	 ypR : reverse(ypR),
	 if wfdebug then print(" yR(1) = ", first(yR)),
	 
	 fac : first(yR)/last(yL2),
	 if wfdebug then print(" fac = ",fac),
	 
	 yL1 : fac*yL1,
	 ypL1 : fac*ypL1,
	 yL2 : fac*yL2,
	 ypL2 : fac*ypL2,
	 
	 done)$
	 
/*
(%i44) load(FW);
 gam2 =  2500 
 h =  0.01 ,   xdecay =  2 ,  ypleft =  1.0E-28   ypright =  1.0E-28 
(%o44) "c:/k3/FW.mac"
(%i45) wf(0.5);
(%o45) done

(%i47) fll(xL1);
(%o47) [-2.0,0.0,201]
(%i48) fll(yL1);
(%o48) [0.0,-5.0917742,201]
(%i49) fll(ypL1);
(%o49) [-7.08074296E-29,-180.0214,201]
(%i50) fll(xL2);
(%o50) [0.0,1.0,101]
(%i51) fll(yL2);
(%o51) [-5.0917742,7.1910168,101]
(%i52) fll(ypL2);
(%o52) [-180.0214,-2.0395497,101]
(%i53) fll(xR);
(%o53) [1.0,3.0,201]
(%i54) fll(yR);
(%o54) [7.1910168,0.0,201]
(%i55) fll(ypR);
(%o55) [-254.24084,-1.0E-28,201]
*/


/*  count the number of nodes in yL2  */

num_nodes() :=
block([ n, numer], numer:true,   		
	n : 0,
	for j thru (length(yL2) - 1) do 
	    if yL2[ j ] * yL2[ j + 1 ] < 0 then n : n + 1,
	n)$
	
	

 
/*  dy_diff() uses global ypL2, ypR, and yR, 
     returns a normalized difference of derivatives
	 (yL'(1) - yR'(1)/ yR(1)
*/

dy_diff() :=
block([dy_left, dy_right,  numer],numer:true,
     dy_left : last(ypL2),
	 dy_right : first(ypR),
	 (dy_left - dy_right)/abs(first(yR)) )$


/* dy_diff2() prints out values of dy_left and
      dy_right, returns same as dy_diff()  
*/	  

dy_diff2() :=
block([dy_left, dy_right,  numer],numer:true,
     dy_left : last(ypL2),
	 dy_right : first(ypR),
	 print(" dy_left = ", dy_left, ",  dy_right = ", dy_right),
	 (dy_left - dy_right)/abs(first(yR)) )$



/* energy eigenvalue if global function F(E) = 0 .
  F(E) calls wf(E) then returns dy_diff(), but
    returns false if E < 0 or > 1 .  
 */	
	
	
F(E) :=
block( [ numer],numer:true,
    if E < 0 or E > 1 then (
         print(" in F(E), E = ",E,"  should be between 0 and 1 "),
  	    return(false)),
    wf(E),	
	dy_diff())$
    
	
/*
(%i65) F(0.5);
(%o65) 35.071714
(%i68) EL : makelist(e,e,0.1,0.9,0.01)$
(%i69) FL : map(F, EL)$	
(%i71) plot2d([discrete,EL,FL],[xlabel,"E"],[ylabel,"F(E)"],
          [gnuplot_preamble,"set grid"])$
		  
*/

/* 
     F1(k)
    energy eigenvalue if global function F1(k) = 0 .
    F1(k) calls wf(E) then returns dy_diff(), but
         returns false if k <= 0 or k >=  gam .  
 */	
	
	
F1(k) :=
block( [ numer],numer:true,
    if k <= 0 or k >= gam then (
	     print(" in F1(k), k = ",k,"  k should be greater than 0 and less than ",gam),
		 return(false) ),    
    wf(k^2/gam2),	
	dy_diff())$
	
/*
(%i10) F1(2.9);
(%o10) 34.29503
(%i12) F1(3.05);
(%o12) -49.882755
(%i13) [ka,kb] : bracket(F1,2.9,0.05,0.02);
(%o13) [3.0125,3.025]
(%i15) k0 : find_root(F1,ka,kb);
(%o15) 3.0206914
(%i16) E0 : ktoE(k0);
(%o16) 0.00364983


(%i23) kL : makelist(k,k,1,49,0.5)$
(%i24) F1L : map(F1, kL)$
(%i25) time(%);
(%o25) [6.14]
(%i26) plot2d([discrete,kL,F1L],[xlabel,"k"],[ylabel,"F1(k)"])$
*/

   


/*  bracket is a modified version of bracket_basic, designed to work with
    the functions F(E) or F1(k) 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(func,xx,dxx,xacc) :=
block([f1,f2, x:xx, dx:dxx,xx1,xx2,it:0,itmax:1000],        	
    
    do (
	  it : it + 1,	  
	  if debug then print(it),
	  if it > itmax then (
          print(" can't find change in sign "),    
	      return([0, 0 ])),
	  xx1 : x,
	  xx2 : x + dx,
	  
	  if debug then print(" it = ",it," xx1 = ",xx1," xx2 = ",xx2," dx = ", dx),
	  f1 : func(xx1),
	  if not f1 then (
	      print(" in bracket, f1 = false , xx1 = ",xx1, "  dx = ", dx),
		  return(f1)),
		  
	  f2 : func(xx2),
	  if not f2 then (
	      print(" in bracket, f2 = false ,  xx2 = ",xx2, "  dx = ", dx),
		  return(f2)),
	  if debug then print (" f1 = ",f1,"  f2 = ",f2),  
	  if f1*f2 < 0 then (
	       if abs(dx) < xacc then return([xx1,xx2]),
		   x : x - dx,
		   dx : dx/2)
	  else x : xx2))$
	  
/*
zero node solution

(%i36) [ea,eb] : bracket(F,0.0005,0.0001,0.00005);
(%o36) [0.003625,0.00365]
(%i37) e : find_root(F,ea,eb);
(%o37) 0.00364983
(%i38) wf(e);
(%o38) done
(%i39) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$



-----------------------------------------
1 node solution

(%i26) load(FW);
 gam2 =  2500 
 h =  0.01 ,   xdecay =  0.5 ,  ypleft =  1.0E-8   ypright =  1.0E-8 
(%o26) "c:/k3/FW.mac"
(%i27) [ea,eb] : bracket(F,0.01,0.005,0.001);
(%o27) [0.014375,0.015]
(%i28) e : find_root(F,ea,eb);
(%o28) 0.0145973
(%i29) wf(e);
(%o29) done
(%i30) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
-------------------------------------------------------


(%i26) load(FW);
 gam2 =  2500 
 h =  0.01 ,   xdecay =  0.5 ,  ypleft =  1.0E-8   ypright =  1.0E-8 
(%o26) "c:/k3/FW.mac"

3 node solution

(%i4) [ea,eb] : bracket(F,0.01,0.01,0.005);
(%o4) [0.0575,0.06]
(%i5) e : find_root(F,ea,eb);
(%o5) 0.0583554
(%i6) wf(e);
(%o6) done
(%i7) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
			 
			 
*/

 /* normalize() uses the current global xL1,yL1, xL2,yL2, xR, yR and
  the  utility functions  merge and simp  to define global
  xn and yn, with the latter being normalized.
  */
  
  

normalize() :=
block ( [AA,x_mean,x2_mean,delx,delx2, numer ], numer:true,     
     
	 xn :  merge( xL1, merge( rest(xL2), rest(xR))),
	 yn :  merge( yL1, merge( rest(yL2), rest(yR))),
	 
	 /* we need xn to have odd # of elements to use simp */
	 if evenp ( length (xn) ) then (
	     xn : rest (xn),
		 yn : rest (yn)),	 
	 if debug then print ( " fll(xn) = ", fll(xn) ),
	 if debug then print( " fll(yn) = ", fll(yn) ),
	 AA : simp(xn,yn^2),			 
	 print( " AA = ",AA),	
	 yn : yn/sqrt(AA),
	 if debug then print( " fll(yn) = ", fll(yn) ),
	 
	 x_mean : simp(xn, xn * yn^2),
     print(" x_mean = ", x_mean),							 

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

     delx2 : x2_mean - x_mean^2,   /* this should be positive!  */

     delx : sqrt(delx2),
	 
	 print(" delx = ", delx),
	 done)$
	 
/*

(%i52) [ka,kb] : bracket(F1,3.03,0.01,0.005);
(%o52) [3.0775,3.08]
(%i53) kv : find_root(F1,ka,kb);
(%o53) 3.0799546
(%i54) Ev : ktoE(kv);
(%o54) 0.00379445
(%i55) wf(Ev);
(%o55) done
(%i56) num_nodes();
(%o56) 0
(%i57) dy_diff();
(%o57) -1.60000886E+16
(%i58) normalize();
 AA =  6.40428117E+32 
 x_mean =  0.489981 
 delx =  0.184389 
(%o58) done
(%i59) fll(xn);
(%o59) [-0.5,1.5,201]
(%i60) fll(yn);
(%o60) [0.0,0.0,201]
(%i61) head(xn);
(%o61) [-0.5,-0.49,-0.48,-0.47,-0.46,-0.45]
(%i62) simp(xn,yn^2);
(%o62) 1.0

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

*/

	 
yn_plot_current() :=
block([ymn, ymx, numer],numer:true,
	 ymn : float(floor( lmin(yn) )),
	 ymx : float( 1 + floor( lmax (yn) ) ),
	 print(" ymax = ", lmax(yn) ),
	 plot2d( [discrete, xn, yn], ['y,ymn,ymx],
        	 [ylabel,"y"],  [xlabel,"x"],
		 [style,[lines,3]], [legend, false],  [gnuplot_preamble,"set grid"]))$
		 
/* numerical ps help:
     yn_plot_ps() creates global ymn, ymx assoc
          with current **numerical** yn  
*/	 
		   
yn_plot_ps() :=
block([numer],numer:true,
     ymn : float(floor( lmin(yn) )),
	 ymx : float( 1 + floor( lmax (yn) ) ))$
	 
     
	 
/*
(%i41) yn_plot_ps()$
(%i42) ymn;
(%o42) 0.0
(%i43) ymx;
(%o43) 2.0
(%i45) plot2d( [discrete, xn, yn], ['y,ymn,ymx],
        	 [ylabel,"y"],  [xlabel,"x"],
		 [style,[lines,3]], [legend, false],  [gnuplot_preamble,"set grid"],
             [psfile,"c:/k3/fw7.eps"]);
(%o45) "c:/k3/fw7.eps"
*/	 		 
		 

/*   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 lists xn and yn. Finally makes a plot
	 of yn over only the region (xmn, xmx)   */
					  
yn_plot(E,xmn,xmx) :=
block([ymn, ymx, numer],numer:true,
     wf(E),
	 print(" E = ",E ),
	 print(" number of nodes = ",num_nodes(),",       dy_diff = ",dy_diff() ),
	 normalize(),
	 print(" normalized ymax = ", lmax(yn) ),
	 ymn : floor( lmin(yn) ),
	 ymx : 1 + floor( lmax (yn) ),
	 plot2d( [discrete, xn, yn], ['x,xmn, xmx],
				 ['y,ymn,ymx],  [ylabel,"y"],  [xlabel,"x"],
				  [style, [lines, 3] ], [legend, false],  [gnuplot_preamble,"set grid"]))$
				  
/*

(%i1) load(cp3);
(%o1) "c:/k3/cp3.mac"
(%i2) load(FW);
  gam =  50     gam2 =  2500 
 h =  0.01 ,   xdecay =  0.5 ,  ypleft =  1.0E-8   ypright =  1.0E-8 
(%o2) "c:/k3/FW.mac"
(%i3) [ea,eb] : bracket(F,0.0005,0.0001,0.00005);
(%o3) [0.003625,0.00365]
(%i4) E0 : find_root(F,ea,eb);
(%o4) 0.00364983

(%i15) yn_plot(E0,-0.5,1.5)$
 E =  0.00364983 
 number of nodes =  0 ,       dy_diff =  -2.0630713E-12 
 AA =  6653.4608 
 x_mean =  0.5 
 delx =  0.188088 
 normalized ymax =  1.3866201
 
 help to make eps plot
 
 (%i14) ymn : floor( lmin(yn) );
(%o14) 0
(%i15) ymx : 1 + floor( lmax (yn) );
(%o15) 2
(%i16) plot2d( [discrete, xn, yn], ['x,-0.5, 1.5],
	['y,ymn,ymx],  [ylabel,"y"],  [xlabel,"x"],
	[style, [lines, 5] ], [legend, false],  [gnuplot_preamble,"set grid"],
        [psfile,"c:/k3/fw3.eps"]);
(%o16) "c:/k3/fw3.eps"

*/

/* levels(kmin,kmax,dk, kacc ) returns a list  [Ea, Eb,...] of energy levels with
     increasingly larger number of nodes in energy range (Emin, Emax)
	 according to Emin = kmin^2/gam^2, and Emax = kmax^2/gam^2.
	 uses F1(k) (inside bracket) to find  roots, and calls wf(E), num_nodes() and dy_diff() for each root found,
    Uses bracket and find_root.
	The arguments (dk, kacc) 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.
	Searching for the k eigenvalues via F1(k) is easier than searching 
	  directly for the E eigenvalues via F(E) for the case gam = 50 we
	  consider in this example.
*/
	

levels(kmin,kmax,dk, kacc ) := 
block([ k,knext, kroot,eroot, eL, ka, kb, nn, nlast : -1, r, numer], numer:true,
     k : kmin,	      
	 eL : [ ],     /*    list eL will hold energy eigenvalues found  */
	 
	 do (   	 if k > kmax then return(),    /* exit do loop  */
	         print("---------------- levels -------------------"),
	        print(" nlast = ", nlast),
	        print(" kstart = ", k," dk = ", dk ),			
            [ka, kb] : bracket(F1,k, dk, kacc),	 
			print(" ka = ",ka," kb = ",kb),
			if float(ka) = 0.0 then (
			    print(" can't find bracket interval "),
				print(" k = ",k),
				return() ),
			kroot : find_root(F1, ka, kb),
            print(" kroot = ", kroot),
			eroot : kroot^2/gam2,
			print(" eroot = ", eroot),
			wf(eroot),
			nn : num_nodes(),
			print(" number of nodes = ", nn),
			print(" dy_diff at x = 1 is    ", dy_diff() ),
			eL : cons(eroot, eL),
            nlast : nn,				
			r : read (" input c; or s; "),
			if string(r) = "s" then return(),	/* exit do loop  */		

			/* search for a k value greater than kb which produces
			      a wave function with nn + 1 nodes */
			knext : kb + dk,
			do (
			        wf(knext^2/gam2),
					if num_nodes() > nlast then (
					     k : knext,
						 return() )
					else knext : knext + dk)),
			
	 reverse(eL) )$
	 
/*
(%i5) levels(1,8,0.05,0.02);
---------------- levels ------------------- 
 nlast =  -1 
 kstart =  1  dk =  0.05 
 ka =  3.0125  kb =  3.025 
 kroot =  3.0206914 
 eroot =  0.00364983 
 number of nodes =  0 
 dy_diff at x = 1 is     -2.49565076E-14 
 input c; or s;  
c;
---------------- levels ------------------- 
 nlast =  0 
 kstart =  3.125  dk =  0.05 
 ka =  6.0375  kb =  6.05 
 kroot =  6.040956 
 eroot =  0.0145973 
 number of nodes =  1 
 dy_diff at x = 1 is     2.18273877E-13 
 input c; or s;  
c;
---------------- levels ------------------- 
 nlast =  1 
 kstart =  6.2  dk =  0.05 
 ka =  9.05  kb =  9.0625 
 kroot =  9.0603555 
 eroot =  0.032836 
 number of nodes =  2 
 dy_diff at x = 1 is     7.10273285E-14 
 input c; or s;  
c;
(%o5) [0.00364983,0.0145973,0.032836]

*/


/*
Examples which use F(E) instead of F1(k):

   /* zero node solution */

(%i36) [ea,eb] : bracket(F,0.0005,0.0001,0.00005);
(%o36) [0.003625,0.00365]
(%i37) e : find_root(F,ea,eb);
(%o37) 0.00364983
(%i38) wf(e);
(%o38) done
(%i39) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$



-----------------------------------------
    /* 1 node solution */

(%i26) load(FW);
 gam2 =  2500 
 h =  0.01 ,   xdecay =  0.5 ,  ypleft =  1.0E-8   ypright =  1.0E-8 
(%o26) "c:/k3/FW.mac"
(%i27) [ea,eb] : bracket(F,0.01,0.005,0.001);
(%o27) [0.014375,0.015]
(%i28) e : find_root(F,ea,eb);
(%o28) 0.0145973
(%i29) wf(e);
(%o29) done
(%i30) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
-------------------------------------------------------

spurious 1 node soln

(%i40) [ea,eb] : bracket(F,0.015,0.0001,0.00005);
(%o40) [0.015175,0.0152]
(%i41) e : find_root(F,ea,eb);
(%o41) 0.0151767
(%i42) wf(e);
(%o42) done
(%i43) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
----------------------------------------------------
/*  2 node solution  */

(%i45) [ea,eb] : bracket(F,0.017,0.0001,0.00005);
(%o45) [0.032825,0.03285]
(%i46) e : find_root(F,ea,eb);
(%o46) 0.032836
(%i47) wf(e);
(%o47) done
(%i48) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
---------------------------------------------------
/*  3 node solution */

(%i4) [ea,eb] : bracket(F,0.01,0.01,0.005);
(%o4) [0.0575,0.06]
(%i5) e : find_root(F,ea,eb);
(%o5) 0.0583554
(%i6) wf(e);
(%o6) done
(%i7) plot2d([[discrete,xL1,yL1],[discrete,xL2,yL2],[discrete,xR,yR]],
             [legend,false])$
			 
*/