; $Id: diatomic_rotor.pro,v 1.15 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  DIATOMIC_ROTOR
;
; 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.
;
;###############################################################################
; DIATOMIC_ROTOR.PRO
;
; Application that calculates the eigenvalues
; and transitions for a quantum diatomic rotor
; hindered by a 2-fold symmetric orientational
; potential.
;
; Written by R.M.Dimeo (02/21/02)
; Modification history:
;
;
;;;;;;;;;;;;;;;;;;;;;;;
pro DRCleanup,tlb
widget_control,tlb,get_uvalue = pState
wdelete,(*pState).winPix
ptr_free,(*pState).evalsPtr,(*pState).vPtr,(*pState).etPtr
drInfo = {drotorEvent,$
          id:(*pState).notifyIds[0],$
          top:(*pState).notifyIds[1],$
          handler:0l}

if widget_info((*pState).notifyIds[0],/valid_id) then begin
  widget_control,(*pState).notifyIds[0],send_event = drInfo
endif


ptr_free,pState
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRQuit,event
widget_control,event.top,/destroy
end
;;;;;;;;;;;;;;;;;;;;;;;
function kronecker,i,j
return,float(i eq j)
end
;;;;;;;;;;;;;;;;;;;;;;;
function nFact,n
; Recursive factorial function written by J.Copley
if (n gt 1) then begin
   z=double(n)*nFact(double(n-1))
endif else begin
   z=1
endelse
return,z
end
;;;;;;;;;;;;;;;;;;;;;;;
function RMDthreeJ,j1,j2,j3,m1,m2,m3
; Routine which calculates the Wigner 3-J numbers based on
; eqs (3.6.10) and (3.7.3) in Edmond's book "Angular Momentum
; in Quantum Mechanics".
;
;         / j1  j2  j3 \
; coef = |              |
;         \ m1  m2  m3 /
;
; Written by R.M.Dimeo (8/18/2000)
; Modified (02/20/02) by RMD.
;

jj1 = j1 & jj2 = j2 & jj3 = j3
mm1 = m1 & mm2 = m2 & mm3 = m3
t1 = (m1+m2+m3) ne 0
t2 = (j1+j2-j3) lt 0
t3 = (j1+j3-j2) lt 0
t4 = (j3+j2-j1) lt 0
t5 = 0;fix(J1+J2+J3) NE (J1+J2+J3)
if (t1 or t2 or t3 or t4 or t5) then return,0.D0

m1 = double(m1) & m2 = double(m2) & m3 = -1.0*double(m3)
j1 = double(j1) & j2 = double(j2) & j3 = double(j3)
term1 = kronecker(m1+m2,m3)
term2 = (2.0*j3+1)*nfact(j1+j2-j3)*nfact(j1-m1)*$
        nfact(j2-m2)*nfact(j3+m3)*nfact(j3-m3)
term2 = term2/(nfact(j1+j2+j3+1.0)*nfact(j1-j2+j3)*$
        nfact(-j1+j2+j3)*nfact(j1+m1)*nfact(j2+m2))
coef = term1*sqrt(term2)
sum = 0.0
nbig = 200
for s = 0,nbig-1 do begin
  z = double(s)
  den1 = j1-m1-z & den2 = j3-m3-z & den3 = j2-j3+m1+z
  if den1 ge 0.0 and den2 ge 0.0 and den3 ge 0.0 then begin
    num = nfact(j1+m1+z)*nfact(j2+j3-m1-z)
    den = nfact(z)*nfact(den1)*nfact(den2)*nfact(den3)
    sum = sum+((-1.0)^(z+j1-m1))*num/den
  endif
endfor
coef = coef*sum
newcoef = ((-1.0)^(j1-j2-m3))*(1.0/sqrt(2.0*j3+1))*coef
coef = newcoef
j1 = jj1 & j2 = jj2 & j3 = jj3
m1 = mm1 & m2 = mm2 & m3 = mm3
return,coef
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DREPlot,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).evalsPtr) lt 2 then return
eSize = size(*(*pState).evalsPtr)
nen = eSize[1]
nv = eSize[2]
v2 = *(*pState).vPtr
elevels = *(*pState).evalsPtr
widget_control,(*pState).B,get_value = B & B = double(B[0])
widget_control,(*pState).pltGroup,get_value = pltValue
pltValue = fix(pltValue[0])


