; $Id: rescal2_calculator__define.pro,v 1.5 2007/02/17 00:10:33 kneller Exp $
;###############################################################################
;
;
;CLASS NAME:
;    Rescal2_calculator
;
;PURPOSE:
;
;CATEGORY:
;
;
;SUPERCLASSES:
;
;
;METHODS:
;    twoslitfwhm
;    fwhmellips
;    Rescal2_calculator::proj_scatplane
;    tasin
;    angle
;    Rescal2_calculator::monc
;    pa
;    phonint
;    oneddisp
;    eintqres
;    Rescal2_calculator::reselps
;    Rescal2_calculator::calculate
;    Rescal2_calculator::getproperty
;    Rescal2_calculator::setproperty
;    Rescal2_calculator::cleanup
;    Rescal2_calculator::init
;    Rescal2_calculator__define
;
;
; AUTHOR:
; Larry Kneller
; NIST Center for Neutron Research
; 100 Bureau Drive, Gaithersburg, MD 20899
; United States
; kneller@nist.gov  301-975-8839
; Thu May 18 15:20:57 2006
;
; 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.
;
;###############################################################################
;###############################################################################
; 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.
;
;###############################################################################
;###############################################################################
;
;NAME:
;        twoslitfwhm
;
;PURPOSE:
;           ;12/11/06   THIS SEEMS TO BE UNUSED.
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;RETURN VALUE:
;
;###############################################################################
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;twoslitfwhm


;###############################################################################
;
;NAME:
;        fwhmellips
;
;PURPOSE:
;           ;12/11/06   THIS SEEMS TO PRODUCE THE x,y OF THE RESOLUTION ELLIPSE GIVEN M.
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
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;fwhmellips



;###############################################################################
;
;NAME:
;        Rescal2_calculator::proj_scatplane
;
;PURPOSE:
;           ;12/11/06   DOES THIS ROTATE rm INTO THE SCATTERING PLANE????
;                       NO.  IT GIVES THE PROJECTION OF THE RES ELLIPSE ONTO THE SCATTERING PLANE.
;
;                       THIS DOES NOT APPEAR TO BE USED.
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::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

    rm = self.rm

    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;proj_scatplane


;###############################################################################
;
;NAME:
;        tasin
;
;PURPOSE:
;           ;12/11/06
;           ;DOES THIS READ AN ICP TAS FILE AND GET THE PARAMETERS?
;           ;I DON'T THINK SO.  PERHAPS IT READS IN A SPECIALIZED FILE
;           ;FOR INPUTTING RESOLUTION FUNCTION PARAMETERS.
;           ;BUT THEN IT DOES NOT SEEM TO DO ANYTHING WITH THEM.
;
;           ;AH!  THIS USED TO READ THE DATA AND PLACE IT INTO A COMMON BLOCK.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
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;tasin


;###############################################################################
;
;NAME:
;        angle
;
;PURPOSE:
;           ;12/11/06
;           ;GET THE ANGLE IN REGULAR RADIANS BETWEEN A RISE AND A RUN.
;           ;THIS DOES NOT APPEAR TO BE USED.
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;RETURN VALUE:
;
;###############################################################################
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;angle


;###############################################################################
;
;NAME:
;        Rescal2_calculator::monc
;
;PURPOSE:
;           ;12/11/06
;           ;DOES A MONITOR CORRECTION SPECIFIC TO BT4.
;           ;THIS IS CURRENTLY USED IN reselps.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;RETURN VALUE:
;
;###############################################################################
function Rescal2_calculator::monc,emon,feief,source
;****************************************************************
;* 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  = self.feief
;    source = self.source

    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;monc


;###############################################################################
;
;NAME:
;        pa
;
;PURPOSE:
;           ;12/11/06
;           ;THIS IS CURRENTLY USED IN reselps.  IT JUST RETURNS 1.
;           ;IT SEEMS THAT IT HAD AN INTENDED PURPOSE THAT IT WOULD
;           ;EVENTUALLY MEET, BUT IT WAS NEVER UPDATED.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;RETURN VALUE:
;
;###############################################################################
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;pa


;APPEARS TO BE UNUSED
;###############################################################################
;
;NAME:
;        phonint
;
;PURPOSE:
;           ;12/11/06
;           ;THIS PROCEDURE IS NOT USED.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
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;phonint



;APPEARS TO BE UNUSED
;###############################################################################
;
;NAME:
;        oneddisp
;
;PURPOSE:
;           ;12/11/06
;           ;THIS SEEMS TO CALCULATE THE 1D WIDTH GIVEN rm.
;           ;CURRENTLY IT DOES NOT APPEAR TO BE USED.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
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;oneddisp



