; $Id: methylse.pro,v 1.3 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  METHYLSE
;
; PURPOSE:
;  Function that solves the 3-fold rotor Schrodinger equation
;
; 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.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;
function methylSE,numpts = numpts,$
                  B = B,$
                  v3 = v3,$
                  nsym1 = nsym1,$
                  nsym2 = nsym2,$
                  v6 = v6,$
                  k = k
; This function solves the Schrodinger equation for a
; particular value of v3
; numpts:	size of each dimension of Hamiltonian array
; B:		rotational constant (meV) of rotor
; V3:		barrier height (meV)
; nsym:		symmetry number of hindering barrier
; nsym2:	symmetry number of additive barrier
; v6:		strength of additive barrier
; k:		sign for additive barrier
if n_elements(numpts) eq 0 then numpts = 51
if n_elements(B) eq 0 then B = 0.650
if n_elements(v3) eq 0 then v3 = 41.0
if n_elements(nsym1) eq 0 then nsym1 = 3
if n_elements(v6) eq 0 then v6 = 0.0
if n_elements(nsym2) eq 0 then nsym2 = 6
if n_elements(k) eq 0 then k = 1.0

m = lonarr(numpts)
mlo = -(numpts-1)/2
m = mlo+dindgen(numpts)
np = (1.0+bytarr(numpts))#m
mp = m#(1.0+bytarr(numpts))
h = B*(np^2)*kronecker(np,mp) + 0.25*v3*(2.0*kronecker(np,mp)-$
              kronecker(np-mp+nsym1,0)-kronecker(np-mp-nsym1,0))+$
              0.25*v6*(2.0*kronecker(np,mp)+$
              ((-1.0)^k)*(kronecker(np-mp+nsym2,0)+kronecker(np-mp-nsym2,0)))

trired,h,evals,e,/double
triql,evals,e,h,/double
esort = sort(evals)
evals = evals[sort(evals)]
return,evals[0:numpts-1]
end