if pltValue eq 0 then begin
  (*pState).xtitle = '!3V!D2!N (B)'
  (*pState).ytitle = '!3Eigenvalues (B)'
  (*pState).title = '!3Hindered Rotor Eigenvalues'
  if nv gt 1 then begin

    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
    if (*pState).autoscale eq 1 then begin
      dx = 0.0*(max(v2/B)-min(v2/B))
      (*pState).xrange = [min(v2/B)-dx,max(v2/B)+dx]
      dy = 0.1*(max(elevels/B))
      (*pState).yrange = [0.0,max(elevels/B)+dy]
    endif
    plot,(v2/B),elevels[0,*]/B,psym = 0,yrange = (*pState).yrange,ystyle = 1,$
       xtitle = (*pState).xtitle,xrange = (*pState).xrange,xstyle = 1,$
       ytitle = (*pState).ytitle,title = (*pState).title
    for i = 1,nen-1 do begin
      oplot,(v2/B),elevels[i,*]/B,psym = 0
    endfor
    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
  endif
endif else begin
  (*pState).xtitle = '!3V!D2!N (meV)'
  (*pState).ytitle = '!3Eigenvalues (meV)'
  (*pState).title = '!3Hindered Rotor Eigenvalues'
  if nv gt 1 then begin

    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
    if (*pState).autoscale eq 1 then begin
      dx = 0.0*(max(v2)-min(v2))
      (*pState).xrange = [min(v2)-dx,max(v2)+dx]
      dy = 0.1*(max(elevels))
      (*pState).yrange = [0.0,max(elevels)+dy]
    endif
    plot,v2,elevels[0,*],psym = 0,yrange = (*pState).yrange,ystyle = 1,$
       xtitle = (*pState).xtitle,xrange = (*pState).xrange,xstyle = 1,$
       ytitle = (*pState).ytitle,title = (*pState).title
    for i = 1,nen-1 do begin
      oplot,v2,elevels[i,*],psym = 0
    endfor
    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
  endif
endelse
(*pState).plotX = !x
(*pState).plotY = !y
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRTPlot,event
; Plot the transitions from a user-selected initial state.
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).evalsPtr) lt 2 then return
eSize = size(*(*pState).evalsPtr)
nen = eSize[1]
nv = eSize[2]
v2 = *(*pState).vPtr
elevels = *(*pState).evalsPtr
widget_control,(*pState).B,get_value = B & B = double(B[0])
widget_control,(*pState).pltGroup,get_value = pltValue
pltValue = fix(pltValue[0])
widget_control,(*pState).group,get_value = grpVal & grpVal = fix(grpVal[0])
widget_control,(*pState).tslider,get_value = slideValue
sv = fix(slideValue[0])

if pltValue eq 0 then begin
  (*pState).xtitle = '!3V!D2!N (B)'
  (*pState).ytitle = '!3Rotational Transitions (B)'
  (*pState).title = '!3Hindered Rotor Transitions'
  if nv gt 1 then begin
   ntot = nen - sv
   et = dblarr(ntot,nv)
;  print,'ntot: ',ntot
;  print,'nen: ',nen
;  print,'sv: ',sv

   ; Create a unit vector to multiply by the initial state
   ; so that we can avoid a for loop to generate the transitions
;  u = 1+bytarr(ntot)
   for i = 0,ntot-1 do begin
     et[i,0:nv-1] =  elevels[sv+i,0:nv-1] - elevels[sv,0:nv-1]
   endfor
   ;et = elevels[sv:ntot-1,0:nv-1] - (elevels[sv,0:nv-1])#u
    ;print,size(et)
   *(*pState).etPtr = et
    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
    if (*pState).autoscale eq 1 then begin
      dx = 0.0*(max(v2/B)-min(v2/B))
      (*pState).xrange = [min(v2/B)-dx,max(v2/B)+dx]
      dy = 0.1*(max(et/B))
      (*pState).yrange = [0.0,max(et/B)+dy]
    endif
    plot,(v2/B),et[0,*]/B,psym = 0,yrange = (*pState).yrange,ystyle = 1,$
       xtitle = (*pState).xtitle,xrange = (*pState).xrange,xstyle = 1,$
       ytitle = (*pState).ytitle,title = (*pState).title
    for i = 1,ntot-1 do begin
      oplot,(v2[0:nv-1]/B),et[i,0:nv-1]/B,psym = 0
    endfor
    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
  endif
