;+
; NAME:
;       SEUTIL
;
; PURPOSE:
;
;       This simple widget program solves the one-dimensional Schrodinger
;	equation using the discrete variable approximation.  User can type in
;	a function of x that represents the potential, in IDL syntax, and the
;	eigenvalues and probability density are calculated for that potential.
;	The accuracy increases as the size of the Hamiltonian increases but at
;	the cost of processing speed.
;
; AUTHOR:
;
;       Robert M. Dimeo, Ph.D.
;	NIST Center for Neutron Research
;       100 Bureau Drive
;	Gaithersburg, MD 20899
;       Phone: (301) 975-8135
;       E-mail: robert.dimeo@nist.gov
;       http://www.ncnr.nist.gov/staff/dimeo
;
; CATEGORY:
;
;       Widgets, mathematics
;
; CALLING SEQUENCE:
;
;       SEUTIL
;
;
; INPUT FIELDS:
;
;	Note that a carriage return in some of the text fields below
;	begins the calculation
;
;	HAMILTONIAN SIZE:	size of the matrix to be diagonalized (default: 200)
;				Increase this to improve accuracy.
;	# EIGENFUNCTIONS:	Number of eigenfunctions to plot to the screen.
;	XLO:	lower bound for evaluation of the wavefunctions and potential
;	XHI:	upper bound for evaluation of the wavefunctions and potential
;	YLO:	lower bound for vertical display
;	YHI:	upper bound for vertical display
;	MASS:	mass of particle in atomic mass units (AMU)
;	POTENTIAL:	function describing the potential (meV)
;
; EXAMPLES FOR POTENTIALS (must use IDL syntax as shown below):
;
;	0.5*x^2				; simple harmonic oscillator
;	(abs(2.0*x))^15			; approximation to the infinite square well
;	abs(x)				; linear potential
;	20.0*((x/2.0)^2-1.0)^2		; quartic potential (exhibits tunneling)
;
;	I have also included a step function in this program so that you
;	can see the effects of confinement in a square well.  An example
;	is listed below and is the default when the program is launched.
;
;	10.0+10.0*(step(x-1.0)-step(x+1.0))	; finite square well
;
; DISCLAIMER
;
;	This software is provided as is without any warranty whatsoever.
;	Permission to use, copy, modify, and distribute modified or
;	unmodified copies is granted, provided this disclaimer
;	is included unchanged.
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, September 30, 2002.
;
;	-Added built-in potential functions as menu items under 
;	 POTENTIALS (10/01/02)
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function step,x
a = 1D & b = 0D & c = 0.5D
y = (x gt 0)*a + ((x lt 0))*b + c*(x eq 0)
return,y
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function kronMat,n,i,j
i = fix(i) & j = fix(j)
k = 0*indgen(n,n)
ni = n_elements(i)
for ii = 0,ni-1 do begin
    if (i[ii] lt n) and (i[ii] ge 0) and $
       (j[ii] lt n) and (j[ii] ge 0) then begin
       k[i[ii],j[ii]] = 1
    endif
endfor
return,k
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro seCleanup,tlb
widget_control,tlb,get_uvalue = pState
tvlct,*(*pState).r,*(*pState).g,*(*pState).b
ptr_free,(*pState).r,(*pState).g,(*pState).b
ptr_free,(*pState).evalsPtr,(*pState).probPtr,(*pState).xPtr,(*pState).vPtr
wdelete,(*pState).winPix
ptr_free,pState
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro seutilQuit,event
widget_control,event.top,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro sePlot,event
widget_control,event.top,get_uvalue = pState
if (n_elements(*(*pState).vPtr) eq 0) or $
   (n_elements(*(*pState).xPtr) eq 0) then return
evals = *(*pState).evalsPtr
evecs = *(*pState).probPtr
x = *(*pState).xPtr
v = *(*pState).vPtr
if (*pState).autoscale eq 1 then begin
  widget_control,(*pState).yhi,get_value = yhi
  widget_control,(*pState).ylo,get_value = ylo
  xlo = min(x) & xhi = max(x)
  (*pState).xrange = [xlo,xhi]
  (*pState).yrange = [float(ylo[0]),float(yhi[0])]
endif

plot,x,v,xrange = (*pState).xrange,yrange = (*pState).yrange, $
     xstyle = 1,ystyle = 1,xtitle = '!6x ( !6!sA!r!u!9 %!6!n )', $
     ytitle = '!6P(x)', title = '!6Probability Density',linestyle = 0, $
     thick = 3.0

