; $Id: resolutiongui.pro,v 1.2 2005/02/09 14:02:27 azuah Exp $
;###############################################################################
; 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 or if the code in this file is
;  included in another product.
;
;###############################################################################
function twoslitfwhm,w1,w2,l
;****************************************************************
;* return  effective gaussian FWHM corresponding to a two slit  *
;* collimator (in radians)            *
;****************************************************************

w1  =    max([w1,l*1.e-6])
w2  =    max([w2,l*1.e-6])

return,sqrt(alog(2.)/12/l^2/w1/w2*((w1+w2)^4-(w1-w2)^4))

end






pro fwhmellips,m,x,y
;****************************************************************
;*                    *
;* return x and y coordinates for which :      *
;*  exp(-1/2*(m(0,0)*x^2+m(1,1)*y^2+2*m(0,1)*x*y))=0.5   *
;* modified 8/13/97 to be more robust. Assumes positive definite*
;* matrix else we do not have an ellips          *
;*                    *
;****************************************************************


disk    =  (m(0,0)+m(1,1))^2-4.*(m(0,0)*m(1,1)-m(0,1)^2)
if m(0,1) ne 0 then begin
    lam0   = (m(0,0)+m(1,1)+sqrt(disk))/2.
    lam1   = (m(0,0)+m(1,1)-sqrt(disk))/2.
    tanv0  =    (lam0-m(0,0))/m(0,1)
    tanv1  =    (lam1-m(0,0))/m(0,1)
    v0 =   [1,tanv0]/sqrt(1.+tanv0^2)
    v1 =   [1,tanv1]/sqrt(1.+tanv1^2)
end else begin
    lam0   = m(0,0)
    lam1   = m(1,1)
    v0 =   [1,0]
    v1 =   [0,1]
end

tln2    =  2.*alog(2.)

v   = findgen(200)/199.*2.01*!pi
xt  =    sqrt(tln2/abs(lam0))*cos(v)
yt  =    sqrt(tln2/abs(lam1))*sin(v)

x   = xt*v0(0)+yt*v1(0)
y   = xt*v0(1)+yt*v1(1)


end






pro proj_scatplane,alfa1,dq1,dq2, event
;************************************************************************
;*                      *
;* returns the projection of the resolution function defined by rm  *
;* on the scattering  plane  in terms of the principal lengths dq1 and  *
;* and dq2 and the angle alfa1 clockwise from scattering vector to  *
;* axis of ellipse with length dq1            *
;*                      *
;* Formula taken from the Brookhaven version of the reselps program *
;* By Collin Broholm 8/17/92                 *
;*                      *
;************************************************************************
;common resout,r0,rm
widget_control, event.top, get_uvalue = ptrAll
      rm=(*ptrAll).rmatrix
      fwhmt1e = 1./sqrt(8.*alog(2.))
      ap=rm(0,0)-rm(0,2)*rm(0,2)/rm(2,2)
      b=rm(1,1)-rm(1,2)*rm(1,2)/rm(2,2)
      c=rm(0,2)*rm(1,2)/rm(2,2)-rm(0,1)
      v2=atan(2*c/(b-ap))
      v=v2/2
      dq1=1/sqrt(ap*cos(v)^2+b*sin(v)^2-c*sin(v2))/fwhmt1e
      dq2=1/sqrt(ap*sin(v)^2+b*cos(v)^2+c*sin(v2))/fwhmt1e
      alfa1=v

end



pro tasin,file
;common resinp,tm,ta,epm,ep,etam,etaa,alz,alm,ala,al3,vbetz,vbet1,vbet2,vbet3,$
;              etash,etasv,efix,feief,source
;     list of symbols
;     aki incident neutron wavevector  (aa-1)
;     tm monochromator tau (aa-1)
;     ta analyzer tau (aa-1)
;     epm +1. if sample scattering direction opposite of monochr. sc. di
;     ep +1. if analyzer scatt.dirct. same as monochr.sc.dirct.
;     if ep or epm is not +1, it is -1.
;     alz etc. horizontal divergences fwhm in minutes
;     vbet1 etc. vertical divergences fwhm in minutes
;     for ei fixed ein is positive.
;     for ef fixed ein is negative.  ef=abs(ein)

print,'Note that input file has been modified to require source specification'
print,'after the horizontal AND vertical sample mosaic:'
print,'0 for cold source, 1 for thermal source'
openr,lun,file,/get_lun
      readf,lun,epm,ep,alz,alm,ala,al3,tm,ta
      readf,lun,etam,etaa,efix,etash,etasv,source
      readf,lun,vbetz,vbet1,vbet2,vbet3
free_lun,lun

;     **********************************
;     *figure out fixed ei or fixed ef *
;     **********************************
      if (efix gt 0) then feief =   0 else feief = 1
      efix  =    abs(efix)

;     converts fwhm minutes to 1/e radians
      degtrad=3.14159/180./sqrt(8.*alog(2.))/60.
      alz=alz*degtrad
      etam=etam*degtrad
      alm=alm*degtrad
      ala=ala*degtrad
      etaa=etaa*degtrad
      al3=al3*degtrad
      vbetz=vbetz*degtrad
      vbet1=vbet1*degtrad
      vbet2=vbet2*degtrad
      vbet3=vbet3*degtrad
      etash=degtrad*etash
      etasv=degtrad*etasv


      end


function angle,x,y
;****************************************************************
;*                    *
;* regarding x,y as the coordinates of a vector, angle returns  *
;* the angle from the x-axis to that vector          *
;*                    *
;****************************************************************
return,atan(y,x)
end

function monc,emon, event
;****************************************************************
;* monitor correction factor for BT4 at NIST can be used as *
;* typical result for thermal neutron source although best  *
;* obtain actual result for source and monochromator in question*
;****************************************************************
;common resinp,tm,ta,epm,ep,etam,etaa,alz,alm,ala,al3,vbetz,vbet1,vbet2,vbet3,$
;              etash,etasv,efix,feief,source
widget_control, event.top, get_uvalue = ptrAll
feief=(*ptrAll).feief
if feief eq 0 then return,1.    ;fixed Ei : no monitor correction
if source eq 0 then return,1.   ;cold source, no correction yet