endif else begin
  (*pState).xtitle = '!3V!D2!N (meV)'
  (*pState).ytitle = '!3Rotational Transitions (meV)'
  (*pState).title = '!3Hindered Rotor Transitions'
  if nv gt 1 then begin

   ntot = nen - sv
   et = dblarr(ntot,nv)
;  print,'ntot: ',ntot
;  print,'nen: ',nen
;  print,'sv: ',sv

   ; Create a unit vector to multiply by the initial state
   ; so that we can avoid a for loop to generate the transitions
   u = 1+bytarr(ntot)
;  for i = 0,ntot-1 do begin
;   et[i,0:nv-1] =    elevels[sv+i,0:nv-1] - elevels[sv,0:nv-1]
;  endfor
   et = elevels[sv:ntot-1,0:nv-1] - u#(elevels[sv,0:nv-1])
    ;print,size(et)
    *(*pState).etPtr = et

    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
    if (*pState).autoscale eq 1 then begin
      dx = 0.0*(max(v2)-min(v2))
      (*pState).xrange = [min(v2)-dx,max(v2)+dx]
      dy = 0.1*(max(et))
      (*pState).yrange = [0.0,max(et)+dy]
    endif
    plot,v2,et[0,*],psym = 0,yrange = (*pState).yrange,ystyle = 1,$
       xtitle = (*pState).xtitle,xrange = (*pState).xrange,xstyle = 1,$
       ytitle = (*pState).ytitle,title = (*pState).title
    for i = 1,ntot-1 do begin
      oplot,v2,et[i,*],psym = 0
    endfor
    !x.margin = [10.0,3.0]
    !y.margin = [4.0,2.0]
    !p.charsize = 1.0
  endif
endelse
(*pState).plotX = !x
(*pState).plotY = !y
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRCalc,event
widget_control,event.top,get_uvalue = pState

progressBar = Obj_New("SHOWPROGRESS",event.top,$
              title = "Calculation progress",message = "Thinking...")
progressBar->Start


; First get the values the user entered.
widget_control,(*pState).vlo,get_value = vlo & vlo = double(vlo[0])
widget_control,(*pState).vhi,get_value = vhi & vhi = double(vhi[0])
widget_control,(*pState).B,get_value = B & B = double(B[0])
widget_control,(*pState).nv,get_value = nv & nv = fix(nv[0])
widget_control,(*pState).jmax,get_value = jmax & jmax = fix(jmax[0])

; Now set up the Hamiltonian matrix, diagonalize it, and get the
; eigenvalues out.
if nv gt 1 then begin
  dv = (vhi-vlo)/(nv-1.0)
  v2 = vlo+dv*dindgen(nv)
endif else begin
  v2 = [vlo]
endelse
sh = (jmax+1L)*(jmax+1L)+1

elevels = dblarr(sh-1,nv)

for k = 0,nv - 1 do begin
   progressBar->Update, 100.0*(k+1)/nv
    strout = strtrim(string(k+1),2)+"/"+strtrim(string(nv),2)
    widget_control,(*pState).progress,set_value = strout
   h = dblarr(sh,sh)
   ix = 0
   for l = 0L,jmax do begin        ; l quantum number
      m = -l+indgen(2*l+1)         ; m quantum number with index j
      for j = 0L,2*l do begin
       iy = 0
      for lp = 0L,jmax do begin   ; lp quantum number
         mp = -lp+indgen(2*lp+1)     ; mp quantum number with
                       ; index jp (-lp,-(lp-1),..,0,..,(lp-1),lp)
         for jp = 0L,2*lp do begin
             l2 = 2.0 & m2 = 0.0
             wignerTerm=   (sqrt((2.0*l+1.0)*(2.0*l2+1.0)*(2.0*lp+1.0)/(4.0*!dpi)))* $
                   (rmdthreej(lp,l2,l,-mp[jp],m2,m[j]))* $
                   (rmdthreej(lp,l2,l,0,0,0))* $
                   (-1.0)^mp[jp]

             h[ix,iy]   = (B*lp*(lp+1.0)+(2.0*v2[k]/3.0))*kronecker(l,lp)* $
                     kronecker(m[j],mp[jp]) - $
                     sqrt(16.*!dpi/45.0)*v2[k]*wignerTerm
           iy = iy + 1L
         endfor
      endfor
      ix = ix + 1L
      endfor
   endfor
   h = h[0:ix-1,0:iy-1]
   trired,h,evals,e,/double
   triql,evals,e,h,/double
   esort = sort(evals)
   elevels[*,k] = evals[esort]