;APPEARS TO BE UNUSED
;###############################################################################
;
;NAME:
;        eintqres
;
;PURPOSE:
;           ;12/11/06
;           ;THIS IS NOT USED.
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
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;eintqres


;###############################################################################
;
;NAME:
;        Rescal2_calculator::reselps
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::reselps,closed
;************************************************************************
;* 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

;print,'###############################################'
;print,'IN reselps:'

    closed = 1  ;ASSUME SCATTERING TRIANGLE CLOSES  ---- CHECK IS BELOW

    source = self.source
    ;widget_control, event.top, get_uvalue = ptrAll
    horifoc=self.horizontal;(*ptrAll).horizontal
    ;print,'horifoc=',horifoc
    alz=self.collimatorh[0];(*ptrAll).collimatorh[0]
    ;print,'alz=',alz
    alm=self.collimatorh[1];(*ptrAll).collimatorh[1]
    ;print,'alm=',alm
    ala=self.collimatorh[2];(*ptrAll).collimatorh[2]
    ;print,'ala=',ala
    al3=self.collimatorh[3];(*ptrAll).collimatorh[3]
    ;print,'al3=',al3
    vbetz=self.collimatorv[0];(*ptrAll).collimatorv[0]
    ;print,'vbetz=',vbetz
    vbet1=self.collimatorv[1];(*ptrAll).collimatorv[1]
    ;print,'vbet1=',vbet1
    vbet2=self.collimatorv[2];(*ptrAll).collimatorv[2]
    ;print,'vbet2=',vbet2
    vbet3=self.collimatorv[3];(*ptrAll).collimatorv[3]
    ;print,'vbet3=',vbet3
    etash=self.mosaic[0];(*ptrAll).mosaic[0]
    ;print,'etash=',etash
    etasv=self.mosaic[1];(*ptrAll).mosaic[1]
    ;print,'etasv=',etasv
    etam=self.mosaic[2];(*ptrAll).mosaic[2]
    ;print,'etaa=',etam
    etaa=self.mosaic[3];(*ptrAll).mosaic[3]
    ;print,'etaa=',etaa
    tm=2.0*!pi/self.dspacing[0];(*ptrAll).dspacing[0]
    ;print,'dm???=',(*ptrAll).dspacing[0]
    ta=2.0*!pi/self.dspacing[1]
    ;print,'da???=',(*ptrAll).dspacing[1]
    epm=self.instrument_orient[0];(*ptrAll).orient[0]
    ;print,'epm=',epm
    ep =self.instrument_orient[1]
    ;print,'ep=',ep


    ;widget_control,(*ptrall).a,get_value=a


    q=sqrt((self.qetrans[0])^2+(self.qetrans[1])^2+(self.qetrans[2])^2)
    ;q in inverse Angstroms
;;;;;;    print,'orig q='
;;;;;;    print,q
    ;print, q


;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################

    ;THIS IS CHANGED SO R IS IN X,Y,w,Z COORDINATES
;NOTE:  THIS MEANS THAT q IS IN THE MEASUREMENT PLANE IN THE COORDINATES
;       OF THE CRYSTAL.  SO THE ONLY DIFFERENCE BETWEEN THIS AND THE
;       AXES OF THE CRYSTAL IS AN OVERALL ROTATION OF THE CRYSTAL
;       INDEPENDENT OF THE SPECTROMETER.


    ;q is in x,y,z
    ;q=self.mylattice->modvec([self.qetrans[0],self.qetrans[1],self.qetrans[2]])
    ;q now in Recip Lattice Units

;;;;;;    print,'new q='
;;;;;;    print,q
    ;print, 'new q='
    ;print, q

;IS THIS DONE CORRECTLY??????????????????????????????????????????????????????????
;IF NOT, THEN THE REST OF THE WELL-TESTED CODE BELOW IS INCORRECT AS WELL!!!!!!!!!!!


;################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################
;#################################################################################

    ;print,'qh???=',(*ptrAll).qetrans[0]
    ;print,'qk???=',(*ptrAll).qetrans[1]
    ;print,'ql???=',(*ptrAll).qetrans[2]

    etran=self.qetrans[3]
    ;print,'etran=',etran
    efix=self.qetrans[4]
    ;print,'efix???=',efix

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

          self.feief = feief

    ;     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

;    print,'ki=',ki
;    print,'kf=',kf
;    print,'Q= ',Q'

;print,'THIS IS A TEST'

    if be^2 gt 1 then begin
        ;------------------------------------
        ;| Scattering triangle cannot close |
        ;------------------------------------
        ;void = dialog_message('rescal2_calculator::reselps  Scattering triangle did not close!')

