; $Id: spurion.pro,v 1.14 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  SPURION
;
; 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.
;
;###############################################################################
; SPURION.PRO
;
; This application calculates the possible elastic incoherent scattering
; processes as well as the inelastic scattering processes of higher
; harmonics in a TAS measurement.
;
; Written by R.M.Dimeo (02/08/02)
;
;;;;;;;;;;;;;
;
pro eqCleanup,tlb
widget_control,tlb,get_uvalue = state,/no_copy
; In order to make sure that the coordinate system is
; correct when we exit this procedure we have to restore
; the system variables !x and !y which define the
; coordinate system.
wset,(*state.upState).winVis
!x = *(*state.upState).xPtr
!y = *(*state.upState).yPtr
!x.margin = [10.,3.]
!y.margin = [4.,2.]
return
end
;;;;;;;;;;;;;;
pro spurionEqDispQuit,event
widget_control,event.top,/destroy
return
end
;;;;;;;;;;;;;;
pro spurionEqDisp,event
; This widget program consists of a plot window and a
; dismiss button.  The plot window displays the equations
; that are used in the calculations here.
widget_control,event.top,get_uvalue = upState
tlb = widget_base(group_leader = event.top,/modal,$
      title = 'Spurion equations',/col)
;filename = file_which('speqs.png',/include_current_dir)
filename = !DAVE_AUXILIARY_DIR+'speqs.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 = "spurionEqDispQuit")

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

state = {upState:upState}
widget_control,tlb,set_uvalue = state,/no_copy
xmanager,'spurionEqDisp',tlb,cleanup = "eqCleanup"
return
end
;;;;;;;;;;;;;;
pro spurionCleanup,tlb
widget_control,tlb,get_uvalue = pState,/no_copy

s = size((*pState).notifyIDs)
if n_elements((*pState).notifyIDs) gt 1 then begin
  if s[0] eq 1 then count = 0 else count = s[2] - 1
  for j = 0,count do begin
    pseudoEvent = {spurionEvent,$
                   ID:(*pState).notifyIDs[0,j],$
                   top:(*pState).notifyIDs[1,j],$
                   handler:0l}
    if widget_info((*pState).notifyIDs[0,j],/valid_id) then $
       widget_control,(*pState).notifyIDs[0,j],send_event = pseudoEvent,/no_copy
  endfor
endif


tvlct,(*pState).r,(*pState).g,(*pState).b
!x.margin = [10.0,3.0]
!y.margin = [4.0,2.0]
wdelete,(*pState).winPix
ptr_free,(*pState).ePtr,(*pState).dePtr,(*pState).eaPtr
ptr_free,(*pState).xPtr,(*pState).yPtr
ptr_free,pState
return
end
;;;;;;;;;;;;;;
pro spurionQuit,event
widget_control,event.top,/destroy
return
end
;;;;;;;;;;;;;;
pro spurionPlot,event
widget_control,event.top,get_uvalue = pState,/no_copy
if n_elements(*(*pState).dePtr) gt 2 then begin

  ymax = 30.0
  ymin = 0.0
  de = *(*pState).dePtr
  e = *(*pState).ePtr
  s = size(de)
  n = s[1] & npts = s[3]
  if (*pState).autoscale eq 1 then begin
    ymax = max(abs(de[0:n-2,0:n-2,*])) & ymin = min(abs(de[0:n-2,0:n-2,*]))
    xmin = min(e) & xmax = max(e)
    dy = ymax-ymin & dx = xmax-xmin
    (*pState).xviewrange = [xmin,xmax]
    (*pState).yviewrange = [0.0,200.0]
  endif

  wset,(*pState).winPix
  for i = 0,n-1 do begin
    for j = 0,n-1 do begin
      if i eq 0 and j eq 0 then begin
        plot,e,de[i,j,*],psym = 0,yrange = (*pState).yviewrange,ystyle = 1,$
        ytitle = (*pState).ylabel,charsize = 1.5,xtitle = (*pState).xlabel,$
        title = (*pState).title,xrange = (*pState).xviewrange,xstyle = 1
      endif
      if total(de[i,j,*]) gt 0 then begin
        oplot,e,de[i,j,*],linestyle = 0
      endif else begin
        oplot,e,abs(de[i,j,*]),linestyle = 2
      endelse
    endfor
  endfor
  wset,(*pState).winVis
  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