widget_control,(*pState).nvals,get_value = nvals
nvals = fix(nvals[0]) < n_elements(v)
;nvals = (*pState).numVals2Plot
for i = 0,nvals-1 do begin
  oplot,x,evals[i]+evecs[*,i]
endfor
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro seSolve,event
!except = 0
widget_control,event.top,get_uvalue = pState
strout = ['Solving the Schrodinger equation','Please wait...']
widget_control,(*pState).info,set_value = strout

widget_control,(*pState).hsize,get_value = nx & nx = fix(nx[0])
widget_control,(*pState).xlo,get_value = xlo & xlo = double(xlo[0])
widget_control,(*pState).xhi,get_value = xhi & xhi = double(xhi[0])
widget_control,(*pState).pot,get_value = pot & pot = string(pot[0])
widget_control,(*pState).mass,get_value = m & m = double(m[0])
hbarc = 1973d3 & mc2 = m*931.5d9
beta = -(hbarc^2)/(2.0*mc2)

dx = (xhi-xlo)/(nx-1.0)
x = xlo+dx*dindgen(nx)
*(*pState).xPtr = x
vstring = 'v='+strtrim(pot,2)
r = execute(vstring,1)
if r ne 1 then begin
  strout = 'Input error'
  void = dialog_message(dialog_parent = event.top,strout)
  return
endif

*(*pState).vPtr = v
; Build the hamiltonian
h = dblarr(nx,nx)
ux = 1+bytarr(nx)
ix = indgen(nx)

