; $Id: methyl_dyn.pro,v 1.25 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  METHYL_DYN
;
; PURPOSE:
;  See description below.
;
; CATEGORY:
;  DAVE, HFBS
;
; 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
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;
; METHYL_DYN.PRO
;
; Calculates the rotational dynamics of methyl
; group under the influence of a three-fold symmetric
; potential barrier.  This program uses CW_TABBASE.PRO.
;
; Modification history:
; --------------------
; - Written by R.M.Dimeo (02/26/02)
; - Added numerous enhancements in functionality such as
;   dynamic position sensitivity displaying the cursor
;   position and the associated eigenvalues and transitions.
;  RMD (02/27/02)
; - Incorporated the notifyIds construct so that this program
;   will work properly within the DAVE environment, RMD (02/28/02)
;
;
;;;;;;;;;;;;;;;;;;;
pro mdCleanup,tlb
widget_control,tlb,get_uvalue = pState
help,/heap,/brief
; Restore the colors
tvlct,(*pState).red,(*pState).green,(*pState).blue
; Restore the margins
!x.margin = [10.0,3.0]
!y.margin = [4.0,2.0]
; Delete all pixmaps
wdelete,(*pState).win1State.winPix
wdelete,(*pState).win2State.winPix
; Free up all of the pointers
ptr_free,(*pState).win1State.xPtr,(*pState).win1State.yPtr
ptr_free,(*pState).win2State.xPtr,(*pState).win2State.yPtr
ptr_free,(*pState).phiPtr,(*pState).probPtr
ptr_free,(*pState).potPtr,(*pState).v3Ptr
ptr_free,(*pState).eigPtr,(*pState).hPtr
ptr_free,(*pState).mPtr,(*pState).dynPtr
ptr_free,(*pState).evalsPtr,(*pState).eigFuncs
ptr_free,(*pState).psi_imPtr,(*pState).psi_rePtr
ptr_free,pState
help,/heap,/brief
;return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro mdEqsQuit,event
widget_control,event.top,get_uvalue = pState
tvlct,(*pState).ri,(*pState).gi,(*pState).bi
widget_control,event.top,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro mdEqs,group_leader = group_leader
tvlct,ri,gi,bi,/get
tlb = widget_base(group_leader = group_leader,floating = 1,$
      title = 'Description of Methyl Rotor Calculation',/col,map = 0,$
      tlb_frame_attr = 1)
;filename = file_which('mrotor.png',/include_current_dir)
filename = !DAVE_AUXILIARY_DIR+'mrotor.png'
image = read_png(filename,r,g,b)
s = size(image)
sx = s[1] & sy = s[2]

winxsize = sx & winysize = sy
plotWin = widget_draw(tlb,xsize = winxsize,ysize = winysize)
void = widget_button(tlb,value = 'DISMISS',event_pro = "mdEqsQuit")

centertlb,tlb
widget_control,tlb,/realize
widget_control,plotWin,get_value = eqVis
wset,eqVis

tvlct,r,g,b ; load the proper color tables as found in
; read_png above.
tvimage,image
widget_control,tlb,map = 1
state = {top:group_leader,ri:ri,gi:gi,bi:bi}
pState = ptr_new(state,/no_copy)
widget_control,tlb,set_uvalue = pState,/no_copy
xmanager,'mdEqs',tlb
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro mdHelp,event
mdEqs,group_leader = event.top
return
end
;;;;;;;;;;;;;;;;;;;
function resetColors,dummy
device,decomposed = 0
loadct,0,/silent
tvlct,255,255,255,1  ;white
tvlct,0,255,0,2     ;green
tvlct,255,0,0,3     ;red
tvlct,255,255,0,4 ;yellow
tvlct,122,203,255,5  ; blue
tvlct,255,152,77,6
gold = 6
red = 3
green = 2
yellow = 4
blue = 5
white = 1
colors = {red:red,green:green,white:white,yellow:yellow,blue:blue,gold:gold}
return,colors
end
;;;;;;;;;;;;;;;;;;;
function kronecker,i,j
return,double(i eq j)
end
;;;;;;;;;;;;;;;;;;;
pro mdQuit,event
widget_control,event.top,/destroy
return
end
;;;;;;;;;;;;;;;;;;;
pro mdPlotProb,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top,$
         'Perform calculation first!',/information)
  return
endif
widget_control,(*pState).nstate,get_value = nstate

widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])
nstate = fix(nstate[0])
eigvals = *(*pState).EigPtr
if nsym*nstate gt numpts then begin
  void = dialog_message('State out of range of calculation',/information)
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

Estate = eigvals[npts-1,nstate]
widget_control,(*pState).Estate,set_value = Estate

pi = !pi
nphi = (*pState).nphi
phi = *(*pState).PhiPtr

h = *(*pState).hPtr
m = *(*pState).mPtr
evals = *(*pState).evalsPtr