;;;        print, 'Rescal2_calculator::reselps     DID Not Close'
;;;print,'THIS IS A TEST'

         closed = 0
         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)/self->monc(ei,feief,source) $    ;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,'###########'
;    print,'In RESELPS:'
;    print, rm


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

end;reselps
;###############################################################################
;
;NAME:
;        Rescal2_calculator::calculate
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::calculate,rm,r0,$
                                  horizontal=horizontal,$
                                  collimatorh=collimatorh,$
                                  collimatorv=collimatorv,$
                                  mosaic=mosaic,$
                                  dspacing=dspacing,$
                                  instrument_orient=instrument_orient,$
                                  qetrans=qetrans,$
                                  lattice=lattice,$
                                  closed=closed

    if n_elements(horizontal) eq 4 then self.horizontal = horizontal
    if n_elements(collimatorh) eq 4 then self.collimatorh = collimatorh
    if n_elements(collimatorv) eq 4 then self.collimatorv = collimatorv
    if n_elements(mosaic) eq 4 then self.mosaic = mosaic
    if n_elements(dspacing) eq 2 then self.dspacing = dspacing
    if n_elements(instrument_orient) eq 2 then self.instrument_orient=instrument_orient
    if n_elements(qetrans) eq 5 then self.qetrans = qetrans
    if n_elements(lattice) eq 1 then self.mylattice = lattice


;    print,'RESCAL2_CALCULATOR::CALCULATE'
;    print,self.horizontal
;    print,self.collimatorh
;    print,self.collimatorv
;    print,self.mosaic
;    print,self.dspacing
;    print,self.instrument_orient
;    print,self.qetrans
;    help,self.mylattice

    self->reselps,closed

    if closed eq 0 then begin
        self.r0 = 1.0
        self.rm = double([[1,0,0,0],$
                          [0,1,0,0],$
                          [0,0,1,0],$
                          [0,0,0,1]])
        void = dialog_message('Scattering triangle did not close!')

    endif;closed

    rm = self.rm
    r0 = self.r0

end;calculate

;###############################################################################
;
;NAME:
;        Rescal2_calculator::getproperty
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::getproperty,_ref_extra=extra

    if arg_present(mylattice) ne 0 then mylattice = self.mylattice

end;getProperty
;###############################################################################
;
;NAME:
;        Rescal2_calculator::setproperty
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::setproperty,_ref_extra=extra


    if n_elements(lattice) gt 0 then begin
        obj_destroy,self.mylattice
        self.mylattice = lattice
    endif

end;getProperty


;###############################################################################
;
;NAME:
;        Rescal2_calculator::cleanup
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator::cleanup
    ;print,'Rescal2_calculator::cleanup'
    if obj_valid(self.mylattice) gt 0 then obj_destroy,self.mylattice
end;cleanup
;###############################################################################
;
;NAME:
;        Rescal2_calculator::init
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;RETURN VALUE:
;
;###############################################################################
function Rescal2_calculator::init,$
                                lattice,$;a, b, c, alpha, beta, gamma,$           ;LATTICE PARAMETERS
                                horizontal,collimatorh,collimatorv,$    ;RESOLUTION PARAMETERS
                                mosaic,dspacing,instrument_orient,$
                                qetrans,feief,source,report=report


