;###############################################################################
function hkl_to_Qxy,h,k,l,a,b,c,alpha,beta,gamma,o1,o2,$
                    zhat_hkl=zhat_hkl,zhat_orient=zhat_orient,$
                    invang=invang,hklunits=hklunits



    if n_elements(invang) eq 0 then invang = 0
    if n_elements(hklunits) eq 0 then hklunits = 0

    ;DEFAULT TO Å^-1
    if invang eq 0 and hklunits eq 0 then invang = 1

;1) ;GET RECIPROCAL LATTICE VECTORS
    val = recip_lattice(a,b,c,alpha,beta,gamma)
    b1    = val.b123.b1
    b2    = val.b123.b2
    b3    = val.b123.b3

;2) ;GET Q IN Å^-1
    Q = h*b1 + k*b2 + l*b3

;3) ;GET o1vec,o2vec IN 1/ang
;THE FOLLOWING SETS THE AXES ALONG o1,o2
    o1vec = (o1[0]*b1 + o1[1]*b2 + o1[2]*b3)
    o2vec = (o2[0]*b1 + o2[1]*b2 + o2[2]*b3)
    mago1vec = sqrt(dotp(o1vec,o1vec))
    mago2vec = sqrt(dotp(o2vec,o2vec))

;4) ; GET UNIT VECTORS ALONG o1,o2
    ;CREATE UNIT VECTORS ALONG THE o1,o2 DIRECTIONS ---  THIS SHOULD NOT BE USED IN
    xhat_o1 = o1vec/mago1vec                            ;GENERAL SINCE THEY ARE ALONG
    yhat_o2 = o2vec/mago2vec                            ;THE ORIENTATION VECTOR DIRECTIONS
                                                        ;AND MAY NOT BE ORTHOGONAL.

;5) DEFINE  xhat,yhat, UNIT VECTORS IN THE SCATTERING PLANE
;           zhat       PERPENDICULAR TO THE SCATTERING PLANE
    xhat = xhat_o1
    xhat = xhat/sqrt(dotp(xhat,xhat))
    zhat = crossp(xhat_o1,yhat_o2)
    zhat = zhat/sqrt(dotp(zhat,zhat))
    yhat = crossp(zhat,xhat)
    yhat = yhat/sqrt(dotp(yhat,yhat))

;    Qz = 0.0
;;6) CHECK THAT THE INPUT Q IS IN THE SCATTERING PLANE
;    ;NOTE: THERE IS NO REASON THAT I COULDN'T SIMPLY USE ALL 3 COMPONENTS
;    ;JUST PROVIDE A WARNING HERE.
;    testvol = clearjunk(dotp(Q,crossp(xhat,yhat)))
;    if testvol ne 0.0 then begin
;
;        ;if h eq 3 and k eq 3 and l eq 3 then print,'Problem at 1193'
;;        print,'Q=',Q
;;        print,'o1vec=',o1vec
;;        print,'o2vec=',o2vec
;        ;void = dialog_message('Q IS NOT IN THE SCATTERING PLANE. ONLY x,y COMPONENTS CORRECT');,/warning)
;        Qz = nan()
;    endif


;RETURN THE Qz COMPONENT AS WELL.  LET THE NEXT GUY CHECK THAT IT IS
;IN THE SCATTERING PLANE WITH Qz eq 0!!!!!!!!!!!!
    Qz = dotp(Q,zhat)

    dummy = calc_Q(1,1,1)
;Q
    array = [dotp(Q,xhat),dotp(Q,yhat),Qz]
    Qxy = clearjunk(array)
xhat_hkl = LK_cramer(b1,b2,b3,xhat)
yhat_hkl = LK_cramer(b1,b2,b3,yhat)
zhat_hkl = LK_cramer(b1,b2,b3,zhat)
zhat_orient = round(zhat_hkl/min(abs(zhat_hkl[where(zhat_hkl ne 0)])))
    if invang eq 1 then begin
        return,Qxy
    endif else begin
        ;hkl units
        ;THESE AXES ARE THE PRIMARY AXES IN spurion
        ;THEY ARE IN UNITS OF hkl, i.e. DIMENSIONLESS!!!
        Qxy[0] = Qxy[0]/mago1vec

        ;GET ANGLE OF o2 WRT X AXIS
        o2theta = acos(dotp(o2vec,xhat)/mago2vec)
        ;CALCULATE THE COMPONENT OF o2 ALONG Y
        yscale  = mago2vec/sin(o2theta)
        ;CREATE ALTERNATE Y AXIS
        Qxy[1] = Qxy[1]/yscale

        Qxy[2] = Qz  ;UNITS NOT CLEAR TO ME


        return,Qxy
    endelse

end;hkl_to_Qxy