wfe = fltarr(nphi)
wfo = fltarr(nphi)
prob = fltarr(nphi)
for i=0,nphi-1 do begin
  wfe[i] = ((cos(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
  wfo[i] = ((sin(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
endfor

prob = (1.0/(2.0*pi))*(wfe^2+wfo^2)+0.0*evals[nstate]
phihi = max(phi) & phimin = min(phi)
thickness = 2.0

plot,phi,prob,xrange=[phimin,phihi],title='!6Eigenstate Probability Density',$
      color = (*pState).colors.white,$
      xstyle = 1,xtitle='!4 u !6 (rads)',yrange=[0.0,1.5],ystyle = 1,/nodata
oplot,phi,prob, color = (*pState).colors.green,thick = thickness
oplot,phi,*(*pState).PotPtr,linestyle = 0,color = (*pState).colors.yellow,$
      thick = thickness

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

widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top,$
                        'Perform calculation first!',/information)
  return
endif
widget_control,(*pState).nstate,get_value = nstate
widget_control,(*pState).nsym,get_value = nsym
nsym = nsym[0]

nstate = fix(nstate[0])
eigvals = *(*pState).EigPtr
if nsym*nstate gt numpts then begin
  void = dialog_message(dialog_parent = event.top,$
         'State out of range of calculation',/information)
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])
Estate = eigvals[npts-1,nstate]
widget_control,(*pState).Estate,set_value = Estate

pi = !pi
nphi = (*pState).nphi
phi = *(*pState).PhiPtr

h = *(*pState).hPtr
m = *(*pState).mPtr
evals = *(*pState).evalsPtr

wfe = fltarr(nphi)
wfo = fltarr(nphi)
prob = fltarr(nphi)

for i=0,nphi-1 do begin
  wfe[i] = 0.15*((cos(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
  wfo[i] = 0.15*((sin(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
endfor

minval = min([wfe[*],wfo[*],*(*pState).PotPtr])
maxval = 1.1*max([wfe[*],wfo[*],*(*pState).PotPtr])

phihi = max(phi) & phimin = min(phi)
thickness = 2.0
thickness1 = 3.0

plot,phi,wfe,xrange=[phimin,phihi],title='!6Eigenvectors',$
      color = (*pState).colors.white,$
      xstyle = 1,xtitle='!4 u !6 (rads)',yrange=[minval[0],maxval[0]],ystyle = 1,/nodata
oplot,phi,wfe, color = (*pState).colors.green,thick = thickness
oplot,phi,wfo, color = (*pState).colors.red,thick = thickness
oplot,phi,*(*pState).PotPtr,linestyle = 2,color = (*pState).colors.yellow,$
      thick = thickness1
plots,[0.6,0.7],[0.9,0.9],color = (*pState).colors.red,/normal,thick = thickness
plots,[0.6,0.7],[0.85,0.85],color = (*pState).colors.green,/normal,thick = thickness
plots,[0.6,0.7],[0.8,0.8],color = (*pState).colors.yellow,/normal,thick = thickness1
xyouts,0.75,0.9,'Imaginary',/normal
xyouts,0.75,0.85,'Real',/normal
xyouts,0.75,0.8,'Potential',/normal

return
end
;;;;;;;;;;;;;;;;;;;
pro mdPlotTrans,event

widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top,$
         'Solve Schrodinger equation first!',/information)
  return
endif

widget_control,(*pState).jinit,get_value = initState
initialState = fix(initState[0])

widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

V3 = *(*pState).V3Ptr
vmin = min(V3) & vmax = max(V3)
eigvals = *(*pState).EigPtr

widget_control,(*pState).nstate,get_value = nstate
nstate = fix(nstate[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

if nsym gt numpts then begin
  void = dialog_message(dialog_parent = event.top,$
         'State out of range of calculation',/information)
  return
endif
Estate = eigvals[npts-1,nstate]

widget_control,(*pState).Estate,set_value = Estate

if (*pState).win2State.autoscale eq 1 then begin
  (*pState).win2State.yrange = [0.0,max(eigvals[*,*])]
  (*pState).win2State.xrange = [vmin,vmax]
endif
thickness = 1.0
plot,V3,(eigvals[0:npts-1,initialState]),xtitle='!6V (meV)',ytitle='!6E (meV)',$
     yrange=(*pState).win2State.yrange,xrange=(*pState).win2State.xrange,$
     xstyle=1,ystyle=1,title='!6Rotational Transitions',/nodata,$
     color = (*pState).colors.white
for i=initialState+1,numpts-1 do begin
  oplot,V3,(eigvals[0:npts-1,i])-(eigvals[0:npts-1,initialState]),$
     color = (*pState).colors.white,$
     thick = thickness
endfor
*(*pState).win2State.xPtr = !x
*(*pState).win2State.yPtr = !y
return
end
;;;;;;;;;;;;;;;;;;;
pro mdPlotEig,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message('Solve Schrodinger equation first',/information)
  return
endif

widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

V3 = *(*pState).V3Ptr
vmin = min(V3) & vmax = max(V3)
eigvals = *(*pState).EigPtr

widget_control,(*pState).nstate,get_value = nstate
nstate = fix(nstate[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

if nsym gt numpts then begin
  void = dialog_message('State out of range of calculation',/information)
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])
Estate = eigvals[npts-1,nstate]

widget_control,(*pState).Estate,set_value = Estate

if (*pState).win1State.autoscale eq 1 then begin
  (*pState).win1State.yrange = [0.0,max(eigvals[*,*])]
  (*pState).win1State.xrange = [vmin,vmax]
endif

plot,V3,(eigvals[0:npts-1,0]),xtitle='!6V (meV)',ytitle='!6E (meV)',$
     yrange=(*pState).win1State.yrange,xrange=(*pState).win1State.xrange,$
     xstyle=1,ystyle=1,title='!6Energy Levels',/nodata;,color = (*pState).colors.white

*(*pState).win1State.xPtr = !x
*(*pState).win1State.yPtr = !y

number = 1.d10
roundedvalue = fltarr(numpts)

eigval_size = size(eigvals)
textIndex = 2
thickness = 1.25
for i = 0,eigval_size[2]-1 do begin
  oplot,V3,(eigvals[0:npts-1,i]),thick = thickness
endfor

xyouts,V3[textIndex],(eigvals[textIndex,0]),strtrim(string(0),2),charsize = 1.0,/data,$
       alignment = 0.0
for i = 1,numpts-1 do begin
  if (eigvals[textIndex,i] - eigvals[textIndex,i-1]) le 1.d-6 then begin
    xyouts,1.1*V3[textIndex],(eigvals[textIndex,i]),' ,'+strtrim(string(i),2),$
           charsize = 1.0,/data,$
           alignment = 0.0

  endif else begin
    xyouts,V3[textIndex],(eigvals[textIndex,i]),strtrim(string(i),2),$
           charsize = 1.0,/data,$
           alignment = 0.0
  endelse
endfor

return
end
;;;;;;;;;;;;;;;;;;;
pro mdCalc,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).label,set_value = 'Solving...'
widget_control,(*pState).B,get_value = B
B = double(B[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

widget_control,(*pState).v3hi,get_value = v3hi
v3hi = double(v3hi[0])


; additive barrier parameters
widget_control,(*pState).v6,get_value = v6
v6 = double(v6[0])
;widget_control,(*pState).nsym1,get_value = nsym1
;nsym1 = fix(nsym1[0])
;widget_control,(*pState).alpha,get_value = alpha
;alpha = !dtor*double(alpha[0])

widget_control,(*pState).nstate,set_value = 0
nstate = 0
pi = !pi
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

V3 = dblarr(npts)
energy1 = dblarr(npts)
energy2 = dblarr(npts)
energy3 = dblarr(npts)
energy4 = dblarr(npts)

;numpts = 31   ; dimension of hamiltonian matrix
energylevels = dblarr(npts,numpts)
cplxelevels = complexarr(npts,numpts)

n = lonarr(numpts)
m = lonarr(numpts)
mlo = -(numpts-1)/2
m = mlo+dindgen(numpts)
*(*pState).mPtr = m

n = m

V3lo = 0.0
dV = (V3hi-V3lo)/(npts-1)
V3 = V3lo+dV*dindgen(npts)

E3 = B   ; energy in meV for CH3 rotor

VDO = 0.0
nsym2 = 0.0

hre = dblarr(numpts,numpts)
him = dblarr(numpts,numpts)
h = complexarr(numpts,numpts)

for k = 0,npts-1 do begin
  strout = strtrim(string(k+1),2)+'/'+strtrim(string(npts),2)
  widget_control,(*pState).label,set_value = 'Progress...'
  widget_control,(*pState).solve,set_value = strout
  np = (1.0+bytarr(numpts))#n
  mp = m#(1.0+bytarr(numpts))
;;  hre = E3*(np^2)*kronecker(np,mp) + 0.25*V3[k]*(2.0*kronecker(np,mp)-$
;;               kronecker(np-mp+nsym,0)-kronecker(np-mp-nsym,0))+$
;;               0.25*v6*(2.0*kronecker(np,mp)-$
;;               kronecker(np-mp+6,0)-kronecker(np-mp-6,0))

  hre = E3*(np^2)*kronecker(np,mp) + 0.25*V3[k]*(2.0*kronecker(np,mp)-$
               kronecker(np-mp+nsym,0)-kronecker(np-mp-nsym,0))+$
               0.25*v6*(2.0*kronecker(np,mp)-$
               kronecker(np-mp+6,0)-kronecker(np-mp-6,0))

  h = complex(hre,him)
  trired,h,evals,e,/double
  triql,evals,e,h,/double
  esort = sort(evals)
  evals = evals[esort]

  cplxelevels[k,0:numpts-1] = evals
  energylevels[k,0:numpts-1] = evals[0:numpts-1]
endfor
widget_control,(*pState).label,set_value = 'Calculating energy levels'
;;;;;;;;;Calculate the positional probability density;;;;;;;;;;;;;;
nphi = (*pState).nphi
dphi = 2.0*pi/(nphi-1)
phi = fltarr(nphi)
phi = dphi*findgen(nphi)
wfo = fltarr(nphi)
wfe = fltarr(nphi)
prob = fltarr(nphi)

E0 = energylevels[npts-1,0*nsym]
E1 = energylevels[npts-1,1*nsym]
E2 = energylevels[npts-1,2*nsym]
Estate = energylevels[npts-1,nstate]
widget_control,(*pState).Estate,set_value = Estate

*(*pState).evalsPtr = evals
h = h[*,esort]
*(*pState).hPtr = h
*(*pState).ProbPtr = prob
*(*pState).PhiPtr = phi
*(*pState).EigPtr = energylevels[0:npts-1,0:numpts-1]
*(*pState).V3Ptr = V3
;if v6 gt 0.0 then begin
;  v3max = max(v3)
  ;ratio = v3max/v6
;  *(*pState).PotPtr = 0.3*(0.5*(1.0-cos(nsym*phi)));(0.0/ratio)*(1.0-cos(nsym1*phi+alpha)))
;endif else begin
  *(*pState).PotPtr = 0.3*0.5*(1.0-cos(nsym*phi))
;endelse
;(*pState).npts = npts
widget_control,(*pState).nstate,Set_Slider_Max = fix(numpts/nsym)
widget_control,(*pState).jinit,Set_Slider_Max = fix(numpts/nsym)
widget_control,(*pState).jfinal,Set_Slider_Max = fix(numpts/nsym)
; Calculate the wavefunctions
psi_im = dblarr(numpts,nphi)
psi_re = dblarr(numpts,nphi)

phivec = 1.0+bytarr(nphi)
numvec = 1.0+bytarr(numpts)
widget_control,(*pState).label,set_value = 'Calculating wavefunctions'
for i = 0,numpts-1 do begin
  psi_im[i,*] = ((sin((phivec#(fix(m)))*(phi#numvec))))#(((h[*,i])))
  psi_re[i,*] = ((cos((phivec#(fix(m)))*(phi#numvec))))#(((h[*,i])))
endfor

eigfuncs = dcomplexarr(numpts,nphi)
eigfuncs = dcomplex(psi_re[0:numpts-1,0:nphi-1],psi_im[0:numpts-1,0:nphi-1])
*(*pState).eigfuncs = eigfuncs
*(*pState).psi_ImPtr = psi_im
*(*pState).psi_RePtr = psi_re

wset,(*pState).win1State.winPix
mdPlotEig,event
wset,(*pState).win1State.winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
*(*pState).win1State.xPtr = !x
*(*pState).win1State.yPtr = !y
widget_control,(*pState).funcgroup,set_value = 2
empty
widget_control,(*pState).win1State.win,draw_motion_events = 1

wset,(*pState).win2State.winPix
mdPlotTrans,event
wset,(*pState).win2State.winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
*(*pState).win2State.xPtr = !x
*(*pState).win2State.yPtr = !y
empty
widget_control,(*pState).label,set_value = 'Finished'
widget_control,(*pState).solve,set_value = 'Solve Schrodinger Equation'
return
end
;;;;;;;;;;;;;;;;;;;
pro mdCalcSingle,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).v3Disp,get_value = v3
v3 = double(v3[0])

widget_control,(*pState).B,get_value = B
B = double(B[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

;; additive barrier parameters
widget_control,(*pState).v6,get_value = v6
v6 = double(v6[0])
;widget_control,(*pState).nsym1,get_value = nsym1
;nsym1 = fix(nsym1[0])
;widget_control,(*pState).alpha,get_value = alpha
;alpha = !dtor*double(alpha[0])

widget_control,(*pState).nstate,set_value = 0
nstate = 0
pi = !pi
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])
energylevels = dblarr(numpts)
cplxelevels = complexarr(numpts)
n = lonarr(numpts)
m = lonarr(numpts)
mlo = -(numpts-1)/2
m = mlo+dindgen(numpts)
*(*pState).mPtr = m
n = m
E3 = B   ; energy in meV for CH3 rotor
VDO = 0.0
nsym2 = 0.0
h = dblarr(numpts,numpts)

np = (1.0+bytarr(numpts))#n
mp = m#(1.0+bytarr(numpts))


;  hre = E3*(np^2)*kronecker(np,mp) + 0.25*V3*(2.0*kronecker(np,mp)-$
;               kronecker(np-mp+nsym,0)-kronecker(np-mp-nsym,0))+$
;               0.25*v6*(2.0*kronecker(np,mp)-$
;               kronecker(np-mp+6,0)-kronecker(np-mp-6,0))

  hre = E3*(np^2)*kronecker(np,mp) + 0.25*V3*(2.0*kronecker(np,mp)-$
               kronecker(np-mp+nsym,0)-kronecker(np-mp-nsym,0))+$
               0.25*v6*(2.0*kronecker(np,mp)-$
               kronecker(np-mp-6,0)+kronecker(np-mp-6,0))

;  him = -0.25*v6*sin(alpha)*(kronecker(np-mp+nsym1,0)-kronecker(np-mp-nsym1,0))
  h = hre
;  h = complex(hre,him)


trired,h,evals,e,/double
triql,evals,e,h,/double
esort = sort(evals)
evals = evals[esort]

cplxelevels[0:numpts-1] = evals
energylevels[0:numpts-1] = evals[0:numpts-1]

eigVals = energylevels[0:numpts-1]

if n_elements(*(*pState).v3Ptr) gt 2 then begin
 !x = *(*pState).win1State.xPtr
 !y = *(*pState).win1State.yPtr
 wset,(*pState).win1State.winPix
 mdPlotEig,event
 wset,(*pState).win1State.winVis
 device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
 plots,[v3,v3],!y.crange,/data

 !x = *(*pState).win2State.xPtr
 !y = *(*pState).win2State.yPtr
 wset,(*pState).win2State.winPix
 mdPlotTrans,event
 wset,(*pState).win2State.winVis
 device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
 plots,[v3,v3],!y.crange,/data
endif


number = 1.d10
roundedvalue = fltarr(numpts)
roundedvalue = round(number*(eigvals[*]))/number
whereunique = uniq(roundedvalue)
nunique = n_elements(whereunique)

noccur = intarr(nunique)
for i = 0,nunique-1 do begin
  wheresame = where(roundedvalue eq roundedvalue[whereunique[i]],ncount)
  noccur[i] = ncount
endfor
evals = eigvals[whereunique]
eSize = size(evals)
nen = eSize[1]
strout = strarr(nen+2)
strout[0] = "Eigenvalues (meV)"
strout[1] = " "
for i = 0,nen-1 do begin
     strout[i+2] = strtrim(string(eigVals[i]),2)
endfor
widget_control,(*pState).eigTextCol,set_value = strout

; Now calculate the transition values
widget_control,(*pState).jinit,get_value = initState
initialState = fix(initState[0])
strout1 = ["Transitions (meV)"," "]
for i=initialState+1,numpts-1 do begin
  thisString = strtrim(string(eigvals[i])-(eigvals[initialState]),2)
  strout1 = [strout1,thisString]
endfor
widget_control,(*pState).transTextCol,set_value = strout1
return
end
;;;;;;;;;;;;;;;;;;;
pro mdDrawEig,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).v3Ptr) lt 2 then return
widget_control,(*pState).funcgroup,get_value = val
if fix(val[0]) ne 2 then return
!x = *(*pState).win1State.xPtr
!y = *(*pState).win1State.yPtr

case event.type of
0: begin    ; button press
     (*pState).win1State.mouse = event.press
     if (*pState).win1State.mouse eq 4 then begin
       (*pState).win1State.autoscale = 1
       return
     endif
     if (*pState).win1State.mouse eq 1 then begin
       (*pState).win1State.xbox[0] = event.x
       (*pState).win1State.ybox[0] = event.y
       wset,(*pState).win1State.winVis
       device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
       empty
       ;widget_control,(*pState).win1State.win,draw_motion_events = 1
       (*pState).win1State.autoscale = 0
     endif
   end
1: begin ; button release
     xll = (*pState).win1State.xbox[0] < (*pState).win1State.xbox[1]
     yll = (*pState).win1State.ybox[0] < (*pState).win1State.ybox[1]
     w = abs((*pState).win1State.xbox[1] - (*pState).win1State.xbox[0])
     h = abs((*pState).win1State.ybox[1] - (*pState).win1State.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).win1State.xrange = [ll[0],ur[0]]
     (*pState).win1State.yrange = [ll[1],ur[1]]
     wset,(*pState).win1State.winPix
     mdPlotEig,event
     wset,(*pState).win1State.winVis
     device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
     ;widget_control,(*pState).win1State.win,draw_motion_events = 0
     (*pState).win1State.mouse = 0B
   end
2: begin ; mouse motion
     if (*pState).win1State.mouse eq 1 then begin
        (*pState).win1State.xbox[1] = event.x
        (*pState).win1State.ybox[1] = event.y
        xc = [(*pState).win1State.xbox[0],event.x,event.x,(*pState).win1State.xbox[0],$
              (*pState).win1State.xbox[0]]
        yc = [(*pState).win1State.ybox[0],(*pState).win1State.ybox[0],event.y,event.y,$
              (*pState).win1State.ybox[0]]
        wset,(*pState).win1State.winVis
        device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
        plots,xc,yc,/device
        empty
     endif else begin


       wset,(*pState).win1State.winVis
       coords = convert_coord(event.x,event.y,/device,/to_data)
       if (coords[0] gt !x.crange[0]) and (coords[0] lt !x.crange[1]) and $
       (coords[1] gt !y.crange[0]) and (coords[1] lt !y.crange[1]) then begin

         device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win1State.winPix]
         plots,[coords[0],coords[0]],!y.crange,/data
         widget_control,(*pState).v3Disp,set_value = strtrim(string(coords[0]),2)
         widget_control,(*pState).numpts,get_value = numpts
         numpts = fix(numpts[0])
         widget_control,(*pState).npts,get_value = npts
         npts = fix(npts[0])
       V3 = *(*pState).V3Ptr
       eigvals = *(*pState).EigPtr
       number = 1.d10
       roundedvalue = fltarr(numpts)
       roundedvalue = round(number*(eigvals[npts-1,*]))/number
       esize = size(eigvals)
       whereunique = indgen(esize[2]);uniq(roundedvalue)
       nunique = n_elements(whereunique)

       noccur = intarr(nunique)
       for i = 0,nunique-1 do begin
        wheresame = where(roundedvalue eq roundedvalue[whereunique[i]],ncount)
        noccur[i] = ncount
       endfor
       evals = eigvals[0:npts-1,whereunique]
       eSize = size(evals)
       nv = eSize[1]
       nen = eSize[2]

         strout = strarr(nen+2)

         strout[0] = "Eigenvalues (meV)"
           strout[1] = " "
           for i = 0,nen-1 do begin
             eig = interpol(evals[*,i],V3,coords[0])
             strout[i+2] = strtrim(string(eig),2)
           endfor
           widget_control,(*pState).eigTextCol,set_value = strout
       endif

     endelse
   end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;
pro mdDrawTrans,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).v3Ptr) lt 2 then return
!x = *(*pState).win2State.xPtr
!y = *(*pState).win2State.yPtr

case event.type of
0: begin    ; button press
     (*pState).win2State.mouse = event.press
     if (*pState).win2State.mouse eq 4 then begin
       (*pState).win2State.autoscale = 1
       return
     endif
     if (*pState).win2State.mouse eq 1 then begin
       (*pState).win2State.xbox[0] = event.x
       (*pState).win2State.ybox[0] = event.y
       wset,(*pState).win2State.winVis
       device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
       empty
       ;widget_control,(*pState).win2State.win,draw_motion_events = 1
       (*pState).win2State.autoscale = 0
     endif
   end
1: begin ; button release
     xll = (*pState).win2State.xbox[0] < (*pState).win2State.xbox[1]
     yll = (*pState).win2State.ybox[0] < (*pState).win2State.ybox[1]
     w = abs((*pState).win2State.xbox[1] - (*pState).win2State.xbox[0])
     h = abs((*pState).win2State.ybox[1] - (*pState).win2State.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).win2State.xrange = [ll[0],ur[0]]
     (*pState).win2State.yrange = [ll[1],ur[1]]
     wset,(*pState).win2State.winPix
     mdPlotTrans,event
     wset,(*pState).win2State.winVis
     device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
     ;widget_control,(*pState).win2State.win,draw_motion_events = 0
     (*pState).win2State.mouse = 0B
   end
2: begin ; mouse motion
     if (*pState).win2State.mouse eq 1 then begin
        (*pState).win2State.xbox[1] = event.x
        (*pState).win2State.ybox[1] = event.y
        xc = [(*pState).win2State.xbox[0],event.x,event.x,$
              (*pState).win2State.xbox[0],(*pState).win2State.xbox[0]]
        yc = [(*pState).win2State.ybox[0],(*pState).win2State.ybox[0],$
              event.y,event.y,(*pState).win2State.ybox[0]]
        wset,(*pState).win2State.winVis
        device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
        plots,xc,yc,/device
        empty
     endif else begin
     widget_control,(*pState).jinit,get_value = initState
     initialState = fix(initState[0])
     wset,(*pState).win2State.winPix     ; *
     mdPlotTrans,event         ; *
       wset,(*pState).win2State.winVis
       coords = convert_coord(event.x,event.y,/device,/to_data)
       if (coords[0] gt !x.crange[0]) and (coords[0] lt !x.crange[1]) and $
       (coords[1] gt !y.crange[0]) and (coords[1] lt !y.crange[1]) then begin

         device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
         plots,[coords[0],coords[0]],!y.crange,/data
         widget_control,(*pState).v3Dispa,set_value = strtrim(string(coords[0]),2)
         widget_control,(*pState).numpts,get_value = numpts
         numpts = fix(numpts[0])
         widget_control,(*pState).npts,get_value = npts
         npts = fix(npts[0])
       V3 = *(*pState).V3Ptr
       eigvals = *(*pState).EigPtr
       number = 1.d10
       roundedvalue = fltarr(numpts)
       roundedvalue = round(number*(eigvals[npts-1,*]))/number
       whereunique = uniq(roundedvalue)
       nunique = n_elements(whereunique)

       noccur = intarr(nunique)
       for i = 0,nunique-1 do begin
        wheresame = where(roundedvalue eq roundedvalue[whereunique[i]],ncount)
        noccur[i] = ncount
       endfor
       evals = eigvals[0:npts-1,whereunique]
       eSize = size(evals)
       nen = eSize[2]
         strout = ["Transitions (meV)"," "]
         for i = initialState+1,numpts-1 do begin
             trans = interpol((eigvals[0:npts-1,i])-(eigvals[0:npts-1,initialState]),$
                             V3,coords[0])
             strout = [strout,strtrim(string(trans),2)]
         endfor
         widget_control,(*pState).transTextCol,set_value = strout
       endif

     endelse

   end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;
pro mdDynamics,event
widget_control,event.top,get_uvalue = pState
psi_offset = 0.15
pot_offset = 0.15

nphi = (*pState).nphi
nframes = (*pState).nframes
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])
widget_control,(*pState).anSpeed,get_value = speed
slo = 0.0 & shi = 0.15
ds = (shi-slo)/100.0

sp = speed[0]*ds
duration = shi-sp
if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top,$
         'Perform calculation first!',/information)
  return
endif
widget_control,(*pState).jinit,get_value = jinit & jinit = fix(jinit[0])
widget_control,(*pState).jfinal,get_value = jfinal & jfinal = fix(jfinal[0])
if (jinit eq jfinal) then begin
  void = dialog_message(dialog_parent = event.top,$
         'Invalid transition!',/information)
  return
endif
if (jinit gt jfinal) then begin
  jtemp = jfinal
  jfinal = jinit
  jinit = jtemp
endif

psi_re = *(*pState).psi_RePtr
psi_im = *(*pState).psi_ImPtr
gauss_func = 0.0125*dcomplex(psi_re[jinit,0:nphi-1]+psi_re[jfinal,0:nphi-1],$
             psi_im[jinit,0:nphi-1]+psi_im[jfinal,0:nphi-1])

energylevels = *(*pState).EigPtr
eigvals = dblarr(numpts)
eigvals[0:numpts-1] = energylevels[npts-1,0:numpts-1]

if eigvals[jfinal] eq eigvals[jinit] then begin
  void = dialog_message(dialog_parent = event.top,$
         'Not an allowed transition!',/information)
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif

phi = *(*pState).PhiPtr

if abs(eigvals[jfinal]-eigvals[jinit]) gt 1.e-10 then begin
  t = dblarr(nframes) & tlo = 0.0 & thi = 2.0*!dpi/abs(eigvals[jfinal]-eigvals[jinit])
  dt = (thi-tlo)/(nframes-1)
  t = tlo+dt*dindgen(nframes)
endif else begin
  t = 0.0*dindgen(nframes)
endelse
eigfuncs = *(*pState).eigfuncs
wavepacket = dblarr(nframes,nphi)

; Here we shall assume definite transitions between eigenstates
psi_init  = dcomplexarr(nframes,nphi)
psi_final = dcomplexarr(nframes,nphi)
psi_real  = dblarr(nframes,nphi)
psi_imag  = dblarr(nframes,nphi)
einit = double(eigvals[jinit])
efinal = double(eigvals[jfinal])

if abs(eigvals[jfinal]-eigvals[jinit]) gt 1.e-10 then begin
  for i = 0,nframes-1 do begin
    psi_real[i,*] = psi_re[jinit,*]*cos(einit*t[i])+psi_im[jinit,*]*sin(einit*t[i])+$
                    psi_re[jfinal,*]*cos(efinal*t[i])+psi_im[jfinal,*]*sin(efinal*t[i])
    psi_imag[i,*] = psi_im[jinit,*]*cos(einit*t[i])-psi_re[jinit,*]*sin(einit*t[i])+$
                    psi_im[jfinal,*]*cos(efinal*t[i])-psi_re[jfinal,*]*sin(efinal*t[i])
    wavepacket[i,*] = 0.15*sqrt(psi_real[i,*]^2+psi_imag[i,*]^2)
  endfor
endif else begin
  for i = 0,nframes-1 do begin
    wavepacket[i,*] = 0.15*sqrt((psi_re[jinit,*]+psi_re[jfinal,*])^2+$
                      (psi_im[jinit,*]+psi_im[jfinal,*])^2)
  endfor
endelse

plotpot = *(*pState).PotPtr
scalefactor = 50.0*max(plotpot)/(abs(0.5*(eigvals[jfinal]+eigvals[jinit])))
thickness = 2
count = 0.0
xviewhi = 0.65
for k = 0,5 do begin
  for j = 0,nframes-1 do begin
    wugga = '!4 u !3 (rads)'
    xlabel = '!6t='+strtrim(string(count/(6.0*nframes-1.0),format = '(F6.3)'),2)
    wset,(*pState).win2State.winPix
    plot,/polar,psi_offset+wavepacket[j,*],phi[*],$
         psym=0,title = '!6Probability Time-Evolution',$
         color = (*pState).colors.white,thick = thickness,xrange=[-xviewhi,xviewhi],$
         yrange=[-xviewhi,xviewhi],xstyle=4,ystyle=4,/nodata
    oplot,/polar,psi_offset+wavepacket[j,*],phi[*],$
         psym=0,color = (*pState).colors.yellow,thick = thickness
    oplot,/polar,scalefactor*(pot_offset+plotpot),phi[*],psym=0,$
          color = (*pState).colors.green,thick = thickness
   oplot,/polar,psi_offset+0.0*findgen(nphi),phi,linestyle=2,$
         color = (*pState).colors.yellow
   oplot,/polar,scalefactor*pot_offset+0.0*findgen(nphi),phi,linestyle=2,$
         color = (*pState).colors.green
   xyouts,0.1,0.1,xlabel,/normal,color = (*pState).colors.white
    wset,(*pState).win2State.winVis
    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
    count = count+1.0
    wait,duration
  endfor
endfor

; create a pointer to the latest animated data set
dynamics = fltarr(nframes,nphi)
dynamics[0:nframes-1,0:nphi-1] = psi_offset+wavepacket[0:nframes-1,0:nphi-1]
*(*pState).dynPtr = dynamics
return
end
;;;;;;;;;;;;;;;;;;;
pro mdTransSlider,event
widget_control,event.top,get_uvalue = pState
wset,(*pState).win2State.winPix
mdPlotTrans,event
wset,(*pState).win2State.winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).win2State.winPix]
empty
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro methylDynHelp,event
widget_control,event.top,get_uvalue = pState
;pdf_file = file_which('methylcalc.pdf',/include_current_dir)
pdf_file = !DAVE_AUXILIARY_DIR+'methylcalc.pdf'
void = launch_help(pdf_file,tlb = event.top)
return
end
;;;;;;;;;;;;;;;;;;;
pro methyldyn_write_eigs,event
; Writes the eigenvalues out to a text file
widget_control,event.top,get_uvalue = pstate

if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top, $
   'Solve Schrodinger equation first',/information)
  return
endif

widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

V3 = *(*pState).V3Ptr
vmin = min(V3) & vmax = max(V3)
eigvals = *(*pState).EigPtr

widget_control,(*pState).nstate,get_value = nstate
nstate = fix(nstate[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

if nsym gt numpts then begin
  void = dialog_message('State out of range of calculation', $
   dialog_parent = event.top, /information)
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])
eigval_size = size(eigvals)

dir = (*!dave_defaults).workDir
format = '('+strtrim(string(eigval_size[2]),2)+'e15.5)'
file = dialog_pickfile(dialog_parent = event.top, path = dir, $
   filter = ['*.txt','*.TXT'],file = 'eigenvalues.txt')
if file eq "" then return
if strpos(strupcase(file),'.TXT') then file = file+'.txt'
openw,lun,file,/get_lun
printf,lun,'Barrier height,delta_E1,delta_E2,....'
printf,lun,'*************************************'
for i = 0,npts-1 do begin
   printf,lun,v3[i],eigvals[i,0:eigval_size[2]-1],format = format
endfor
free_lun,lun,/force

end
;;;;;;;;;;;;;;;;;;;
pro methyldyn_write_trans,event
; Writes the eigenvalues out to a text file
widget_control,event.top,get_uvalue = pstate


;***********
widget_control,(*pState).jinit,get_value = initState
initialState = fix(initState[0])

widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

V3 = *(*pState).V3Ptr
vmin = min(V3) & vmax = max(V3)
eigvals = *(*pState).EigPtr

widget_control,(*pState).nstate,get_value = nstate
nstate = fix(nstate[0])
widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])

if nsym gt numpts then begin
  void = dialog_message(dialog_parent = event.top,$
         'State out of range of calculation',/information)
  return
endif

eigval_size = size(eigvals)

dir = (*!dave_defaults).workDir
format = '('+strtrim(string(eigval_size[2]),2)+'e15.5)'
file = dialog_pickfile(dialog_parent = event.top, path = dir, $
   filter = ['*.txt','*.TXT'],file = 'transitions.txt')
if file eq "" then return
if strpos(strupcase(file),'.TXT') then file = file+'.txt'
openw,lun,file,/get_lun
printf,lun,'Initial state: '+strtrim(string(initialState),2)
printf,lun,'Barrier height, eig1,eig2,....'
printf,lun,'******************************'
for i = 0,npts-1 do begin
   printf,lun,v3[i],eigvals[i,0:eigval_size[2]-1]-eigvals[i,initialState], $
      format = format
endfor
free_lun,lun,/force

end
;;;;;;;;;;;;;;;;;;;
pro methyldyn_write_prob,event
; Writes the probability densities out to a text file
widget_control,event.top,get_uvalue = pstate

if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top, $
   'Solve Schrodinger equation first',/information)
  return
endif
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

widget_control,(*pState).nstate,get_value = nstate

widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])
nstate = fix(nstate[0])
eigvals = *(*pState).EigPtr
if nsym*nstate gt numpts then begin
  void = dialog_message(dialog_parent = event.top, $
   'State out of range of calculation',/information)
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

Estate = eigvals[npts-1,nstate]
widget_control,(*pState).Estate,set_value = Estate

pi = !pi
nphi = (*pState).nphi
phi = *(*pState).PhiPtr

h = *(*pState).hPtr
m = *(*pState).mPtr
evals = *(*pState).evalsPtr

wfe = fltarr(nphi)
wfo = fltarr(nphi)
prob = fltarr(nphi)
for i=0,nphi-1 do begin
  wfe[i] = ((cos(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
  wfo[i] = ((sin(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
endfor

prob = (1.0/(2.0*pi))*(wfe^2+wfo^2)+0.0*evals[nstate]
phihi = max(phi) & phimin = min(phi)
thickness = 2.0

nphi = n_elements(phi)
v3 = max(*(*pState).V3Ptr)

dir = (*!dave_defaults).workDir

file = dialog_pickfile(dialog_parent = event.top, path = dir, $
   filter = ['*.txt','*.TXT'],file = 'probability_density.txt')
if file eq "" then return
if strpos(strupcase(file),'.TXT') then file = file+'.txt'
openw,lun,file,/get_lun
printf,lun,'State index: '+strtrim(string(nstate),2)
printf,lun,'Eigenvalue: '+strtrim(string(estate),2)
printf,lun,'Barrier height: '+strtrim(string(v3),2)
printf,lun,'Angle (rads),probability density'
printf,lun,'********************************'
for i = 0,nphi-1 do begin
   printf,lun,phi[i],prob[i],format = '(2e15.4)'
endfor
free_lun,lun,/force

end
;;;;;;;;;;;;;;;;;;;
pro methyldyn_write_wf,event
; Writes the probability densities out to a text file
widget_control,event.top,get_uvalue = pstate

if n_elements(*(*pState).V3Ptr) lt 2 then begin
  void = dialog_message(dialog_parent = event.top, $
   'Solve Schrodinger equation first',/information)
  return
endif
widget_control,(*pState).numpts,get_value = numpts
numpts = fix(numpts[0])

widget_control,(*pState).nstate,get_value = nstate

widget_control,(*pState).nsym,get_value = nsym
nsym = fix(nsym[0])
nstate = fix(nstate[0])
eigvals = *(*pState).EigPtr
if nsym*nstate gt numpts then begin
  void = dialog_message(dialog_parent = event.top, $
   'State out of range of calculation',/information)
  return
endif
widget_control,(*pState).npts,get_value = npts
npts = fix(npts[0])

Estate = eigvals[npts-1,nstate]
widget_control,(*pState).Estate,set_value = Estate

pi = !pi
nphi = (*pState).nphi
phi = *(*pState).PhiPtr

h = *(*pState).hPtr
m = *(*pState).mPtr
evals = *(*pState).evalsPtr

wfe = fltarr(nphi)
wfo = fltarr(nphi)
prob = fltarr(nphi)
for i=0,nphi-1 do begin
  wfe[i] = ((cos(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
  wfo[i] = ((sin(fix(m[0:numpts-1])*phi[i])))##(transpose((h[0:numpts-1,nstate])))
endfor



nphi = n_elements(phi)
v3 = max(*(*pState).V3Ptr)

dir = (*!dave_defaults).workDir

file = dialog_pickfile(dialog_parent = event.top, path = dir, $
   filter = ['*.txt','*.TXT'],file = 'wavefunction.txt')
if file eq "" then return
if strpos(strupcase(file),'.TXT') then file = file+'.txt'
openw,lun,file,/get_lun
printf,lun,'State index: '+strtrim(string(nstate),2)
printf,lun,'Eigenvalue: '+strtrim(string(estate),2)
printf,lun,'Barrier height: '+strtrim(string(v3),2)
printf,lun,'Angle (rads),Re(psi),Im(psi)'
printf,lun,'*****************************'
for i = 0,nphi-1 do begin
   printf,lun,phi[i],wfe[i],wfo[i],format = '(3e15.4)'
endfor
free_lun,lun,/force

end
;;;;;;;;;;;;;;;;;;;
pro md_event,event
widget_control,event.top,get_uvalue = pState
if dave_set_focus(event) then return
case event.id of
(*pState).funcgroup: begin
              widget_control,event.id,get_value = val
              val = fix(val[0])
              case val of
              0:  begin
                    wset,(*pState).win1State.winPix
                  mdPlotProb,event
                  wset,(*pState).win1State.winVis
                  device,copy = [0,0,!d.x_size,!d.y_size,0,0,$
                         (*pState).win1State.winPix]
                  empty
                  widget_control,(*pState).win1State.win,$
                                  draw_motion_events = 0
                  end
              1:  begin
                    wset,(*pState).win1State.winPix
                  mdPlotWF,event
                  wset,(*pState).win1State.winVis
                  device,copy = [0,0,!d.x_size,!d.y_size,0,0,$
                         (*pState).win1State.winPix]
                  empty
                    widget_control,(*pState).win1State.win,$
                                    draw_motion_events = 0
                  end
              2:  begin
                    wset,(*pState).win1State.winPix
                  mdPlotEig,event
                  wset,(*pState).win1State.winVis
                  device,copy = [0,0,!d.x_size,!d.y_size,0,0,$
                         (*pState).win1State.winPix]
                  empty
                    widget_control,(*pState).win1State.win,$
                           draw_motion_events = 1
                  end
              endcase
            end
(*pState).nstate:    begin
               widget_control,(*pState).funcgroup,get_value = val
               val = fix(val[0])
               if val eq 1 then begin
                  wset,(*pState).win1State.winPix
                 mdPlotWF,event
                 wset,(*pState).win1State.winVis
                 device,copy = [0,0,!d.x_size,!d.y_size,0,0,$
                         (*pState).win1State.winPix]
                 empty
               endif
               if val eq 0 then begin
                  wset,(*pState).win1State.winPix
                 mdPlotProb,event
                 wset,(*pState).win1State.winVis
                 device,copy = [0,0,!d.x_size,!d.y_size,0,0,$
                         (*pState).win1State.winPix]
                 empty
               endif
            end
(*pState).v3Disp: begin
             widget_control,(*pState).v3Disp,get_value = v3
             widget_control,(*pState).v3Dispa,set_value = string(v3[0])
             mdCalcSingle,event
           end
(*pState).v3Dispa:   begin
             widget_control,(*pState).v3Dispa,get_value = v3
             widget_control,(*pState).v3Disp,set_value = string(v3[0])
             mdCalcSingle,event
           end
else:
endcase
return
end
;;;;;;;;;;;;;;;;;;;
pro aboutMethylDyn,event
strout = ['Methyl Rotor Calculator','Written by R.M.Dimeo',$
          'February 6, 2002',$
          'NIST Center for Neutron Research']
void = dialog_message(strout,/information,dialog_parent = event.top,$
                      title = 'About Methyl Rotor Calculator')
return
end
;;;;;;;;;;;;;;;;;;;
pro methyl_dyn,group_leader = group_leader
; Widget definition module
compile_opt idl2,hidden
registerName = "MD"
if xregistered(registerName) then return
!except = 0
device,decomposed = 0
tvlct,red,green,blue,/get
loadct,0,/silent
!x.margin = [12.0,3.0]

if n_elements(group_leader) ne 0 then begin
  tlb = widget_base(group_leader = group_leader,mbar = bar,$
        title = "Hindered Methyl Rotor Dynamics")
endif else begin
  tlb = widget_base(mbar = bar,$
        title = "Hindered Methyl Rotor Dynamics")
endelse
filemenu = widget_button(bar,value = "File",/menu)
;void = widget_button(filemenu,value = "Program description",event_pro = 'mdHelp')
void = widget_button(filemenu,value = "Write eigenvalues to text file",$
   event_pro = "methyldyn_write_eigs")
void = widget_button(filemenu,value = "Write transitions to text file",$
   event_pro = "methyldyn_write_trans")
void = widget_button(filemenu,value = "Write probability density to text file",$
   event_pro = "methyldyn_write_prob")
void = widget_button(filemenu,value = "Write wavefunctions (re and im) to text file",$
   event_pro = "methyldyn_write_wf")
void = widget_button(filemenu,value = "Help",event_pro = 'methylDynHelp')
void = widget_button(filemenu,value = "About",event_pro = 'aboutMethylDyn')

button_labels = ["Input Parameters","Eigenvalues","Transitions"]

; Assuming that this is IDL release 5.6 or greater, we will use the new
; built-in tab-base widget.
if (!version).release ge 5.6 then begin
  tab_base = widget_tab(tlb)
  nlabels = n_elements(button_labels)
  tab = lonarr(nlabels)
  for i = 0,nlabels-1 do begin
    tab[i] = widget_base(tab_base,title = button_labels[i],/col)
  endfor
endif else begin
   tab_base = cw_tabbase(tlb,button_labels,bases = tab,/col,/no_label,/frame)
endelse

; Build the control base
ctrlBase = widget_base(tab[0],/col,/base_align_right)
B  = cw_field(ctrlBase,title='Rotational constant (meV)',value=0.650,xsize=8,/string)
nsym  = cw_field(ctrlBase,title='Potential symmetry',value=3,xsize=8,/string)
v3hi  = cw_field(ctrlBase,title='Barrier height (meV)',value=41.0,xsize=8,/string)

;nsym1  = cw_field(ctrlBase,title='Additive potential symmetry',value=6,xsize=8,/string)
v6  = cw_field(ctrlBase,title='Additive barrier height (meV)',value=0.0,xsize=8,/string)
;alpha  = cw_field(ctrlBase,title='Additive barrier phase (degrees)',$
;                  value=0.0,xsize=8,/string)

numpts  = cw_field(ctrlBase,title='Hamiltonian size',value=51,xsize=8,/string)
npts  = cw_field(ctrlBase,title='# barrier values',value=30,xsize=8,/string)
solve = widget_button(ctrlBase,value = 'Solve Schrodinger Equation',$
                      event_pro = 'mdCalc')
label = widget_label(ctrlBase,value = 'Progress                                  ')

; Build the eigenvalue base
eigBase = widget_base(tab[1],/row)
buttonSize = 100
base = eigBase
ctrlBase1 = widget_base(base,/col,/base_align_left)
ctrlBase2 = widget_base(base,/col,/base_align_left)
plotBase = widget_base(base,/col,/base_align_center)

void = widget_button(filemenu,value = 'Exit',event_pro = 'mdQuit')

functypes = ['Probability', 'Wavefunction','Eigenvalues']
funcgroup = cw_bgroup(ctrlBase1, functypes, /col, /exclusive,$
         label_top='Selection',/no_release,$
         frame=1,set_value=0,/return_index)
nstate = widget_slider(ctrlBase1,title='Free rotor state',value = 0)
v3Disp = cw_field(ctrlBase1,title = 'V3 (meV)',xsize = 8,value = 0.0,/string,/return_events)
winxsize = 450
winysize = 400

win1 = widget_draw(base,xsize=winxsize,ysize=winysize,$
          /button_events,event_pro = 'mdDrawEig')

Estate = cw_field(ctrlBase1,title='Energy (meV)',value=0.0,xsize=8,/string)

eigTextCol = widget_text(eigBase,xsize = 20,ysize = 20,/scroll)

; Build the transitions base
transBase = widget_base(tab[2],/row)
ctrlBase2 = widget_base(transBase,/col)

void = widget_label(ctrlBase2,value='Mixture of free rotor states')
jinit = widget_slider(ctrlBase2,title='Initial state',value = 0,$
        event_pro = 'mdTransSlider')
jfinal = widget_slider(ctrlBase2,title='Final state',value = 1)
void = widget_button(ctrlBase2,value = 'Animate',$
       event_pro = 'mdDynamics',$
       xsize = buttonSize,/align_center)
anSpeed = widget_slider(ctrlBase2,title='Animation Speed (%)',$
          maximum = 100,minimum = 1,value = 100)
v3Dispa = cw_field(ctrlBase2,title = 'V3 (meV)',xsize = 8,value = 0.0,$
          /string,/return_events)
win2 = widget_draw(transBase,xsize=winxsize,ysize=winysize,$
          /button_events,/motion_events,event_pro = 'mdDrawTrans')

newbuttonSize = 100
colors = resetColors(0)

nphi = 60
transTextCol = widget_text(transBase,xsize = 20,ysize = 20,/scroll)

centertlb,tlb
widget_control,tlb,/realize

widget_control,win1,get_value = winVis1
window,/free,/pixmap,xsize = winxsize,ysize = winysize
winPix1 = !d.window

widget_control,win2,get_value = winVis2
window,/free,/pixmap,xsize = winxsize,ysize = winysize
winPix2 = !d.window


win1State = {winPix:winPix1,$
             winVis:winVis1,$
             win:win1,$
             xbox:[0.0,1.0],$
             xrange:[0.0,1.0],$
             ybox:[0.0,1.0],$
             yrange:[0.0,1.0],$
             xPtr:ptr_new(/allocate_heap),$
             yPtr:ptr_new(/allocate_heap),$
             mouse:0B,$
             autoscale:1}

win2State = {winPix:winPix2,$
             winVis:winVis2,$
             win:win2,$
             xbox:[0.0,1.0],$
             xrange:[0.0,1.0],$
             ybox:[0.0,1.0],$
             yrange:[0.0,1.0],$
             xPtr:ptr_new(/allocate_heap),$
             yPtr:ptr_new(/allocate_heap),$
             mouse:0B,$
             autoscale:1}

state = {B:B,$
         nsym:nsym,$
         solve:solve,$
         label:label,$
         v3Disp:v3Disp,$
         v3Dispa:v3Dispa,$
         anSpeed:anSpeed,$
         eigTextCol:eigTextCol,$
         transTextCol:transTextCol,$
;         nsym1:nsym1,$
;         alpha:alpha,$
         v6:v6,$
         v3hi:v3hi,$
         numpts:numpts,$
         win1State:win1State,$
         win2State:win2State,$
         red:red,green:green,blue:blue,$
         PhiPtr:ptr_new(/allocate_heap),$
         ProbPtr:ptr_new(/allocate_heap),$
         PotPtr:ptr_new(/allocate_heap),$
         V3Ptr:ptr_new(/allocate_heap),$
         EigPtr:ptr_new(/allocate_heap),$
         nstate:nstate,$
         nphi:300,$
         hPtr:ptr_new(/allocate_heap),$
         mPtr:ptr_new(/allocate_heap),$
         evalsPtr:ptr_new(/allocate_heap),$
         Estate:Estate,$
         jinit:jinit,$
         jfinal:jfinal,$
         eigfuncs:ptr_new(/allocate_heap),$
         colors:colors,$
         psi_ImPtr:ptr_new(/allocate_heap),$
         psi_RePtr:ptr_new(/allocate_heap),$
         nframes:15,$
         npts:npts,$
         funcgroup:funcgroup,$
         konstant:1,$
         dynPtr:ptr_new(/allocate_heap)}



pState = ptr_new(state,/no_copy)
widget_control,tlb,set_uvalue = pState
ret = dave_set_focus(tlb)
xmanager,registerName,tlb,cleanup = "mdCleanup",/no_block
end
;;;;;;;;;;;;;;;;;;;