endfor
progressBar->Destroy
Obj_Destroy, progressBar
*(*pState).evalsPtr = elevels
*(*pState).vPtr = v2
widget_control,(*pState).group,get_value = grpVal & grpVal = fix(grpVal[0])
widget_control,(*pState).progress,set_value = 'IDLE'
widget_control,(*pState).tslider,set_slider_max = sh - 2
wset,(*pState).winPix
if grpVal eq 0 then DREPlot,event else DRTPlot,event
wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
empty
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DROneVCalc,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).group,get_value = grp & grp = fix(grp[0])
  ; First get the values the user entered.

  widget_control,(*pState).B,get_value = B & B = double(B[0])
  widget_control,(*pState).v,get_value = v & v = double(v[0])
  widget_control,(*pState).jmax,get_value = jmax & jmax = fix(jmax[0])

  ; Now set up the Hamiltonian matrix, diagonalize it, and get the
  ; eigenvalues out.

  sh = (jmax+1L)*(jmax+1L)+1

  elevels = dblarr(sh-1)

    widget_control,(*pState).progress,set_value = "Thinking...."
   h = dblarr(sh,sh)
   ix = 0
   for l = 0L,jmax do begin        ; l quantum number
      m = -l+indgen(2*l+1)         ; m quantum number with index j
      for j = 0L,2*l do begin
       iy = 0
      for lp = 0L,jmax do begin   ; lp quantum number
         mp = -lp+indgen(2*lp+1)     ; mp quantum number with
                       ; index jp (-lp,-(lp-1),..,0,..,(lp-1),lp)
         for jp = 0L,2*lp do begin
             l2 = 2.0 & m2 = 0.0
             wignerTerm=   (sqrt((2.0*l+1.0)*(2.0*l2+1.0)*(2.0*lp+1.0)/(4.0*!dpi)))* $
                   (rmdthreej(lp,l2,l,-mp[jp],m2,m[j]))* $
                   (rmdthreej(lp,l2,l,0,0,0))* $
                   (-1.0)^mp[jp]

             h[ix,iy]   = (B*lp*(lp+1.0)+(2.0*v/3.0))*kronecker(l,lp)*$
                     kronecker(m[j],mp[jp]) - $
                     sqrt(16.*!dpi/45.0)*v*wignerTerm
           iy = iy + 1L
         endfor
      endfor
      ix = ix + 1L
      endfor
   endfor
   h = h[0:ix-1,0:iy-1]
   trired,h,evals,e,/double
   triql,evals,e,h,/double
   esort = sort(evals)
   elevels[*] = evals[esort]
   nen = n_elements(elevels)

   widget_control,(*pState).pltGroup,get_value = pltype & pltype = fix(pltype[0])

   strout = strarr(nen+2)
   if (pltype eq 1) and (grp eq 0) then begin
     strout[0] = "Eigenvalues (meV)"
     strout[1] = " "
     for i = 0,nen-1 do begin
        strout[i+2] = strtrim(string(elevels[i]),2)
     endfor
   endif
   if (pltype eq 0) and (grp eq 0) then begin
     strout[0] = "Eigenvalues (B)"
     strout[1] = " "
     for i = 0,nen-1 do begin
        strout[i+2] = strtrim(string(elevels[i]/B),2)
     endfor
   endif
   widget_control,(*pState).textCol,set_value = strout
  if (n_elements(*(*pState).vPtr) gt 1) and (grp eq 0) then begin
    wset,(*pState).winVis
    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
    empty
    if pltype eq 0 then v = v/B
    coords1 = convert_coord([v,(!y.crange[0])],/data,/to_device)
    coords2 = convert_coord([v,(!y.crange[1])],/data,/to_device)
    plots,[coords1[0],coords1[0]],[coords1[1],coords2[1]],/device

  endif