h = (ux#(v-2.0*beta/dx^2))*kronMat(nx,ix,ix) + $
       (ux#(beta/dx^2)#ux)*kronMat(nx,ix,ix+1) + $
       (ux#(beta/dx^2)#ux)*kronMat(nx,ix,ix-1)


trired,h,evals,e,/double
triql,evals,e,h,/double

esort = sort(evals)
evals = evals[esort]
evecs = h[*,esort]
*(*pState).evalsPtr = evals
prob = evecs*evecs

widget_control,(*pState).ylo,get_value = ylo
widget_control,(*pState).yhi,get_value = yhi
ylo = float(ylo[0]) & yhi = float(yhi[0])

yfactor = yhi-ylo
nvals = (*pState).numVals2Plot
for i = 0,nvals-1 do begin
  sf = 0.05*yfactor/(max(prob[*,i])-min(prob[*,i]))
  prob[*,i] = sf*prob[*,i]
endfor

*(*pState).probPtr = prob
trans = string(evals[1:nx-1]-evals[0:nx-1])
feig = '(f10.2)'
f = '(e10.3)'
trans = ['--',trans]
strout = strarr(nx+2)
strout[0] = 'Eigenvalue (meV), Transition (meV)'
strout[1] = '----------------------------------'
strout[2] = strtrim(string(evals[0],format = feig),2)+ $
              ',  ------'
for i = 1,nx-1 do begin
  strout[2+i] = strtrim(string(evals[i],format = feig),2)+ $
              ',  '+strtrim(string(trans[i],format = f),2)
endfor
widget_control,(*pState).info,set_value = strout

wset,(*pState).winPix
sePlot,event
wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro sedraw,event
widget_control,event.top,get_uvalue = pState

case event.type of
0:	begin		; button press
	  (*pState).mouse = event.press
	  if (*pState).mouse eq 4 then begin
	    (*pState).autoscale = 1

	    wset,(*pState).winPix
	    sePlot,event
	    wset,(*pState).winVis
	    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
	  endif
	  if (*pState).mouse eq 1 then begin
	    (*pState).xbox[0] = event.x
	    (*pState).ybox[0] = event.y
	    wset,(*pState).winVis
	    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
	    empty
	    (*pState).autoscale = 0
	    widget_control,(*pState).win,/draw_motion_events
	  endif
	end
1:	begin	; button release
	 if (*pState).mouse eq 1 then begin
	  xll = (*pState).xbox[0] < (*pState).xbox[1]
	  yll = (*pState).ybox[0] < (*pState).ybox[1]
	  w = abs((*pState).xbox[1] - (*pState).xbox[0])
	  h = abs((*pState).ybox[1] - (*pState).ybox[0])
	  xur = xll + w
	  yur = yll + h
	  ll = convert_coord(xll,yll,/device,/to_data)
	  ur = convert_coord(xur,yur,/device,/to_data)
	  (*pState).xrange = [ll[0],ur[0]]
	  (*pState).yrange = [ll[1],ur[1]]
	  wset,(*pState).winPix
	  sePlot,event
	  wset,(*pState).winVis
	  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
	  (*pState).mouse = 0B
	  widget_control,(*pState).win,draw_motion_events = 0
	 endif
	 if (*pState).mouse eq 4 then begin
	  wset,(*pState).winPix
	  sePlot,event
	  wset,(*pState).winVis
	  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
	  (*pState).mouse = 0B
	  widget_control,(*pState).win,draw_motion_events = 0
	 endif
	end
2:	begin	; mouse motion
	  if (*pState).mouse eq 1 then begin
	  	(*pState).xbox[1] = event.x
	  	(*pState).ybox[1] = event.y
	  	xc = [(*pState).xbox[0],event.x,event.x,$
	  	      (*pState).xbox[0],$
	  	      (*pState).xbox[0]]
	  	yc = [(*pState).ybox[0],(*pState).ybox[0],$
	  	      event.y,event.y,$
	  	      (*pState).ybox[0]]
	  	wset,(*pState).winVis
	  	device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
	  	plots,xc,yc,/device
	  	empty
	  endif
	end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro se_sizeEvents,event
widget_control,event.top,get_uvalue = pState
if (event.id eq (*pState).hsize) or $
   (event.id eq (*pState).xlo) or $
   (event.id eq (*pState).xhi) or $
   (event.id eq (*pState).mass) or $
   (event.id eq (*pState).pot) then begin
   seSolve,event
   return
endif
if (event.id eq (*pState).ylo) or $
   (event.id eq (*pState).nvals) or $
   (event.id eq (*pState).yhi) then begin
   wset,(*pState).winPix
   sePlot,event
   wset,(*pState).winVis
   device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
   return
endif
ctrlgeom = widget_info((*pState).base,/geometry)
tlbgeom = widget_info(event.top,/geometry)
xsize = event.x
ysize = event.y

newxsize = xsize-ctrlgeom.xsize
newysize = ysize

widget_control,(*pState).win,draw_xsize = newxsize, $
               draw_ysize = newysize
wdelete,(*pState).winPix
window,/free,/pixmap,xsize = newxsize,ysize = newysize
(*pState).winPix = !d.window

wset,(*pState).winPix
sePlot,event
wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEDefaults,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).xlo,set_value = '-5.0'
widget_control,(*pState).xhi,set_value = '5.0'
widget_control,(*pState).hsize,set_value = '200'
widget_control,(*pState).nvals,set_value = '10'
widget_control,(*pState).ylo,set_value = '0.0'
widget_control,(*pState).yhi,set_value = '20.0'
widget_control,(*pState).mass,set_value = '1.0'
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEloadSHO,event
widget_control,event.top,get_uvalue = pState
seDefaults,event
strout = '0.5*x^2'
widget_control,(*pState).pot,set_value = strout
seSolve,event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEloadSquareWell,event
widget_control,event.top,get_uvalue = pState
seDefaults,event
strout = '10.0+10.0*(step(x-1.0)-step(x+1.0))'
widget_control,(*pState).pot,set_value = strout
seSolve,event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEloadLinearWell,event
widget_control,event.top,get_uvalue = pState
seDefaults,event
strout = 'abs(4.0*x)'
widget_control,(*pState).pot,set_value = strout
seSolve,event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEloadQuarticWell,event
widget_control,event.top,get_uvalue = pState
seDefaults,event
strout = '10.0*((x/2.0)^2-1.0)^2'
widget_control,(*pState).pot,set_value = strout
seSolve,event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro SEloadDoubleWell,event
widget_control,event.top,get_uvalue = pState
seDefaults,event
widget_control,(*pState).mass,set_value = '5.0'
strout = '10.0+10.0*(step(x-2.0)-step(x+2.0)) + 5.0*(step(x+0.5)-step(x-0.5))'
widget_control,(*pState).pot,set_value = strout
seSolve,event
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro seutil,group_leader = group_leader
tvlct,r,g,b,/get
loadct,0,/silent
; Widget definition module
if n_elements(group_leader) eq 0 then begin
  tlb = widget_base(/row,title = 'One Dimensional Schrodinger Equation Solver', $
        /tlb_size_events,mbar = bar)
endif else begin
  tlb = widget_base(/row,title = 'One Dimensional Schrodinger Equation Solver', $
        group_leader = group_leader,/tlb_size_events,mbar = bar)
endelse

filemenu = widget_button(bar,value = 'FILE',/menu)
potMenu = widget_button(bar,value = 'POTENTIALS',/menu)

void = widget_button(potMenu,value = 'Simple harmonic oscillator', $
       event_pro = 'SEloadSHO')
void = widget_button(potMenu,value = 'Finite square well', $
       event_pro = 'SEloadSquareWell')
void = widget_button(potMenu,value = 'Linear well', $
       event_pro = 'SEloadLinearWell')
void = widget_button(potMenu,value = 'Quartic well', $
       event_pro = 'SEloadQuarticWell')
void = widget_button(potMenu,value = 'Double square well', $
       event_pro = 'SEloadDoubleWell')

if !d.name eq 'WIN' then begin
  thisFont = "Comic Sans MS*16"
  bigFont = "Comic Sans MS*20*Bold"
endif else begin
  thisFont = -1
  bigFont = -1
endelse

base = widget_base(tlb,/col)
newBase = widget_base(base,/row)
thisbase = widget_base(newbase,/col,/grid_layout)
hsize = cw_field(thisbase,/col,title = 'Hamiltonian size',value = 200,$
        /string,font = thisFont,/return_events)
nvals = cw_field(thisbase,/col,title = '# of eigenfunctions',value = 10,$
        /string,font = thisFont,/return_events)
xlo = cw_field(thisbase,/col,title = 'XLO',value = -5.0,/string,$
      font = thisFont,/return_events)
xhi = cw_field(thisbase,/col,title = 'XHI',value = 5.0,/string,font = thisFont, $
      /return_events)
ylo = cw_field(thisbase,/col,title = 'YLO',value = 0.0,/string,$
      font = thisFont,/return_events)
yhi = cw_field(thisbase,/col,title = 'YHI',value = 20.0,/string,font = thisFont, $
      /return_events)
mass = cw_field(thisbase,/col,title = 'MASS (AMU)',value = '1.0',/string,$
       font = thisFont,/return_events)

void = widget_button(filemenu,value = 'QUIT',event_pro = 'seutilQuit', $
       font = thisFont)
info = widget_text(newbase,value = '',/editable,/scroll,xsize = 30,$
       ysize = 30,font = thisFont)
pot = cw_field(base,/col,title = 'Potential: V(x)', $
      value = '10.0+10.0*(step(x-1.0)-step(x+1.0))',/string,$
      fieldfont = bigFont,xsize = 45,/return_events, $
      font = thisFont)
xsize = 500 & ysize = 400
win = widget_draw(tlb,xsize = xsize,ysize = ysize,/button_events, $
      event_pro = 'sedraw')

widget_control,tlb,/realize
geom = widget_info(base,/geometry)
newysize = geom.ysize
widget_control,win,draw_ysize = newysize

widget_control,win,get_value = winVis
window,/free,/pixmap,xsize = xsize,ysize = newysize
winPix = !d.window

state = {winVis:winVis, $
         winPix:winPix, $
         autoscale:1, $
         xrange:[0.0,1.0], $
         yrange:[0.0,1.0], $
         xbox:[0.0,1.0], $
         ybox:[0.0,1.0], $
         mouse:0B, $
         numVals2Plot:100, $
         win:win, $
         xlo:xlo, $
         xhi:xhi, $
         ylo:ylo, $
         yhi:yhi, $
         pot:pot, $
         nvals:nvals, $
         mass:mass, $
         base:base, $
         hsize:hsize, $
         info:info, $
         r:ptr_new(r), $
         g:ptr_new(g), $
         b:ptr_new(b), $
         xPtr:ptr_new(/allocate_heap), $
         vPtr:ptr_new(/allocate_heap), $
         evalsPtr:ptr_new(/allocate_heap), $
         probPtr:ptr_new(/allocate_heap)}


pState = ptr_new(state,/no_copy)
widget_control,tlb,set_uvalue = pState

seSolve,{event,id:void,top:tlb,handler:0l}

xmanager,'seutil',tlb,event_handler = 'se_sizeEvents',$
         cleanup = 'seCleanup', /no_block

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;