endif
*(*pState).xPtr = !x
*(*pState).yPtr = !y
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;
pro spurionCalc,event
widget_control,event.top,get_uvalue = pState,/no_copy
widget_control,(*pState).group,get_value = value
widget_control,(*pState).order,get_value = order

n = fix(order[0]) ; upper limit on the order of beam contamination
nm = 1 + indgen(n)
na = 1 + indgen(n)
nap = fltarr(n,n)
uvec = 1+bytarr(n)

(*pState).n = n

if value[0] eq 1 then begin   ; Ef is fixed
  (*pState).xlabel = '!3E!Dfinal!n (meV)'
  for i = 0,n-1 do begin
    nap[i,0:n-1] = (shift(na,n))[i]
  endfor
  nanm = (nap/(uvec#nm))^2-1.0
  npts = 100
  elo = 0.0 & ehi = 100.0
  ef = hfbs_makepoints(xlo = elo,xhi = ehi,npts = npts)
  de = fltarr(n,n,npts)
  for i = 0,npts-1 do begin
    de[*,*,i] = nanm[*,*]*ef[i]
  endfor
  *(*pState).ePtr = ef
  *(*pState).dePtr = de
endif
if value[0] eq 0 then begin   ; Ei is fixed
  (*pState).xlabel = '!3E!Dinitial!n (meV)'
  for i = 0,n-1 do begin
    nap[i,0:n-1] = (shift(na,n))[i]
  endfor
  nanm = (nap/(uvec#nm))^2-1.0
  nanm = -1.0*transpose(temporary(nanm))
  npts = 100
  elo = 0.0 & ehi = 100.0
  ei = hfbs_makepoints(xlo = elo,xhi = ehi,npts = npts)
  de = fltarr(n,n,npts)
  for i = 0,npts-1 do begin
    de[*,*,i] = nanm[*,*]*ei[i]
  endfor
  *(*pState).ePtr = ei
  *(*pState).dePtr = de
endif
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;
pro spurionDraw,event
widget_control,event.top,get_uvalue = pState,/no_copy
if n_elements(*(*pState).dePtr) lt 2 then begin
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif
widget_control,(*pState).group,get_value = grp & grp = fix(grp[0])
case event.type of
0: begin    ; button press
     (*pState).mouse = event.press
     if (*pState).mouse eq 4 then begin
       (*pState).autoscale = 1
     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).xviewrange = [ll[0],ur[0]]
     (*pState).yviewrange = [ll[1],ur[1]]
     wset,(*pState).winPix
     widget_control,event.top,set_uvalue = pState,/no_copy
     spurionPlot,event
     widget_control,event.top,get_uvalue = pState,/no_copy
     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
     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
       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).efout,set_value = string(coords[0])
         widget_control,(*pState).deout,set_value = string(coords[1])

         device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
           n = (*pState).n
           nn = n*n
           stroutEL = strarr(nn+3)
           stroutIE = strarr(nn+4)
           de = *(*pState).dePtr
           e = *(*pState).ePtr
           stroutIE[0] = "Inelastic"
           stroutIE[1] = ""
           stroutEL[0] = "Elastic incoherent"
           stroutEL[1] = ""
           stroutEL[2] = "(mon,an): dEapp"
      count = 3
      for i = 0,n-1 do begin
        for j = 0,n-1 do begin
            nm = i+1 & na = j+1
            if grp eq 1 then begin
              eact = (nm^2-na^2)*coords[0]+nm^2*coords[1]
              stroutIE[2] = "(mon,an): dEact"
              stroutIE[3] = ""
           endif else begin
             eact = (nm^2-na^2)*coords[0]+na^2*coords[1]
             stroutIE[2] = "(mon,an): dEact"
           endelse
           yout = interpol(de[i,j,*],e[*],coords[0])
           stroutIE[count+1] = "("+strcompress(string(i+1))+","+$
                             strcompress(string(j+1))+"): "+$
                             strcompress(string(eact,format = fcode))
           if i ne j then begin
             stroutEL[count] = "("+strcompress(string(j+1))+","+$
                             strcompress(string(i+1))+"): "+$
                             strcompress(string(yout,format = fcode))
           endif else begin
             stroutEL[count] = ""
           endelse
          count = count + 1
        endfor
      endfor
      widget_control,(*pState).elText,set_value = stroutEL
      widget_control,(*pState).ieText,set_value = stroutIE
       endif

     endelse
   end
else:
endcase
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;
pro spurion_event,event
ret = dave_set_focus(event)
widget_control,event.top,get_uvalue = pState,/no_copy
if (event.id eq (*pState).deout) or (event.id eq (*pState).efout) then begin
  ; calculate the exact spurions
  widget_control,(*pState).deout,get_value = deout & deout = float(deout[0])
  widget_control,(*pState).efout,get_value = efout & efout = float(efout[0])
  n = (*pState).n
  widget_control,(*pState).group,get_value = grp & grp = fix(grp[0])
  coords = [efout,deout]

  wset,(*pState).winVis
  device,copy = [0,0,!d.x_size,!d.y_size,0,0,(*pState).winPix]
  plots,[coords[0],coords[0]],!y.crange,linestyle = 0,thick = 1.5,/data
  plots,!x.crange,[coords[1],coords[1]],linestyle = 0,thick = 1.5,/data
  n = (*pState).n
  nn = n*n
  stroutEL = strarr(nn+3)
  stroutIE = strarr(nn+4)
  de = *(*pState).dePtr
  e = *(*pState).ePtr
  stroutIE[0] = "Inelastic"
  stroutEL[0] = "Elastic incoherent"
  stroutEL[1] = ""
  stroutEL[2] = "(mon,an): dEapp"
  count = 3
  for i = 0,n-1 do begin
   for j = 0,n-1 do begin
     nm = i+1 & na = j+1
     if grp eq 1 then begin
     eact = (nm^2-na^2)*efout+(nm^2)*deout
     stroutIE[1] = ""
       stroutIE[2] = "(mon,an): dEact"
       stroutIE[3] = ""
     endif else begin
       eact = (nm^2-na^2)*efout+na^2*deout
       stroutIE[1] = ""
       stroutIE[2] = "(mon,an): dEact"
       stroutIE[3] = ""
     endelse
     yout = interpol(de[i,j,*],e[*],coords[0])
     stroutIE[count+1] = "("+strcompress(string(nm))+","+$
                             strcompress(string(na))+"): "+$
                             strcompress(string(eact,format = fcode))
     if i ne j then begin
       stroutEL[count] = "("+strcompress(string(na))+","+$
                             strcompress(string(nm))+"): "+$
                             strcompress(string(yout,format = fcode))
     endif else begin
     stroutEL[count] = ""
     endelse
     count = count + 1
   endfor
endfor
widget_control,(*pState).elText,set_value = stroutEL
widget_control,(*pState).ieText,set_value = stroutIE
endif
if (event.id eq (*pState).group) or (event.id eq (*pState).order) then begin
  widget_control,(*pState).order,get_value = n
  n = fix(n[0])
  if n le 2 then begin
    out = 'n (order) > 2'
    void = dialog_message(dialog_parent = event.top,out)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
  widget_control,event.top,set_uvalue = pState,/no_copy
  spurionCalc,event
  spurionPlot,event
  widget_control,event.top,get_uvalue = pState,/no_copy
endif

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;
pro spurionDesc,event
strout   = ["This program calculates the apparent inelastic scattering",$
       "due to the elastic and inelastic scattering of higher order",$
       "harmonics from the monochromator and analyzer.  Drag the cursor",$
       "around the screen to determine which harmonics result in apparent",$
       "inelastic features.  The dashed lines indicate a negative apparent",$
       "energy transfer.  The order of the transmission is displayed in the",$
       "text box at the right as you drag the cursor around the screen."]
void = dialog_message(dialog_parent = event.top,strout,/information)
return
end
;;;;;;;;;;;;;;
pro spurionAbout,event
strout   = ["Spurion Calculator",$
       "Written by R.M.Dimeo (02/07/02)",$
       "NIST Center for Neutron Research"]
void = dialog_message(dialog_parent = event.top,strout,/information)
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro spurionHelp,event
widget_control,event.top,get_uvalue = pState
;pdf_file = file_which('spurcalc.pdf',/include_current_dir)
pdf_file = !DAVE_AUXILIARY_DIR+'spurcalc.pdf'
void = launch_help(pdf_file,tlb = event.top)
return
end
;;;;;;;;;;;;;;
pro spurion,NOTIFYIDS = notifyIDs,INPTR = inPtr
; Widget definition module
tvlct,r,g,b,/get
loadct,0,/silent

if not keyword_set(notifyIDs) then notifyIDs = [0L,0L]
if not keyword_set(inPtr) then inPtr = 0.0
if notifyIds[1] eq 0L then group_leader = 0L else group_leader = notifyIDS[1]
!y.margin = [5.0,3.0]
tlb = widget_base(/col,title = 'Spurion calculator',/tlb_frame_attr,mbar = bar, $
   group_leader = group_leader)
filemenu = widget_button(bar,value = 'File')
miscmenu = widget_button(bar,value = 'Misc')
void = widget_button(miscmenu,value = 'Help',event_pro = 'spurionHelp')
void = widget_button(miscmenu,value = 'About spurion calculator',event_pro = 'spurionAbout')

rowBase1 = widget_base(tlb,/row)
type = ['Ei','Ef']
group = cw_bgroup(rowBase1,type, /row, /exclusive,$
         label_top='Fixed energy',/no_release,$
         frame=1,set_value=1,/return_index)
order = cw_field(rowBase1,value = '4',/string,xsize = 5,title = 'Max order',$
        /return_events,/col)
efout = cw_field(rowBase1,value = "",/string,title = 'Fixed energy',xsize = 6,$
        /col,/return_events)
deout = cw_field(rowBase1,value = "",/string,title = 'Apparent energy transfer (meV)',$
        xsize = 6,/col,/return_events)

void = widget_button(filemenu,value = 'Quit',event_pro = 'spurionQuit')
rowBase2 = widget_base(tlb,/row)
winxsize = 550 & winysize = 400
plotWin = widget_draw(rowBase2,xsize = winxsize,ysize = winysize,$
          event_pro = 'spurionDraw',/button_events,/motion_events)
elText = widget_text(rowBase2,xsize = 15,ysize = 20,/scroll)
ieText = widget_text(rowBase2,xsize = 15,ysize = 20,/scroll)
centertlb,tlb
widget_control,tlb,/realize

window,/free,/pixmap,xsize = winxsize,ysize = winysize
winPix = !d.window
widget_control,plotWin,get_value = winVis

state = {group:group,$
         plotWin:plotWin,$
         winPix:winPix,$
         winVis:winVis,$
         winxsize:winxsize,$
         winysize:winysize,$
         elText:elText,$
         ieText:ieText,$
         ePtr:ptr_new(/allocate_heap),$
         dePtr:ptr_new(/allocate_heap),$
         ylabel:"!4D!3E!Dapparent!N (meV)",$
         xlabel:"",$
         title:"False Peak Conditions",$
         order:order,$
         autoscale:1,$
         fcode:"(f6.2)",$
         xviewrange:[0.0,1.0],$
         yviewrange:[0.0,1.0],$
         xbox:[0,1],$
         ybox:[0,1],$
         mouse:0B,$,
         r:r,g:g,b:b,$
         efout:efout,$
         deout:deout,$
         n:0,$
         eaPtr:ptr_new(/allocate_heap),$
         notifyIDs:notifyIDs,$
         inPtr:inPtr,$
         xPtr:ptr_new(/allocate_heap),$
         yPtr:ptr_new(/allocate_heap)}

widget_control,tlb,set_uvalue = ptr_new(state,/no_copy)
spurionCalc,{event,id:elText,top:tlb,handler:0l}
spurionPlot,{event,id:ieText,top:tlb,handler:0l}
ret = dave_set_focus(tlb)
xmanager,'spurion',tlb,cleanup = 'spurionCleanup'
end