; Calculate the transitions
   widget_control,(*pState).tslider,get_value = sv & sv = fix(sv[0])

   nen = n_elements(elevels)
   ntot = nen - sv
   et = dblarr(ntot)
   et = elevels[sv:ntot-1] - elevels[sv]
   ;for i = 0,ntot-1 do begin
   ; et[i] =  elevels[sv+i] - elevels[sv]
   ;endfor

   if (pltype eq 1) and (grp eq 1) then begin
     strout[0] = "Transitions (meV)"
     strout[1] = " "
     for i = 0,ntot-1 do begin
        strout[i+2] = strtrim(string(et[i]),2)
     endfor
   endif
   if (pltype eq 0) and (grp eq 1) then begin
     strout[0] = "Transitions (B)"
     strout[1] = " "
     for i = 0,ntot-1 do begin
        strout[i+2] = strtrim(string(et[i]/B),2)
     endfor
   endif
   widget_control,(*pState).textCol,set_value = strout

  if (n_elements(*(*pState).vPtr) gt 1) and (grp eq 1) then begin
    wset,(*pState).winVis
    device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
    empty
    if pltype eq 0 then v = v/B
    coords1 = convert_coord([v,(!y.crange[0])],/data,/to_device)
    coords2 = convert_coord([v,(!y.crange[1])],/data,/to_device)
    plots,[coords1[0],coords1[0]],[coords1[1],coords2[1]],/device

  endif

widget_control,(*pState).progress,set_value = 'IDLE'
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRDraw,event
widget_control,event.top,get_uvalue = pState
if n_elements(*(*pState).vPtr) lt 2 then return
widget_control,(*pState).pltGroup,get_value = pltype & pltype = fix(pltype[0])
widget_control,(*pState).B,get_value = B & B = double(B[0])
widget_control,(*pState).group,get_value = grpVal & grpVal = fix(grpVal[0])

case event.type of
0: begin    ; button press
     (*pState).mouse = event.press
     if (*pState).mouse eq 4 then begin
       (*pState).autoscale = 1
       return
     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
     endif
   end
1: begin ; button release
     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
     if grpVal eq 0 then DREPlot,event else DRTPlot,event
     wset,(*pState).winVis
     device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
     (*pState).mouse = 0B
   end
2: begin ; mouse motion
     !x = (*pState).plotX
     !y = (*pState).plotY
     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 else begin
     if grpVal eq 0 then begin
       wset,(*pState).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
         widget_control,(*pState).v,set_value = string(coords[0])

         device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
         plots,[coords[0],coords[0]],!y.crange,/data
         eSize = size(*(*pState).evalsPtr)
         nen = eSize[1]
         nv = eSize[2]
         strout = strarr(nen+2)

         if pltype eq 1 then begin
           strout[0] = "Eigenvalues (meV)"
           strout[1] = " "
           for i = 0,nen-1 do begin
             eig = interpol((*(*pState).evalsPtr)[i,*],*(*pState).vPtr,coords[0])
             strout[i+2] = strtrim(string(eig),2)
           endfor
         endif else begin
           strout[0] = "Eigenvalues (B)"
           strout[1] = " "
           for i = 0,nen-1 do begin
             eig = interpol(((*(*pState).evalsPtr)[i,*])/B,*(*pState).vPtr/B,coords[0])
             strout[i+2] = strtrim(string(eig/1.0),2)
           endfor
         endelse
         widget_control,(*pState).textCol,set_value = strout
       endif
     endif
     if grpVal eq 1 then begin
       wset,(*pState).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
         widget_control,(*pState).v,set_value = string(coords[0])

         device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
         plots,[coords[0],coords[0]],!y.crange,/data
         eSize = size(*(*pState).etPtr)
         nen = eSize[1]
         nv = eSize[2]
         strout = strarr(nen+2)

         if pltype eq 1 then begin
           strout[0] = "Transitions (meV)"
           strout[1] = " "
           for i = 0,nen-1 do begin
             eig = interpol((*(*pState).etPtr)[i,*],*(*pState).vPtr,coords[0])
             strout[i+2] = strtrim(string(eig),2)
           endfor
         endif else begin
           strout[0] = "Transitions (B)"
           strout[1] = " "
           for i = 0,nen-1 do begin
             eig = interpol(((*(*pState).etPtr)[i,*])/B,*(*pState).vPtr/B,coords[0])
             strout[i+2] = strtrim(string(eig/1.0),2)
           endfor
         endelse
         widget_control,(*pState).textCol,set_value = strout
       endif
     endif
     endelse
   end