;----------------------------------------------------
;| Thermal source use BT4 measurement by D. Neumann |
;----------------------------------------------------
if emon lt 12  then $
  print,'monitor correction factor out of calibration range'

result  =.0479*(33.3+32.87/(emon-7.57)-.66*emon+.00857*emon^2)

return,result
end

function pa,e
;************************************************************************
;*                      *
;* reflectivity of analyzer as a function of energy and as defined  *
;* in the Chesser abd Axe paper. This function is only weakly dependent *
;* on energy decreasing slightly. The form interpolated here was    *
;* measured at BT4 in January 1992 as ratio of intensity of Si(111) *
;* powder peak measure with 20'-40'  collimtion before and after    *
;* PG(002) analyzer and with just 20' collimation before analyzer   *
;*                      *
;************************************************************************

;* for the time being fix to 1 *
return,1

end

pro phonint,result, event
;***************************************************************
;* calculate energy integrated intensity of planar dispersion  *
;* surface normalized to monitor counts. This is formula       *
;* (16) of Chesser & Axe taking monitor sensistivity to be     *
;* proportional to monc(ei)/ki                 *
;* resolution function is assumed previously calculated and    *
;* present in variables r0 and rm of common block resout       *
;* error corrected 3/15/92 Collin Broholm            *
;***************************************************************
;common resout,r0,rm
widget_control, event.top, get_uvalue = ptrAll
rm=(*ptrAll).rm
r0=(*ptrAll).r0
result  =    r0*(2*!pi)^2/sqrt(abs(determ(rm)))

end


pro oneddisp,c,wq,we,odrm, event
;****************************************************************
;* calculate fwhm wq of constant e scan and fwhm we of   *
;* constant q scan for 1D dispersion relation with gradient *
;* c (meV*AA) along the first coordinate direction of the   *
;* resolution matrix currently saved in resout       *
;* Also return the 2x2 matrix characterizing resolution effects *
;* in the e-q plane of dispersion          *
;*                    *
;* formulae modified 12/11/92             *
;* formula for odrm was incorrect until 8/13/97 when it was     *
;*  corrected by Collin Broholm from       *
;*  2*[[alfa,gama],[gama,beta]]  to          *
;*  [[2.*alfa,gama],[gama,2.*beta]]          *
;*                    *
;****************************************************************
;common resout,r0,rm
widget_control, event.top, get_uvalue = ptrAll
rm=(*ptrAll).rm
ln2 =   alog(2.)
alfa    =  0.5*(rm(0,0)-rm(1,0)*rm(1,0)/rm(1,1))
beta    =  0.5*(rm(2,2)-rm(2,1)*rm(2,1)/rm(1,1))
gama    =  rm(2,0)-rm(1,0)*rm(2,1)/rm(1,1)


we  =    2.*sqrt(ln2*4.*(alfa+beta*c^2+gama*c)/(4.*alfa*beta-gama^2))
wq  =    we/c
odrm    =  [[2.*alfa,gama],[gama,2.*beta]]
end

pro eintqres,m, event
;*******************************************************************
;* calculate the 2x2 q resolution matrix characterizing the energy *
;* transfer integrated wave vector resolution         *
;* created 8/13/97 by Collin Broholm             *
;*******************************************************************
;common resout,r0,rm
widget_control, event.top, get_uvalue = ptrAll
rm=(*ptrAll).rm
m   = fltarr(2,2)
m(0,0)  =    rm(0,0)-rm(0,2)^2/rm(2,2)
m(1,1)  =    rm(1,1)-rm(1,2)^2/rm(2,2)
m(0,1)  =    rm(0,1)-rm(0,2)*rm(1,2)/rm(2,2)
m(1,0)  =    m(0,1)

end


pro reselps,event
;************************************************************************
;* PURPOSE :                    *
;*  CALCULATE THE NORMALIZATION CONSTANT AND CHARACTERISTIC MATRIX  *
;*  OF THE GAUSSAN RESOLUTION FUNCTION OF THREE AXIS NEUTRON     *
;*  SPECTROMETER.                   *
;*                      *
;* METHOD  :                    *
;*  THE CALCULATIONS ARE BASED ON REFS. 1-3 WITH THE MODIFICATIONS  *
;*  LISTED BELOW:                *
;*  --       THE INTENSITY PR.MONITOR IS PROPOTIONAL TO  THE   *
;*   CONVOLUTION OF THE RESOLUTION FUNCTION DELIVERED BY     *
;*   THIS ROUTINE FOLDED WITH THE SCATERING LAW S(Q,W),  *
;*      I.E. THE KF/KI FACTOR IS TAKEN INTO THE NORMALIZATION  *
;*   FACTOR OF THE RESOLUTION FUNCTION.         *
;*  THE MONITOR SENSITIVITY IS ASSUMED TO BE INVERSELY PROPOTIONAL  *
;*  TO THE INCIDENT WAVE VECTOR KI.            *
;*                      *
;* USE     :                    *
;*  RESELPS ASUMES THE PARAMETERS CHARACTERIZING THE EXPERIMENTAL   *
;*  CONFIGURATION TO BE AVAILABLE IN THE COMMON BLOCK,B   *
;*  THEY CAN BE READ IN BY THE PROCEDURE TASIN.         *
;*  FEIEF determines whether Ei or Ef are fixed :         *
;*   FEIEF=0     :    Ei fixed          *
;*   FEIEF=1     :    Ef fixed          *
;*                      *
;* RESULT  :                    *
;*  ARE DELIVERED IN R0, THE NORMALIZATION CONSTANT, AND RM(4,4) *
;*  THE CHARACTERISTIC MATRIX OF THE GAUSSAIN RESOLUTION FUNCTION.   *
;*  RM(4,4)'S INDEXES ARE EXPLAINED BELOW:          *
;*   1 : dq parallel to ki-kf=q          *
;*   2 : dq perpendicular to ki-kf=q, counter clockwise be   *
;*       positive.               *
;*   3 : de                 *
;*   4 : dq perpendicular to the scattering plane       *
;*  Wave vectors are measured in AA^-1 and energy in meV.      *