;                           help,lattice,$
;                                horizontal,collimatorh,collimatorv,$
;                                mosaic,dspacing,instrument_orient,$
;                                qetrans,feief,source

    if n_elements(report) eq 0 then report = 0

    if n_elements(horizontal) eq 1 then begin
        self.horizontal = horizontal
    endif else print,'Rescal2_calculator problem: horizontal '

    if n_elements(collimatorh) eq 4 then begin
        self.collimatorh = collimatorh
    endif else print,'Rescal2_calculator problem: collimatorh '

    if n_elements(collimatorv) eq 4 then begin
        self.collimatorv = collimatorv
    endif else print,'Rescal2_calculator problem: collimatorv '

    if n_elements(mosaic) eq 4 then begin
        self.mosaic = mosaic
    endif else print,'Rescal2_calculator problem: mosaic '

    if n_elements(dspacing) eq 2 then begin
        self.dspacing = dspacing
    endif else print,'Rescal2_calculator problem: dspacing '

    if n_elements(instrument_orient) eq 2 then begin
        self.instrument_orient = instrument_orient
    endif else print,'Rescal2_calculator problem: instrument_orient '

    if n_elements(qetrans) eq 5 then self.qetrans = qetrans

    if n_elements(feief) eq 1 then self.feief = feief else feief = 0

    if n_elements(source) eq 1 then self.source = source else source = 0


    self.mylattice = lattice;obj_new('lattice', a, b, c, alpha, beta, gamma)


    ;rm = [ x   ;X ALONG SCATTERING PLANE ORIENTATION VECTOR 1
    ;       y   ;Y ALONG SCATTERING PLANE ORIENTATION VECTOR 2
    ;       E   ;ENERGY TRANSFER
    ;       z]  ;Z ALONG SPECTROMETER z DIRECTION.



    ;ONCE I HAVE rm IN THIS FORMAT, FOR EACH POINT IN MY SCAN, I SIMPLY
    ;CALCULATE

    ;1) FIRST I CALCULATE rm FOR EACH POINT IN SCAN
    ;2) THEN I SELECT A DELTA FUNCTION FEATURE AT SOME POINT IN THE SCAN (i.e. SELECT Xo,Yo,Eo)
    ;3) FINALLY I CALCULATE exp((1/2)*[Xo,Yo,E,0]#(rm#[Xo,Yo,E,0])) FOR EACH RM IN THE SCAN.


    ;HOW ABOUT HALF WIDTH?
    ;1/2 = Ro exp((-1/2)[X,Y,E,0]#(rm#[X,Y,E,0])))
    ;
    ;
    ;WHAT IS UNITS CONVERSION??? TAKEN CARE OF IN RES MATRIX CALCULATION (see reselps)
    ;
    ;DO I CONVERT E TO Q OR VICE-VERSA?



    if report ne 0 then begin

        print,'#################################################################'
        print,'Rescal2_calculator::init ---- REPORTING INITIAL VARIABLE VALUES:'
        print,'horifoc =',self.horizontal
        print,'horcol =',self.collimatorh
        print,'vercol =',self.collimatorv
        print,'mosaic =',self.mosaic
        print,'dspacing=',self.dspacing
        print,'instrumentOrient=',instrument_orient
        print,'qetrans=',self.qetrans
        print,'feief=',self.feief
        print,'source=',self.source
        print,'rm='
        print,self.rm
        print,'r0=',self.r0


    endif;report




    return,1
end;init


;###############################################################################
;
;NAME:
;        Rescal2_calculator__define
;
;PURPOSE:
;
;PARAMETERS:
;
;KEYWORDS:
;
;_EXTRA - is used to pass keywords meant for the superclass.
;
;###############################################################################
pro Rescal2_calculator__define,class

;    print,'RESELPS LINE 337'
;    widget_control, event.top, get_uvalue = ptrAll


    class = {   Rescal2_calculator,$
                    mylattice:obj_new(),$               LatticeObject
                    ;                                   ;ORIGINAL RESELPS VARIABLES:
                    horizontal:0,$                      ;horifoc
                    collimatorh:dblarr(4),$
;                    alz=(*ptrAll).collimatorh[0]       ;alz    ;z,0 = Source-Mono
;                    alm=(*ptrAll).collimatorh[1]       ;alm    ;m,1 = Mono-Sample
;                    ala=(*ptrAll).collimatorh[2]       ;ala    ;a,2 = Sample-Ana
;                    al3=(*ptrAll).collimatorh[3]       ;al3    ;3   = Ana-Det
                    collimatorv:dblarr(4),$
;                    vbetz=(*ptrAll).collimatorv[0]     ;vbetz
;                    vbet1=(*ptrAll).collimatorv[1]     ;vbet1
;                    vbet2=(*ptrAll).collimatorv[2]     ;vbet2
;                    vbet3=(*ptrAll).collimatorv[3]     ;vbet3
                    mosaic:dblarr(4),$
;                    etash=(*ptrAll).mosaic[0]          ;etash  sample horizontal
;                    etasv=(*ptrAll).mosaic[1]          ;etasv  sample vertical
;                    etam=(*ptrAll).mosaic[2]           ;etam   monochromator
;                    etaa=(*ptrAll).mosaic[3]           ;etaa   analyzer
                    dspacing:dblarr(2),$
;                    tm=2.0*!pi/(*ptrAll).dspacing[0]   ;tm = 2*!pi*monoDspacing
;                    ta=2.0*!pi/(*ptrAll).dspacing[1]   ;ta = 2*!pi*anaDspacing
                    instrument_orient:intarr(2),$
;                    epm=(*ptrAll).orient[0]            ;epm
;                    ep=(*ptrAll).orient[1]             ;ep
                    qetrans:dblarr(5),$                 ;Q AND ENERGY TRANSFER (QX,QY,E,QZ)???
                    feief:0,$
                    source:0,$
                    rm:dblarr(4,4),$                    ;rm
                    r0:0.0d}                            ;r0


end;rescal2_calculator__define