else:
endcase

return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRTPlotIt,event
widget_control,event.top,get_uvalue = pState
widget_control,(*pState).group,set_value = 1
wset,(*pState).winPix
DRTPlot,event
wset,(*pState).winVis
device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DREqsQuit,event
;widget_control,event.top,get_uvalue = state,/no_copy
;tvlct,state.ri,state,gi,state.bi
widget_control,event.top,/destroy

return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DREqs,group_leader = group_leader
tvlct,ri,gi,bi,/get
tlb = widget_base(group_leader = group_leader,floating = 1,$
      title = 'Description of Rigid Rotor Calculation',/col,map = 0,$
      tlb_frame_attr = 1)
;filename = file_which('rotor.png',/include_current_dir)
filename = !DAVE_AUXILIARY_DIR+'rotor.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 = "DREqsQuit")

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}
widget_control,tlb,set_uvalue = state,/no_copy
xmanager,'DREqs',tlb

return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro DRDesc,event
;widget_control,event.top,sensitive = 0
DREqs,group_leader = event.top
;widget_control,event.top,sensitive = 1
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro DRHelp,event
widget_control,event.top,get_uvalue = pState
;pdf_file = file_which('diatomiccalc.pdf',/include_current_dir)
pdf_file = !DAVE_AUXILIARY_DIR+'diatomiccalc.pdf'
void = launch_help(pdf_file,tlb = event.top)
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro diatomic_rotor_event,event
if dave_set_focus(event) then return
widget_control,event.top,get_uvalue = pState

c01 = event.id eq (*pState).pltGroup
c02 = event.id eq (*pState).group
widget_control,(*pState).group,get_value = grp
grp = fix(grp[0])

if (c01 or c02) and grp eq 0 then begin
  wset,(*pState).winPix
  DREPlot,event
  wset,(*pState).winVis
  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
  empty
endif
if (c01 or c02) and grp eq 1 then begin
  wset,(*pState).winPix
  DRTPlot,event
  wset,(*pState).winVis
  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
  empty
endif

c1 = event.id eq (*pState).B
c2 = event.id eq (*pState).nv
c3 = event.id eq (*pState).vlo
c4 = event.id eq (*pState).vhi
c5 = event.id eq (*pState).jmax
if (c1 or c2 or c3 or c4 or c5) then DRCalc,event
if event.id eq (*pState).v then DROneVCalc,event

return
end
;;;;;;;;;;;;;;;;;;;
pro DRAbout,event
strout = ['Diatomic Rotor Calculator','Written by R.M.Dimeo',$
          'February 21, 2002',$
          'NIST Center for Neutron Research']
void = dialog_message(strout,/information,dialog_parent = event.top,$
       title = 'About Diatomic Rotor Calculator')
return
end
;;;;;;;;;;;;;;;;;;;;;;;
pro diatomic_rotor,group_leader = group_leader,notifyIds = notifyIds
; Widget definition module
if n_elements(notifyIds) eq 0 then return
loadct,0,/silent
if n_elements(group_leader) gt 0 then begin
  tlb = widget_base(/row,mbar = bar,title = "Hindered Diatomic Quantum Rotor Calculator")
endif else begin
  tlb = widget_base(group_leader = group_leader,/row,mbar = bar,$
        title = "Hindered Diatomic Quantum Rotor Calculator")