;*                      *
;* REFERENCE :                    *
;*      (1) : M.J.COOPER AND R.NATHANS. ACTA CRYST. 23, 357, (1967)    *
;*  (2) : S.A.WERNER AND R.PYNN. J.APPL.PHYS. 42,4736, (1971)    *
;*  (3) : N.J.CHESSER AND J.D.AXE. ACTA CRYST. A29,160,(1972)    *
;*                      *
;* INPUT     :                    *
;*   q   : momentum transfer in AA^-1   *
;*   e   : energy transfer in meV       *
;*                      *
;************************************************************************
;common resout,r0,rm
;common resinp,tm,ta,epm,ep,etam,etaa,alz,alm,ala,al3,vbetz,vbet1,vbet2,vbet3,$
;              etash,etasv,efix,feief,source
; horizontal focusing ? (if so, horifoc have to be 1)
;if n_elements(horifoc) eq 0 then horifoc=0


widget_control, event.top, get_uvalue = ptrAll
horifoc=(*ptrAll).horizontal
alz=(*ptrAll).collimatorh[0]
alm=(*ptrAll).collimatorh[1]
ala=(*ptrAll).collimatorh[2]
al3=(*ptrAll).collimatorh[3]
vbetz=(*ptrAll).collimatorv[0]
vbet1=(*ptrAll).collimatorv[1]
vbet2=(*ptrAll).collimatorv[2]
vbet3=(*ptrAll).collimatorv[3]
etash=(*ptrAll).mosaic[0]
etasv=(*ptrAll).mosaic[1]
etam=(*ptrAll).mosaic[2]
etaa=(*ptrAll).mosaic[3]
tm=2*!pi/(*ptrAll).dspacing[0]
ta=2*!pi/(*ptrAll).dspacing[1]
epm=(*ptrAll).orient[0]
ep=(*ptrAll).orient[1]
q=sqrt(((*ptrAll).qetrans[0])^2+((*ptrAll).qetrans[1])^2+((*ptrAll).qetrans[2])^2)
etran=(*ptrAll).qetrans[3]
efix=(*ptrAll).qetrans[4]

;     **********************************
;     *figure out fixed ei or fixed ef *
;     **********************************
      if (efix gt 0) then feief =   0 else feief = 1
      efix  =    abs(efix)

;     converts fwhm minutes to 1/e radians
      degtrad=!pi/180./sqrt(8.*alog(2.))/60.
      alz=alz*degtrad
      etam=etam*degtrad
      alm=alm*degtrad
      ala=ala*degtrad
      etaa=etaa*degtrad
      al3=al3*degtrad
      vbetz=vbetz*degtrad
      vbet1=vbet1*degtrad
      vbet2=vbet2*degtrad
      vbet3=vbet3*degtrad
      etash=degtrad*etash
      etasv=degtrad*etasv










;*****************
;* define arrays *
;*****************
rm  =    fltarr(4,4)
a   = fltarr(3,3)

;************************************************
;* definition of constants         *
;* E = 2.072 k^2 (E:meV, k:AA^-1)     *
;************************************************
f   = 1./2.072

;*************************
;* calculate ei,ef,ki,kf *
;* based on efix & feief *
;*************************
aom =   f*etran     ; energy transfer in units AA^-2

if feief eq 0 then begin

;   *************************
;   * fixed incident energy *
;   *************************
    ei =   efix
    ef =   efix-etran
endif   else  begin

;   **********************
;   * fixed final energy *
;   **********************
    ei =   efix+etran
    ef =   efix
endelse

ki  =    sqrt(f*ei)
kf  =    sqrt(f*ef)
alam    =  ki/kf

;**********************************************
;* do we have a 58Ni guide before the sample? *
;**********************************************
if alz le 0 then alzv   = 50./ki*(!dtor/sqrt(8.*alog(2.))/60.) $
        else alzv  =    alz

;****************************************************************
;* various trigonometric expressions from        *
;* scattering triangle  :    ki=kf-Q              *
;* (34) be = cos(2theta_s)  ;   al = sin(2theta_s)    *
;*  sints = sin(theta_s)               *
;* (29) psi=angle between ki and -Q ;   aom=ki^2-kf^2     *
;*  sa = sin(psi)       ; sb = cos(psi)     *
;*      b  = cos(2theta_s+psi)  ;    aa = sin(2theta_s+psi)     *
;****************************************************************
be  =    (ki^2+kf^2-q^2)/(2*ki*kf)  ; q^2=ki^2+kf^2-2ki*kf*be

if be^2 gt 1 then begin
;------------------------------------
;| Scattering triangle cannot close |
;------------------------------------
 r0 =   0
 return
endif

al  =    sqrt(1.0-be^2)
sints   = sqrt(0.5*(1.-be))
sb  =    (q*q+aom)/(2.0*q*ki)     ; kf^2=q^2+ki^2-2q*ki*sb
sa  =    sqrt(1.0-sb^2)
b   = (aom-q*q)/(2.0*q*kf)       ; ki^2=q^2+kf^2+2qkf*b
aa  =    sqrt(1.0-b^2)

;***************************************************************************
;* calculation of d (=alp), f (=bet) and h (=gam) and sa, sb (small a,and  *
;* small b) (29) and (39)                  *
;* 09/01/98 NEGATIVE SIGNS ARE ADDED TO sa,sb,aa,b BELOW SO THAT rm IS  *
;* GIVEN IN TERMS OF Q=ki-kf AND NOT Q_{C-N}=kf-ki AS BEFORE.   *
;***************************************************************************
alp =   [-b/al,aa/al,1./(al*2.0*kf)]
bet =   [-sb/al,sa/al,be/(al*2.0*kf)]
gam =   [0.,0.,-1./(2.*kf)]

;***************************************************************************
;* tangent of half scattering angle at monochrometer and analyzer, tom     *
;* and toa respectively.                    *
;* tm = 2pi/dm : reciprocal lattice spacing             *
;* ki = pi/(dm sin(theta_m)) --> sin(theta_m) = tm/(2ki)          *
;***************************************************************************
tom =   tm/(2.0*sqrt(ki^2-(tm/2.0)^2))*epm
toa =   (ep*ta/2.0)/sqrt(kf^2-(ta/2.0)^2)

;***************************************************************************
;* calculate the a parameters (42)               *
;***************************************************************************
a1  =    tom/(ki*etam)
a2  =    1.0/(ki*etam)
a3  =    1.0/(ki*alm)
a4  =    1.0/(kf*ala)
a5  =    toa/(kf*etaa)
a6  =    -(1.0-horifoc)/(kf*etaa)
a7  =    2.0*tom/(ki*alzv)
a8  =    1.0/(ki*alzv)
a9  =    2.0*toa/(al3*kf)*(1.-horifoc)
a10 =   -(1.0-horifoc)/(al3*kf)

;**************************************************************************
;* calculate the b parameters (44)              *
;**************************************************************************
b0  =    a1*a2+a7*a8
b1  =    a2*a2+a3*a3+a8*a8
b2  =    a4*a4+a6*a6+a10*a10
b3  =    a5^2+a9^2
b4  =    a5*a6+a9*a10
b5  =    a1^2+a7^2

;**************************************************************************
;* C,E,and A' (c,e,ap) (40a) (40b) and (45a)            *
;**************************************************************************
c   = -(alam-be)/al
e   = -(be*alam-1.0)/al
ap  =   2.0*b0*c+b1*c^2+b2*e^2+b3*alam^2+2.0*b4*alam*e+b5

;**************************************************************************
;* g parameters (here called d) (49)                *
;**************************************************************************
d0  =    b1-(b0+b1*c)^2/ap
d1  =    b2-(b2*e+b4*alam)^2/ap
d2  =    b3-(b3*alam+b4*e)^2/ap
d3  =    2.0*b4-2.0*(b2*e+b4*alam)*(b3*alam+b4*e)/ap
d4  =    -2.0*(b0+b1*c)*(b2*e+b4*alam)/ap
d5  =    -2.0*(b0+b1*c)*(b3*alam+b4*e)/ap