endelse
filemenu = widget_button(bar,value = "File")
;void = widget_button(filemenu,value = "Program Description",event_pro = "DRDesc")
void = widget_button(filemenu,value = "Help",event_pro = "DRHelp")
void = widget_button(filemenu,value = "About",event_pro = "DRAbout")
quit = widget_button(filemenu,value = "Quit",event_pro = "DRQuit")

base1 = widget_base(tlb,/col,/base_align_left,/frame)
text1Size = 5
butSize = 75
void = widget_label(base1,value = 'B (rotational',/align_left)
void = widget_label(base1,value = 'constant, meV)',/align_left)
B = cw_field(base1,title = ' ',value = '1.0',/string,$
    /return_events,xsize = text1Size)
void = widget_label(base1,value = 'VLO (meV)',/align_left)
vlo = cw_field(base1,title = ' ',value = '0.0',/string,$
      /return_events,xsize = text1Size)
void = widget_label(base1,value = 'VHI (meV)',/align_left)
vhi = cw_field(base1,title = ' ',value = '10.0',/string,$
      /return_events,xsize = text1Size)
void = widget_label(base1,value = 'NV',/align_left)
nv = cw_field(base1,title = ' ',value = '5',/string,$
     /return_events,xsize = text1Size)
void = widget_label(base1,value = 'Lmax',/align_left)
jmax = cw_field(base1,title = ' ',value = '6',/string,$
       /return_events,xsize = text1Size)
base1a = widget_base(base1,/col,/base_align_bottom)

base2 = widget_base(tlb,/col,/base_align_left,/frame)


text2Size = 12
void = widget_label(base2,value = 'STATUS',/align_center)
progress = cw_field(base2,title = ' ',value = 'IDLE',$
           /string,/return_events,xsize = text2Size)

pltypes = ['Reduced (B)','E (meV)']
pltGroup = cw_bgroup(base2, pltypes, /col, /exclusive,$
         label_top='PLOT UNITS',/no_release,$
         frame=1,set_value=0,/return_index)

void = widget_label(base2,value = 'V (meV)',/align_center)
v = cw_field(base2,title = ' ',value = ' ',$
           /string,/return_events,xsize = text2Size)
types = ['Eigenvalues','Transitions']
group = cw_bgroup(base2, types, /col, /exclusive,$
         label_top='PLOT TYPE',/no_release,$
         frame=1,set_value=0,/return_index)

tslider = widget_slider(base2,value = 0,maximum = 10,minimum = 0,$
          event_pro = 'DRTPlotIt',title = 'Initial eigenstate')

calc = widget_button(base1a,value = "CALCULATE",$
       event_pro = "DRCalc",xsize = butSize)


xsize = 400 & ysize = 400
win = widget_draw(tlb,xsize = xsize,ysize = ysize,$
      /button_events,/motion_events,event_pro = "DRDraw")

textCol = widget_text(tlb,xsize = 20,ysize = 20,/scroll)
centertlb,tlb
widget_control,tlb,/realize
widget_control,win,get_value = winVis
window,/free,/pixmap,xsize = xsize,ysize = ysize
winPix = !d.window

state = {vlo:vlo,$
         vhi:vhi,$
         nv:nv,$
         B:B,$
         v:v,$
         progress:progress,$
         textCol:textCol,$
         jmax:jmax,$
         winPix:winPix,$
         winVis:winVis,$
         win:win,$
         group:group,$
         tslider:tslider,$
         pltGroup:pltGroup,$
         vPtr:ptr_new(/allocate_heap),$
         etPtr:ptr_new(/allocate_heap),$
         evalsPtr:ptr_new(/allocate_heap),$
         xtitle:'!3V!D2!N (B)',$
         ytitle:'!3Eigenvalues (B)',$
         title:'!3Hindered Rotor Eigenvalues',$
         autoscale:1,$
         mouse:0B,$
         xbox:[0.0,1.0],$
         ybox:[0.0,1.0],$
         xrange:[0.0,1.0],$
         yrange:[0.0,1.0],$
         plotX:!x,$
         plotY:!y,$
         notifyIds:notifyIds}

pState = ptr_new(state,/no_copy)
widget_control,tlb,set_uvalue = pState
ret = dave_set_focus(tlb)
xmanager,'diatomic_rotor',tlb,cleanup = 'DRCleanup',/no_block
end