;***************************************************************************
;* final calculation of resolution matrix (51)             *
;***************************************************************************
rm(0:2,0:2) =   d0*alp#alp + d1*bet#bet + d2*gam#gam      $
         + 0.5*d3*(bet#gam + gam#bet)        $
         + 0.5*d4*(alp#bet + bet#alp)        $
         + 0.5*d5*(alp#gam + gam#alp)
rm(2,*) =   rm(2,*)*f
rm(*,2) =   rm(*,2)*f

;*****************************************************************************
;* calculate vertical resolution    (55), (57) ref.1,3          *
;* there are some typographical errors in the ref.1 that have been       *
;* fixed in ref.3                     *
;*****************************************************************************
snm =   tm/(2.0*ki)
sna =   ta/(2.0*kf)
a11 =   1/((2.0*snm*etam)^2+vbetz^2)/ki^2+1/(vbet1*ki)^2
a12 =   1/((2.0*sna*etaa)^2+vbet3^2)/kf^2+1/(vbet2*kf)^2
rm(3,3) =   a11*a12/(a11+a12)

;*****************************************************************************
;* calculation of resolution function prefactor and             *
;* modification for finite sample mosaic               *
;* see S.A. Werner and R.pynn J.Appl.phys. 42,4736,1971             *
;* section C. formulae (18), (19)                  *
;* error corrected 3/15/92 : previously multiplied instead of dividing by    *
;*            monc(ei)                 *
;* error corrected 6/26/98 : in r0, previously sin(theta_s) now sin(2theta_s)*
;*****************************************************************************
gamsh=(q*etash)^2
gamsv=(q*etasv)^2
denom=rm(1,1)*gamsh+1
r0  =             pa(ef)/             $
         (ki^3*kf*al)/sqrt(ap*(a11+a12))          $
    *   vbetz/             $
         sqrt(vbetz^2 + (2*etam*snm)^2)          $
    *   vbet3/             $
         sqrt(vbet3^2 + (2*etaa*sna)^2)          $
    *  sqrt(vbetz^2+vbet1^2+(2*etam*snm)^2)/         $
         (vbetz*vbet1)              $
    *  sqrt(alzv^2+alm^2+(2*etam)^2)/          $
               (alzv*alm*etam)             $
    *  abs(tom)/monc(ei, event)              $
    *       1./              $
       sqrt((1.+gamsv*rm(3,3))*(1.+gamsh*rm(1,1)))
if r0 le 0 then print,r0

for i=0,2 do for j=0,2 do a(i,j)=rm(i,j)-rm(i,1)*rm(j,1)*gamsh/denom
rm(0:2,0:2) =   a
rm(3,3) =   rm(3,3)-gamsv*rm(3,3)*rm(3,3)/(rm(3,3)*gamsv+1)

;print, rm
(*ptrAll).rm = rm
(*ptrAll).r0 = r0

end







function resolution_gui_attrib_assign

   attrib = $
      { psym:8 $     ; 0 - 8
      , selsym:0 $   ; 0 - 9
      , linestyle:0 $   ; 0 - 5
      , thick:1. $
      }

   return, attrib

end



pro resolution_gui_plotAll, event


widget_control, event, get_uvalue=ptrAll
wset, (*ptrAll).winpix
erase, (*ptrAll).erase_color
qellipse, (*ptrAll).rm, coordtitle=(*ptrAll).title
wset, (*ptrAll).wID1
;; copy from pixmap to visible window
device,copy=[0,0,!d.x_size,!d.y_size,0,0,(*ptrAll).winpix]

end




pro resolution_gui_event_open, event

   file = dialog_pickfile()
   if file eq '' then return
   ;     list of symbols
;     aki incident neutron wavevector  (aa-1)
;     tm monochromator tau (aa-1)
;     ta analyzer tau (aa-1)
;     epm +1. if sample scattering direction opposite of monochr. sc. di
;     ep +1. if analyzer scatt.dirct. same as monochr.sc.dirct.
;     if ep or epm is not +1, it is -1.
;     alz etc. horizontal divergences fwhm in minutes
;     vbet1 etc. vertical divergences fwhm in minutes
;     for ei fixed ein is positive.
;     for ef fixed ein is negative.  ef=abs(ein)

;print,'Note that input file has been modified to require source specification'
;print,'after the horizontal AND vertical sample mosaic:'
;print,'0 for cold source, 1 for thermal source'
openr,lun,file,/get_lun
      readf,lun,epm,ep,alz,alm,ala,al3,tm,ta
      readf,lun,etam,etaa,efix,etash,etasv,source
      readf,lun,vbetz,vbet1,vbet2,vbet3
free_lun,lun
widget_control, event, get_uvalue = ptrAll
(*ptrAll).collimatorh[0]=alz
(*ptrAll).collimatorh[1]=alm
(*ptrAll).collimatorh[2]=ala
(*ptrAll).collimatorh[3]=al3
(*ptrAll).collimatorv[0]=vbetz
(*ptrAll).collimatorv[1]=vbet1
(*ptrAll).collimatorv[2]=vbet2
(*ptrAll).collimatorv[3]=vbet3
(*ptrAll).mosaic[0]=etash
(*ptrAll).mosaic[1]=etasv
(*ptrAll).mosaic[2]=etam
(*ptrAll).mosaic[3]=etaa
(*ptrAll).dspacing[0]=2*!pi/tm
(*ptrAll).dspacing[1]=2*!pi/ta
(*ptrAll).orient[0]=epm
(*ptrAll).orient[1]=ep
(*ptrAll).qetrans[4]=efix

   wID = widget_info( event, find_by_uname = 'Colh_alz' )
   widget_control, wID, set_value = (*ptrAll).collimatorh[0]
   wID = widget_info( event, find_by_uname = 'Colh_alm' )
   widget_control, wID, set_value = (*ptrAll).collimatorh[1]
   wID = widget_info( event, find_by_uname = 'Colh_ala' )
   widget_control, wID, set_value = (*ptrAll).collimatorh[2]
   wID = widget_info( event, find_by_uname = 'Colh_al3' )
   widget_control, wID, set_value = (*ptrAll).collimatorh[3]

   wID = widget_info( event, find_by_uname = 'Colv_vbetz' )
   widget_control, wID, set_value = (*ptrAll).collimatorv[0]
   wID = widget_info( event, find_by_uname = 'Colv_vbet1' )
   widget_control, wID, set_value = (*ptrAll).collimatorv[1]
   wID = widget_info( event, find_by_uname = 'Colv_vbet2' )
   widget_control, wID, set_value = (*ptrAll).collimatorv[2]
   wID = widget_info( event, find_by_uname = 'Colv_vbet3' )
   widget_control, wID, set_value = (*ptrAll).collimatorv[3]

   wID = widget_info( event, find_by_uname = 'Mosaic_etash' )
   widget_control, wID, set_value = (*ptrAll).mosaic[0]
   wID = widget_info( event, find_by_uname = 'Mosaic_etasv' )
   widget_control, wID, set_value = (*ptrAll).mosaic[1]
   wID = widget_info( event, find_by_uname = 'Mosaic_etam' )
   widget_control, wID, set_value = (*ptrAll).mosaic[2]
   wID = widget_info( event, find_by_uname = 'Mosaic_etaa' )
   widget_control, wID, set_value = (*ptrAll).mosaic[3]

   wID = widget_info( event, find_by_uname = 'Dspacings_dm' )
   widget_control, wID, set_value = (*ptrAll).dspacing[0]
   wID = widget_info( event, find_by_uname = 'Dspacings_da' )
   widget_control, wID, set_value = (*ptrAll).dspacing[1]

   wID = widget_info( event, find_by_uname = 'Instrument_orientation_epm' )
   widget_control, wID, set_value = (*ptrAll).orient[0]
   wID = widget_info( event, find_by_uname = 'Instrument_orientation_ep' )
   widget_control, wID, set_value = (*ptrAll).orient[1]

   wID = widget_info( event, find_by_uname = 'Qetransfer_ein' )
   widget_control, wID, set_value = (*ptrAll).qetrans[4]
   wID = widget_info( event, find_by_uname = 'Source' )
   widget_control, wID, set_value = (*ptrAll).source







   resolution_gui_plotAll, event

end


pro resolution_gui_motion, event
widget_control, event.top, get_uvalue=ptrAll
;if not ptr_valid(state.dataptr) then return
coords=convert_coord(event.x, event.y,/device,/to_data)
;get two element coordinate array (x,y)
x=coords[0]
y=coords[1]
;The following line would print the coordinate of the cursor as the title
;of the resolution ellipse graph
;(*ptrAll).title=strcompress('('+string(x)+','+string(y)+')',/remove_all)

end






pro resolution_resize, event
widget_control, event.top, get_uvalue=ptrAll
;; What is the current size of the TLB?
;; forces the minimum size to 100 x 100 ">" chooses max of the two values
newxsize=event.x > 100  & newysize=event.y > 100
;print, newxsize, newysize
;;Destroy the old pixmap
wdelete, (*ptrAll).winpix
;; Create a new pixmap of the newsize
window,/free,/pixmap, xsize=newxsize, ysize=newysize
(*ptrAll).winpix=!d.window ; update the state structure
;;get the widget id of the draw widget
win_id=widget_info(event.top, find_by_uname='WIN')
widget_control, win_id, draw_xsize=newxsize, draw_ysize=newysize
;; Send the updated state structure back to the TLB UVALUE
(*ptrAll).draw_xsize=newxsize
(*ptrAll).draw_ysize=newysize



;print, 'resize'
end


pro resolution_gui_event, event

;   widget_control, event.id, get_uvalue = uvalue
;
;   if n_elements( uvalue ) ne 0 then message, uvalue, /continue else return

   uname=widget_info(event.id,/uname)
   if n_elements(uname) eq 0 then return
;   if n_elements( uname ) ne 0 then message, uname, /continue else return
   widget_control, event.top, get_uvalue = ptrAll
   case uname of

      ; MENU BAR
      'Open':begin
      resolution_gui_event_open, event.top
      ;print, (*ptrAll).collimatorv[3]
      end
      'Ps': begin
;         set_plot, 'ps'
;         device, /portrait, /color, file = 'idl.ps' $
;            , xsize = 8.5, ysize = 11, /inches, xoffset = 0, yoffset = 0
;         resolution_gui_plotAll, event.top
;         device, /close
;         set_plot, 'win'
     reselps, event
; the below should have made a nice postscript file.  Perhaps a future programmer
; can figure out why the output isn't pretty

;     keywords = PSConfig(Cancel=cancelled)
;   IF NOT cancelled THEN BEGIN
;      thisDevice = !D.Name
;      Set_Plot, "PS"
;      Device, _Extra=keywords
;      qellipseprinter, (*ptrAll).rm, coordtitle=(*ptrAll).title
;      Device, /Close_File
;      Set_Plot, thisDevice
;   ENDIF

      end
      'Print': begin

         if not dialog_printersetup() then return
         set_plot, 'printer'
         device, /portrait $
            , xsize = 8.5, ysize = 11, /inches, xoffset = 0, yoffset = 0
         qellipseprinter, (*ptrAll).rm, coordtitle=(*ptrAll).title
         device, /close_document
         set_plot, 'win'

      end
      'Quit': widget_control, event.top, /destroy

      'Title': begin
         widget_control, event.id, get_value = value
         (*ptrAll).title = value
         resolution_gui_plotAll, event.top
      end
      'Xtitle': begin
         widget_control, event.id, get_value = value
         (*ptrAll).xtitle = value
         resolution_gui_plotAll, event.top
      end
      'Ytitle': begin
         widget_control, event.id, get_value = value
         (*ptrAll).ytitle = value
         resolution_gui_plotAll, event.top
      end

      'Colh_alz': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorh[0] = value
         resolution_gui_plotAll, event.top
      end
      'Colh_alm': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorh[1] = value
         resolution_gui_plotAll, event.top
      end
      'Colh_ala': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorh[2] = value
         resolution_gui_plotAll, event.top
      end
      'Colh_al3': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorh[3] = value
         resolution_gui_plotAll, event.top
      end

      'Colv_vbetz': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorv[0] = value
         resolution_gui_plotAll, event.top
      end
      'Colv_vbet1': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorv[1] = value
         resolution_gui_plotAll, event.top
      end
      'Colv_vbet2': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorv[2] = value
         resolution_gui_plotAll, event.top
      end
      'Colv_vbet3': begin
         widget_control, event.id, get_value = value
         (*ptrAll).collimatorv[3] = value
         resolution_gui_plotAll, event.top
      end

      'Mosaic_etash': begin
         widget_control, event.id, get_value = value
         (*ptrAll).mosaic[0] = value
         resolution_gui_plotAll, event.top
      end
      'Mosaic_etasv': begin
         widget_control, event.id, get_value = value
         (*ptrAll).mosaic[1] = value
         resolution_gui_plotAll, event.top
      end
      'Mosaic_etam': begin
         widget_control, event.id, get_value = value
         (*ptrAll).mosaic[2] = value
         resolution_gui_plotAll, event.top
      end
      'Mosaic_etaa': begin
         widget_control, event.id, get_value = value
         (*ptrAll).mosaic[3] = value
         resolution_gui_plotAll, event.top
      end

      'Dspacings_dm': begin
         widget_control, event.id, get_value = value
         (*ptrAll).dspacing[0] = value
         resolution_gui_plotAll, event.top
      end
      'Dspacings_da': begin
         widget_control, event.id, get_value = value
         (*ptrAll).dspacing[1] = value
         resolution_gui_plotAll, event.top
      end
      'Instrument_orientation_epm': begin
         widget_control, event.id, get_value = value
         (*ptrAll).orient[0] = value
         resolution_gui_plotAll, event.top
      end
      'Instrument_orientation_ep': begin
         widget_control, event.id, get_value = value
         (*ptrAll).orient[1] = value
         resolution_gui_plotAll, event.top
      end
      'Qetransfer_qh': begin
         widget_control, event.id, get_value = value
         (*ptrAll).qetrans[0] = value
         resolution_gui_plotAll, event.top
      end
      'Qetransfer_qk': begin
         widget_control, event.id, get_value = value
         (*ptrAll).qetrans[1] = value
         resolution_gui_plotAll, event.top
      end
      'Qetransfer_ql': begin
         widget_control, event.id, get_value = value
         (*ptrAll).qetrans[2] = value
         resolution_gui_plotAll, event.top
      end


      'Qetransfer_etran': begin
         widget_control, event.id, get_value = value
         (*ptrAll).qetrans[3] = value
         resolution_gui_plotAll, event.top
      end
      'Qetransfer_ein': begin
         widget_control, event.id, get_value = value
         (*ptrAll).qetrans[4] = value
         resolution_gui_plotAll, event.top
      end
      'Source': begin
         widget_control, event.id, get_value = value
         (*ptrAll).source = value
         resolution_gui_plotAll, event.top
      end
      'Horizontal_focusing': begin
         widget_control, event.id, get_value = value
         (*ptrAll).horizontal = value
         resolution_gui_plotAll, event.top
      end
      'Calculate': begin
       reselps, event
       eigenvalues = float(LA_EIGENPROBLEM((*ptrAll).rm, EIGENVECTORS = eigenvectors))
       eigenvectors=float(eigenvectors)
       eiginv=invert(eigenvectors)
        ;diagonal=float(eiginv#f#eigenvectors)
       d=sqrt(2*alog(2.0))
       a=d/sqrt(eigenvalues[0])*100.0
       b=d/sqrt(eigenvalues[1])*100.0
       c=d/sqrt(eigenvalues[2])*10.0
       transmatrix=eigenvectors; note that this is a 4x4 matrix since it comes from the R-matrix
       ; delete the last row
       transmatrix=transmatrix[*,[0,1,2]]
       ; delete the last column
       transmatrix=transmatrix[[0,1,2],*]


       ;print, 'Eigen', transmatrix
       ;print, 'RMatrix', (*ptrAll).rm
;transmatrix=[[1.0,0.0,0.0],[0.0,cos(!dtor*30.0),sin(!dtor*30.0)], [0.0,-sin(!dtor*30.0), cos(!dtor*30.0)]]


       plot_ellipse, a, b, c, transmatrix, (*ptrAll).rm, (*ptrAll).qetrans[0],   $
       (*ptrAll).qetrans[1], exact=[1b,1b,1b],group_leader = event.top
; plot cross sections
    resolution_gui_plotAll, event.top
       end
    'TLB': begin
    resolution_resize, event
    resolution_gui_plotall, event.top
    end
    'WIN' :begin
    reselps, event
    resolution_gui_motion, event
    resolution_gui_plotall, event.top
    end
      else: return
   endcase

end









pro Resolution_Gui_cleanup, tlb

   message, 'Cleanup begins...', /continue
   widget_control, tlb, get_uvalue = ptrAll
   wdelete, (*ptrall).winpix
   heap_free, ptrAll
   help,/heap,/brief
;   heap_gc, /verbose

end



PRO ResolutionGUI,group_leader = group_leader



   ;draw_xsize = 480
   ;draw_ysize = 621
   draw_xsize = 680
   draw_ysize = 621
   color_table = 0
   plot_color = 0
   erase_color = color24([ 255, 255, 255 ])
   position = [ .2, .3, .8, .7 ]
   xrange = [0d, 0d]
   yrange = [0d, 0d]
   collimatorh=[60,60, 60, 60]
   collimatorv=[120,120, 120, 120]
   mosaic=[10,10, 30, 30]
   dspacing=[3.35, 3.35]
   qetrans = [1.5,0.0,0.0, 6.0, -5.0]
   orient = [1, 1]
   source = 0
   horizontal = 0
   charsize = 1.
   symscale = 1.
   rm  =    fltarr(4,4)
   r0 = 10.0
   feief=0
   tab_xsize = 190
   tab_ysize = draw_ysize

   title = ''
   xtitle = 'X Title'
   ytitle = 'Y Title'

   selectionArr = [ 'Select the data set.' ]

   device, decomposed = 1




   ; TOP-LEVEL-BASE
   screen_size = get_screen_size()
   screen_left = ( screen_size[0] - draw_xsize - tab_xsize ) / 2
   screen_top = ( screen_size[1] - draw_ysize ) / 2
   if n_elements(group_leader) eq 0 then group_leader = 0L
   tlb = widget_base( mbar = mbar, /row, xoffset = screen_left, yoffset = screen_top, $
    /tlb_size_events, uname='TLB',group_leader = group_leader)

   ; MENU BAR
   mFile = widget_button( mbar, value = 'File', /menu )
   mOpen = widget_button( mFile, value = 'Open', uname='Open', uvalue = 'open' )
;   mConvert = widget_button( mFile, /separator, /menu,uname='Convert', value = 'Convert' )
;   mPS = widget_button( mConvert, value = 'To PostScript', uname='Ps',uvalue = 'ps' )
   mPrint = widget_button( mFile, /separator, value = 'Print', uname='Print',uvalue = 'print' )
   mQuit = widget_button( mFile, /separator, value = 'Quit',uname='Quit', uvalue = 'quit' )

    ;tld=widget_base (tlb, /row)
    wDraw1 = widget_draw( tlb, xsize = draw_xsize, ysize = draw_ysize, $
    graphics_level=1, renderer=1,/motion_events, uname='WIN' )
   ;wTab = widget_tab( tlb, xsize = tab_xsize, location = 0 )
   wTab=widget_tab(tlb, location=0)
   ;tInstrumentParams = widget_base( wTab, title = 'Instrument Parameters', /column )
   tExperiment = widget_base( wTab, title = 'Experiment', /column )




   tHorizontalParams = widget_base(wTab, title = 'Horizontal Collimators', /column )
   wColh = widget_base( tHorizontalParams, frame = 3, /base_align_right, /column, xpad = 3 )
   ;void = widget_label( wColh, value = 'Horizontal Collimators', /align_center )
   bColh1R = cw_field( wColh, /all_events, title = 'ALZ' $
      , value = collimatorh[0], /integer, uname='Colh_alz', uvalue = 'colh_alz', xsize=4)
   bColh2R = cw_field( wColh, /all_events, title = 'ALM' $
      , value = collimatorh[1], /integer, uname='Colh_alm', uvalue = 'colh_alm', xsize=4 )
   bColh3R = cw_field( wColh, /all_events, title = 'ALA' $
      , value = collimatorh[2], /integer, uname='Colh_ala', uvalue = 'colh_ala', xsize=4 )
   bColh4R = cw_field( wColh, /all_events, title = 'AL3' $
      , value = collimatorh[3], /integer, uname='Colh_al3', uvalue = 'colh_al3', xsize=4 )


   tVerticalParams = widget_base(wTab, title = 'Vertical Collimators', /column )
   wColv = widget_base( tVerticalParams, frame = 3, /base_align_right, /column, xpad = 3 )
   ;void = widget_label( wColv, value = 'Collimation Vertical', /align_center )

   bColv1R = cw_field( wColv, /all_events, title = 'VBETZ', xsize=4 $
      , value = collimatorv[0], /integer, uname = 'Colv_vbetz', uvalue = 'colv_vbetz' )

   bColv2R = cw_field( wColv, /all_events, title = 'VBET1', xsize=4 $
      , value = collimatorv[1], /integer, uname = 'Colv_vbet1', uvalue = 'colv_vbet1' )

   bColv3R = cw_field( wColv, /all_events, title = 'VBET2', xsize=4 $
      , value = collimatorv[2], /integer, uname = 'Colv_vbet2', uvalue = 'colv_vbet2' )

   bColv4R = cw_field( wColv, /all_events, title = 'VBET3', xsize=4 $
      , value = collimatorv[3], /integer, uname = 'Colv_vbet3', uvalue = 'colv_vbet3' )


   tMosaic = widget_base(wTab, title = 'Mosaics', /column )
   wMosaic = widget_base(tMosaic, frame = 3, /base_align_right, /column, xpad = 3 )
   ;void = widget_label( wMosaic, value = 'Mosaics', /align_center )

   bEtash1R = cw_field( wMosaic, /all_events, title = 'ETASH', xsize=4 $
      , value = mosaic[0], /integer, uname = 'Mosaic_etash', uvalue = 'mosaic_etash' )

   bEtasv2R = cw_field( wMosaic, /all_events, title = 'ETASV', xsize=4 $
      , value = mosaic[1], /integer, uname = 'Mosaic_etasv', uvalue = 'mosaic_etasv' )

   bEtam3R = cw_field( wMosaic, /all_events, title = 'ETAM', xsize=4 $
      , value = mosaic[2], /integer, uname = 'Mosaic_etam', uvalue = 'mosaic_etam' )

   bEtaa4R = cw_field( wMosaic, /all_events, title = 'ETAA' , xsize=4$
      , value = mosaic[3], /integer, uname = 'Mosaic_etaa', uvalue = 'mosaic_etaa' )



   tDspacings = widget_base(wTab, title = 'Monochormator D spacings', /column )
   wDspacings = widget_base( tDspacings, frame = 3, /base_align_right, /column, xpad = 3 )
   ;void = widget_label( wDspacings, value = 'Dspacings', /align_center )

   bDm1R = cw_field( wDspacings, /all_events, title = 'dm', xsize=4 $
      , value = dspacing[0], /floating, uname = 'Dspacings_dm', uvalue = 'dspacings_dm' )

   bDa2R = cw_field( wDspacings, /all_events, title = 'da', xsize=4 $
      , value = dspacing[1], /floating, uname = 'Dspacings_da', uvalue = 'dspacings_da' )


   tInstrument_orientation = widget_base(wTab, title = 'Instrument Orientation', /column )
   wInstrument_orientation = widget_base( tInstrument_orientation, frame = 3, /base_align_right, /column, xpad = 3 )
   ;void = widget_label( wInstrument_orientation, value = 'Instrument_orientation', /align_center )

   bEpm1R = cw_field( wInstrument_orientation, /all_events, title = 'epm', xsize=2 $
      , value = orient[0], /integer, uname = 'Instrument_orientation_epm', uvalue = 'Instrument_orientation_epm' )

   bEp2R = cw_field( wInstrument_orientation, /all_events, title = 'ep', xsize=2 $
      , value = orient[1], /integer, uname = 'Instrument_orientation_ep', uvalue = 'Instrument_orientation_ep' )


   wQetransfer = widget_base( tExperiment, frame = 3, /base_align_right, /column, xpad = 3 )
   void = widget_label( wQetransfer, value = 'Qetransfer', /align_center )

   bQh1R = cw_field( wQetransfer, /all_events, title = 'qh' $
      , value = qetrans[0], /floating, uname = 'Qetransfer_qh', uvalue = 'qetransfer_qh' )
   bQk2R = cw_field( wQetransfer, /all_events, title = 'qk' $
      , value = qetrans[1], /floating, uname = 'Qetransfer_qk', uvalue = 'qetransfer_qk' )
   bQl3R = cw_field( wQetransfer, /all_events, title = 'ql' $
      , value = qetrans[2], /floating, uname = 'Qetransfer_ql', uvalue = 'qetransfer_ql' )

   bEtran4R = cw_field( wQetransfer, /all_events, title = 'etran' $
      , value = qetrans[3], /floating, uname = 'Qetransfer_etran', uvalue = 'qetransfer_etran' )
   bEin5R = cw_field( wQetransfer, /all_events, title = 'ein' $
      , value = qetrans[4], /floating, uname = 'Qetransfer_ein', uvalue = 'qetransfer_ein' )


   wSource_value = widget_base( tExperiment, frame = 3, /base_align_right, /column, xpad = 3 )
   void = widget_label( wSource_value, value = 'Source', /align_center )

   bSource1R = cw_field( wSource_value, /all_events, title = 'Source', xsize=2 $
      , value = source, /integer, uname = 'Source', uvalue = 'source' )

   wHorizontal_value = widget_base( tExperiment, frame = 3, /base_align_right, /column, xpad = 3 )
   void = widget_label( wHorizontal_value, value = 'Horizontal focusing', /align_center )

   bSource1R = cw_field( wHorizontal_value, /all_events, title = 'Horizontal focus', xsize=2 $
      , value = horizontal, /integer, uname = 'Horizontal_focusing', uvalue = 'horizontal_focusing' )

    buttbase=widget_base(tlb, row=1)
    buttsize =70
    buttcalculate1=widget_button(buttbase, value='Calculate RM' $
    ,uname='Calculate',uvalue='calculate', xsize=buttsize)


widget_control, tlb, /realize

 widget_control, wDraw1, get_value = wID1
 window, /free,/pixmap,xsize=draw_xsize, ysize=draw_ysize
 winpix=!d.window
 wset, wID1
 erase, erase_color
All = $
      ; Draw Widget Related.
      { wID1:wID1 $
      , winpix:winpix $
      , draw_xsize:draw_xsize $
      , draw_ysize:draw_ysize $
      , color_table:color_table $
      , plot_color:plot_color $
      , title:title $
      , erase_color:erase_color $
      , xrange:xrange $
      , yrange:yrange $
      , position:position $
      , collimatorh:collimatorh $
      , collimatorv:collimatorv $
      , mosaic:mosaic $
      , dspacing:dspacing $
      , qetrans:qetrans $
      , orient:orient $
      , rm:rm $
      , r0:r0 $
      , source:source $
      , horizontal:horizontal $
      , feief:feief $
      , charsize:charsize $
      , symscale:symscale $

      ; Tab Widget Related.
      , tab_xsize:tab_xsize $

      ; Data Tab Related.
      , selection:ptr_new(selectionArr, /no_copy) $


      ; Data Related.
      , data:ptr_new(/allocate_heap) $
      , attrib:ptr_new(/allocate_heap) $
      , n_data:0L $

      }

   ptrAll = ptr_new( All, /no_copy )
widget_control, tlb, set_uvalue = ptrAll


xmanager, 'resolution_gui', No_Block=1, tlb $
      , event_handler = 'resolution_gui_event' $
      , cleanup = 'resolution_gui_cleanup'
end