;####################################################################
pro alewald,g,hklal,cu=cu,mo=mo

    ;BOTH Al AND Cu ARE FCC, SO ALL OF THE CALCULATIONS BELOW APPLY.
    ;THE ONLY CHANGE IS THE LATTICE CONSTANT:
    ;
    ;www.webelements.com
    ;   Al  4.0495Å
    ;   Cu  3.6149Å
    ;   Mo  3.147Å  Mo is BCC AND WILL RAISE THE APPROPRIATE SWITCH BELOW.

    ;g and hklal ARE CALCULATED BELOW AND PASSED BACK TO CALLING PROCEDURE
    ;VIA THE ARGUMENTS ABOVE.
    ;
    ;TRANSLATED DIRECTLY FROM FORTRAN BY LRK 01/06/06
    ;ADDED Cu AND Mo RINGS. 08/16/06
    ;

      n=13
      G = dblarr(n)

      if n_elements(cu) eq 0 then cu = 0
      if n_elements(mo) eq 0 then mo = 0

      if mo ne 0 then bcc = 1 else bcc = 0

;     lattice constant.
      if cu ne 0 then begin
        a = 3.6149      ;Cu
      endif else begin
        if mo ne 0 then begin
            a = 3.147   ;Mo
        endif else begin
            ;DEFAULT TO Al
            a = 4.04096 ;Al   FROM ORIGINAL alewald CODE
        endelse
      endelse

      pi = !PI
      magq = 0.0
      dum = 0.0
      qdiff = dblarr(n)
      i = 0
      h = intarr(3,n)
      hklal = intarr(3,n)
;  Read in the aluminum peaks to be considered.
;
          h[*,0]=[1,1,1]
          h[*,1]=[2,0,0]
          h[*,2]=[2,2,0]
          h[*,3]=[3,1,1]
          h[*,4]=[2,2,2]
          h[*,5]=[4,0,0]
          h[*,6]=[3,3,1]
          h[*,7]=[4,2,0]
          h[*,8]=[4,2,2]
          h[*,9]=[3,3,3]
          h[*,10]=[4,4,0]
          h[*,11]=[5,3,1]
          h[*,12]=[ 6,0,0]

          if bcc eq 1 then begin
              h[*,0]=[1,1,0]
              h[*,1]=[2,0,0]
              h[*,2]=[2,1,1]
              h[*,3]=[2,2,0]
              h[*,4]=[3,1,0]
              h[*,5]=[2,2,2]
              h[*,6]=[3,2,1]
              h[*,7]=[4,0,0]
              h[*,8]=[3,3,0]
              h[*,9]=[4,2,0]
              h[*,10]=[3,3,2]
              h[*,11]=[4,2,2]
              h[*,12]=[4,3,1]
          endif

      pi=!PI;4.0*atan(1.0)

;  Calculate the FCC G's.
      for i=0,n-1 do begin
         G[i]=(2*pi/a)*sqrt(float(h[0,i]^2 + h[1,i]^2 + h[2,i]^2))
         for j=0,2 do begin
            hklal[j,i]=h[j,i]
         endfor
      endfor

      return
end;alewald
;####################################################################
function calc_k,d,a1;,a2,d

    ;CALCULATE THE MAGNITUDE OF k FROM THE MONO ANGLE AND d SPACING

    ;a1 in deg
    ;d in ang


    ;CONVERT TO RADIANS
    a = !dtor*a1

    l = 2.0*d*sin(a)

    k = 2.0*!PI/l
    return,k

end;calc_k
;####################################################################
function calc_kmag_from_E,E
    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)
    kmag = sqrt(E/hbar_sq_over2m)

    return,kmag

end;calc_kmag_from_E
;####################################################################
function calc_lambda_from_e,E

    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)


    lambda = 2*!PI/sqrt(E/hbar_sq_over2m)

    return,lambda

end;calc_lambda_from_e
;####################################################################
function calc_E_from_lambda,lambda

    hbar_sq_over2m = 2.072142

    k = 2*!PI/lambda

    E = hbar_sq_over2m*k^2

    return,E
end;calc_E_from_lambda
;####################################################################
function calc_E_from_kvec,k


    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)



    E = hbar_sq_over2m*(k[0]^2 + k[1]^2 + k[2]^2)
    return,E

end;calc_Ei_from_kvec
;####################################################################
function calc_E_from_kxy,k


    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)



    E = hbar_sq_over2m*(k[0]^2 + k[1]^2)
    return,E

end;calc_Ei_from_kxy
;####################################################################
function calc_E_from_k,k

    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)

    E = hbar_sq_over2m*k^2
    return,E

end;calc_Ei
;####################################################################

function calc_lambdamax_from_d,d
    l = 2*d;*sin(!pi/2)
    return,l
end;calc_lambdamin_from_d
;####################################################################
function calc_Emin_from_d,d


    hbar_sq_over2m = 2.072142 ;meV.Ang^2  (Zheludev value)

    l = 2*d;*sin(!pi/2)

    k = 2*!PI/l

    E = hbar_sq_over2m*k^2
    return,E

end;calc_Emax_from_d
;####################################################################
function calc_lambda_from_theta,d,theta

    return,(d*sin(!dtor*theta))

end;calc_lambda_from_theta
;####################################################################
function calc_k_from_theta,d,theta

    k = 2*!PI/(d*sin(!dtor*theta))
    return,k

end;calc_ki_from_theta
;####################################################################
function calc_kpara,ki

    kpara = [ki,0]

end;calc_kpara
;####################################################################
function calc_kf,kf,A4

    A4 = !dtor*A4

    kperp = kf*[cos(A4),sin(A4)]

end;calc_kf
;####################################################################
function calc_Qmax,ei,ef,A4

    if n_elements(A4) eq 0 then A4 = 180

    ang = !dtor*A4

    ki = calc_kmag_from_E(ei)
    kf = calc_kmag_from_E(ef)



    Qmax = sqrt(ki^2 + kf^2 - 2*ki*kf*cos(ang))

    return,Qmax

end;calc_Qmax

;;####################################################################
;function LK_determ,M
;    dM   = M[0,0]*(M[1,1]*M[2,2]-M[2,1]*M[1,2]) + $
;           M[1,0]*(M[2,1]*M[0,2]-M[0,1]*M[2,2]) + $
;           M[2,0]*(M[0,1]*M[1,2]-M[1,1]*M[0,2])
;    return,DM
;end;LK_determ
;;####################################################################
;function LK_cramer,b1,b2,b3,yhat
;
;;COEFFICIENTS
;M = transpose([[b1],[b2],[b3]])
;;M = [[b1],[b2],[b3]]
;
;;SUB RESULTS (yhat) IN TO COLUMNS
;;print,'A1,A2,A3:'
;MA1 = M
;MA1[0,*] = yhat;BA
;;print,'MA1=',MA1
;MA2 = M
;MA2[1,*] = yhat;BA
;;print,'MA2=',MA2
;MA3 = M
;MA3[2,*] = yhat
;;print,'MA3=',MA3
;
;
;;CALCULATE DETERMINANTS
;dM   = LK_determ(M)
;dMA1 = LK_determ(MA1)
;dMA2 = LK_determ(MA2)
;dMA3 = LK_determ(MA3)
;
;;print,'dA1/dA=',dA1/dM
;;print,'dA2/dA=',dA2/dM
;;print,'dA3/dA=',dA3/dM
;
;A1 = dMA1/dM
;A2 = dMA2/dM
;A3 = dMA3/dM
;
;;print,'[A1,A2,A3]=',A1,A2,A3
;;print,'yhat=',yhat
;;print,'A1*b1+A2*b2+A3*b3=',A1*b1+A2*b2+A3*b3
;
;return,[A1,A2,A3]
;
;
;end;LK_cramer


;;####################################################################
;function clearjunk,x
;
;    print,'clearjunk in file calc_q.pro'
;    for i=0,n_elements(x)-1 do begin
;        if abs(x[i]) lt 1.0e-8 then x[i] = 0.0
;    endfor;i
;
;    return,x
;end;clearjunk
;####################################################################
function get_recip_lattice,a,b,c,alpha,beta,gamma,hkl

    vals = recip_lattice(a,b,c,alpha,beta,gamma)

    return,hkl[0]*vals.b123unit.b1unit + $
           hkl[1]*vals.b123unit.b2unit + $
           hkl[2]*vals.b123unit.b3unit

end;get_recip_lattice
;####################################################################
function make_hkl,d


    if n_elements(d) eq 0 then d = 9

    hkl = intarr(3,d^3)

    for i=0,d-1 do begin
        for j=0,d-1 do begin
            for k=0,d-1 do begin
                hkl[*,i+j*d+k*d^2] = [i-(d-1)/2,j-(d-1)/2,k-(d-1)/2]
            endfor;k
        endfor;j
    endfor;i

    return,hkl
end;make_hkl
;####################################################################
function QE,ki,kf,nvec

        ;CALCULATE ALL OF THE Q AND E VALUES FOR ONE POINT ONLY.
        Q = ki - kf
        if n_elements(nvec) eq 0 then nvec = 5

        kiarr = dblarr(3,nvec)
        kfarr = dblarr(3,nvec)
        Qarr  = dblarr(3,nvec,nvec)
        QMagarr = dblarr(nvec,nvec)
        Earr  = dblarr(nvec,nvec)
        for i=1,nvec do begin
            kiarr[*,i-1] = i*ki
            kfarr[*,i-1] = i*kf
        endfor;i
        for i=0,nvec-1 do begin
            for j=0,nvec-1 do begin
                ;CALCULATE THE Q'S FOR THE HARMONICS
                Qarr[*,i,j] = kiarr[*,i] - kfarr[*,j]
                QMagarr[i,j] = sqrt(dotp(Qarr[*,i,j],Qarr[*,i,j]))
                ;CALCULATE THE ENERGIES FOR THE HARMONICS
                Earr[i,j] = calc_E_from_kvec(kiarr[*,i]) - calc_E_from_kvec(kfarr[*,j])

            endfor;j
        endfor;i

   return,{Q:Q,kiarr:kiarr,kfarr:kfarr,Qarr:Qarr,QMagarr:QMagarr,Earr:Earr}
end;QE
;####################################################################
function getkikf,Ei,Ef,Qxy,closed


    ;RETURNS ki,kf IN x,y COORDINATES IN UNITS OF Å^-1

    ;I WILL GET E, Ei or Ef FIXED, AND DETERMINE WHICH IS FIXED.
    thetaq = atan(Qxy[1]/Qxy[0])

    Q = sqrt(dotp(Qxy,Qxy))

    ki = calc_kmag_from_E(Ei)
    kf = calc_kmag_from_E(Ef)

    ;CHECK THAT SCATTERING TRIANGLE CAN CLOSE WITH THE LAW OF COSINES
    ;Q^2 = ki^2 + kf^2 - 2*ki*kf*cos(angQ)
    check = (ki^2 + kf^2 - Q^2)/(2*ki*kf)
;print,'In getkikf:'
;print,'check=',check
;print,'Ei,f=',Ei,Ef
;help,/traceback
    if abs(check) gt 1 or (finite(check,/nan) ne 0) then begin
        closed = 0
;        print,'SCATTERING TRIANGLE CANNOT CLOSE WITH THE VALUES:'
;        print,'cos(angQ)=',check
;        print,'Ei=',Ei
;        print,'Ef=',Ef
;        print,'ki=',ki
;        print,'kf=',kf
;        print,'Q =['+strtrim(string(Qxy[0]),2)+','+strtrim(string(Qxy[1]),2)+']'
;        help,/traceback
        kivec = double(['NaN','NaN','NaN'])
        kfvec = double(['NaN','NaN','NaN'])

    endif else begin
        closed = 1
        ;kf^2 = ki^2 + Q^2 -2kiQcos(thetaQ-thetai)
        ;thetaQ-thetai = acos((ki^2+Q^2-kf^2)/2kiQ)
        ;thetaQ-acos(arg)=thetai
        arg = (ki^2 + Q^2 - kf^2)/(2*ki*Q)
        thetai = thetaq - acos(arg)

        kivec = ki*[cos(thetai),sin(thetai),0.0]
        kfvec = kivec - Qxy
    endelse

    return,{ki:kivec,kf:kfvec}

end;getkikf
;####################################################################
function getkikf_from_efixed_etran,Efixed,Etran,Qxy,init=init,final=final

    if n_elements(init) eq 0  then init = 0
    if n_elements(final) eq 0 then final = 0

    if init eq 0 and final eq 0 then final = 1

;print,'init=',init
;print,'final=',final
    ;E = Ei - Ef
    E = Etran
    if final eq 1 then begin
        Ef = Efixed
        Ei = E + Ef
    endif
    if init eq 1 then begin
        Ei = Efixed
        Ef = Ei - E
    endif

;    print,'Ei=',Ei
;    print,'Ef=',Ef
;    print,'E=',Ei-Ef
;    print,'Qxy=',Qxy
    return,getkikf(Ei,Ef,Qxy)


end;getkikf_from_efixed_etran
;####################################################################
function getkikf_from_Qxy,Efixed,Etran,Qxy,init=init,final=final
    return,getkikf_from_efixed_etran(Efixed,Etran,Qxy,init=init,final=final)
end;get_kikf_from_Qxy
;####################################################################
function angstrom
    return,'Å'
end;angstrom
;####################################################################
function plusminus
    return,'±'
end;plusminus
;###############################################################################
function recip_kikf,o1,o2,a,b,c,alpha,beta,gamma,Q,Efix,E,init=init,final=final,qobj=qobj

    nan = double('NaN')
    result = {Q:Q,ki:[nan,nan,nan],kf:[nan,nan,nan]}

    if final eq 1 then begin
        Ef = Efix
        Ei = E + Ef
    endif
    if init eq 1 then begin
        Ei = Efix
        Ef = Ei - E
    endif



;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/ang
    Qorig = Q   ;SAVE hkl VALUES FOR Q
    Q = Q[0]*b1 + Q[1]*b2 + Q[2]*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))


;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
    testvol = clearjunk(dotp(Q,crossp(xhat,yhat)))
    if testvol ne 0.0 then begin
        print,'Q=',Q
        print,'o1vec=',o1vec
        print,'o2vec=',o2vec
        print,'Q IS NOT IN THE SCATTERING PLANE!!!!!!!'
        print,'QUITTING.'
        void = dialog_message('Q IS NOT IN THE SCATTERING PLANE.',/error)
        return,result
    endif


;Q
    Qxy = [dotp(Q,xhat),dotp(Q,yhat),0.0]   ;THIS ASSUMES THE SCATTERING PLANE CHECK IS DONE!!!

;;;;;
;;;;;;FINALLY TRANSFER kixy,kfxy TO hkl SPACE

;FIRST FIND yhat IN TERMS OF b1,b2,b3
;
;USE CRAMER'S RULE TO SOLVE:
;[[b1],[b2],[b3]]#[A1,A2,A3] = yhat  --- FIND A1,A2,A3 COEFFICIENTS

;7) GET THE hkl VALUES FOR xhat,yhat SO INDICES DESCRIBE THE UNIT VECTORS
    xhat_hkl = LK_cramer(b1,b2,b3,xhat)
    yhat_hkl = LK_cramer(b1,b2,b3,yhat)
    zhat_hkl = LK_cramer(b1,b2,b3,zhat)

;print,'yhat_hkl=',yhat_hkl
;print,'yhat=',yhat_hkl[0]*b1 + yhat_hkl[1]*b2 + yhat_hkl[2]*b3
;print,'o2=',o2

;THE NEXT CHECK WORKS:
;    print,'IN xyz  |b1| = ',sqrt(total([b1x,b1y,b1z]^2))
;    print,'IN b123 |b1| = ',sqrt(dotp(b1,b1))


;NOW ki,kf IN XY SPACE (Qxy SETS THE AXES.)

;NOTE THAT IMPLICITLY ki_z and kf_z are zero!!!!!!!
    ;GET ki,kf IN XY PLANE

;8) GET ki,kf IN X,Y,Z COORDINATES
    kikf = getkikf_from_efixed_etran(Efix,E,Qxy,init=init,final=final)
    ki = kikf.ki
    kf = kikf.kf




    ;VERIFY THAT SCATTERING TRIANGLE CLOSED
    nani = total(finite(ki,/nan))
    nanf = total(finite(kf,/nan))

    if total([nani,nanf]) ne 0 then begin
        void = dialog_message('SCATTERING TRIANGLE CANNOT CLOSE WITH ENTERED VALUES.',/error)
        return,result
    endif else begin

;9) GET ki,kf IN HKL VALUES
        kix_hkl = ki[0]*xhat_hkl
        kiy_hkl = ki[1]*yhat_hkl
        kiz_hkl = 0.0*zhat_hkl  ;NO z COMPONENT SINCE THERE WAS AN EARLIER CHECK TO SEE Q AND THUS
                                ;ki,kf ARE IN THE SCATTERING PLANE.

        ki_hkl = kix_hkl + kiy_hkl + kiz_hkl

        kfx_hkl = kf[0]*xhat_hkl
        kfy_hkl = kf[1]*yhat_hkl
        kfz_hkl = 0.0*zhat_hkl

        kf_hkl = kfx_hkl + kfy_hkl + kfz_hkl

        result.ki = ki_hkl
        result.kf = kf_hkl
        result.Q  = Qorig

;AT THIS POINT, CREATE, ki,kf, Q objects
;BUT WHICH VALUES DO I WANT TO PASS?  Qx,Qy,kix,kiy,kfx,kfy,Ei,Ef,E


;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
;        qobj = obj_new('spurion2_q_data',$
;                            Qx = Qxy[0],Qy = Qxy[1],$
;                            kix= ki[0], kiy = ki[1],$
;                            kfx= kf[0], kfy = kf[1],$
;                            Energy=E)
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;



        return,result
    endelse



end;recip_kikf_in_hkl
;####################################################################
pro test_recip,o1,o2,a,b,c,alpha,beta,gamma,ki=ki,kf=kf,qwid=qwid,ewid=ewid,scanwid=scanwid,harmonics=harmonics

    ;INPUTS:
    ;a,b,c,alpha,beta,gamma
    ;o1,o2
    ;hkl,strf
    ;ki,kf IN TERMS OF hkl VALUES

    if n_elements(qwid) eq 0 then qwid = 0
    if n_elements(ewid) eq 0 then ewid = 1
    if n_elements(scanwid) eq 0 then scanwid = 2
    if n_elements(harmonics) eq 0 then harmonics = 0


    hkl = make_hkl(31)

;STEP 1 GET RECIP LATTICE VECTORS FROM abc,abg
    val = recip_lattice(a,b,c,alpha,beta,gamma)
    b1    = val.b123.b1
    b2    = val.b123.b2
    b3    = val.b123.b3

    ;CALCULATE ki,kf,Q IN 1/ang UNITS
    if n_elements(ki) eq 3 then begin
        ki = ki[0]*b1 + ki[1]*b2 + ki[2]*b3
        Ei = calc_E_from_kvec(ki)
        ;print,Ei
    endif
    if n_elements(kf) eq 3 then begin
        kf = kf[0]*b1 + kf[1]*b2 + kf[2]*b3
        Ef = calc_E_from_kvec(kf)
        ;print,Ef
    endif

    if n_elements(ki) eq 3 and n_elements(kf) eq 3 then begin
        nvec = 5
        qevals = QE(ki,kf,nvec)

        Q       = qevals.Q
;        print,'_________________________________________________________'
;        print,'In test_recip b1=',b1
;        print,'In test_recip b2=',b2
;        print,'In test_recip b3=',b3
;        print,'In test_recip Q= ',Q
;        print,'Ei=',Ei
;        print,'Ef=',Ef
;        print,'E=',Ei-Ef
;        print,'ki_b123=',ki
;        print,'kf_b123=',kf
;        print,'kimag=',sqrt(dotp(ki,ki))
;        print,'kfmag=',sqrt(dotp(kf,kf))


        kiarr   = qevals.kiarr
        kfarr   = qevals.kfarr
        Qarr    = qevals.Qarr
        QMagarr = qevals.QMagarr
        Earr    = qevals.Earr
    endif


;IT SEEMS LIKE IT SHOULD BE SOMTHING LIKE

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

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


;BUT I REALLY WANT A PAIR OF AXES WHERE
    ;o1 along x
    ;y perp to x regardless of whether o2 perp to o1
    ;z perp to x,y
;SO:
;    xhat = xhat_o1
;    zhat = crossp(xhat_o1,yhat_o2)
;    yhat = crossp(zhat,xhat)
    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))

;print,'xhat=',xhat
;print,'yhat=',yhat
;print,'zhat=',zhat
;HERE THE UNITS ARE SAME AS THOSE OF b1,b2,b3 ---- 1/Angstrom
    ;So adding Al rings is easy now

    szhkl = size(hkl)
    nq = szhkl[2]
    strf  = findgen(nq)

    title = '';'Qxy in o1,o2 plane'
    xtitle = 'X ['+angstrom()+'^-1]  || ['+ $
                strtrim(string(o1[0],format='(f4.1)'),2)+','+ $
                strtrim(string(o1[1],format='(f4.1)'),2)+','+ $
                strtrim(string(o1[2],format='(f4.1)'),2)+']';o1[0]*b1+o1[1]*b2+o1[2]*b3'
    ytitle = 'Y ['+angstrom()+'^-1]  || Z x X';perp o1.[b1,b2,b3],coplanar w/o1.[b1,b2,b3]&o2.[b1,b2,b3]'

    overalltitle = 'Plane perp to Z = O1 x O2'

    ;SET THE MOST USEFUL COLORS
    black = 0L
    red   = 255L
    green = 256L*255L
    blue  = 256L^2*255L
    white = red+green+blue

    if qwid eq 0 then begin
        window,0,xsize=500,ysize=500,xpos=600,ypos=200
    endif else begin
        wset,qwid
    endelse

;    xrange = [-10,10]
;    yrange = [-10,10]
;
;    xrange = [-1,1]
;    yrange = [-1,1]

;GET RANGES
    allgx = [0.0]
    allgy = [0.0]
    for i=0,nq-1 do begin
        h = hkl[0,i]
        k = hkl[1,i]
        l = hkl[2,i]
        g    = h*b1 + k*b2 + l*b3
        gx = dotp(g,xhat)
        gy = dotp(g,yhat)


        if dotp(G,crossp(o1vec,o2vec)) eq 0.0 then begin
            if abs(h) lt 5 and abs(k) lt 5 and abs(l) lt 5 then begin
                allgx = [allgx,gx]
                allgy = [allgy,gy]
            endif

        endif;IN SCATTERING PLANE
    endfor;i

    allgx = [allgx,2*dotp(kiarr[*,0],xhat)]
    allgy = [allgy,2*dotp(kiarr[*,0],yhat)]
    allgx = [allgx,2*dotp(kfarr[*,0],xhat)]
    allgy = [allgy,2*dotp(kfarr[*,0],yhat)]

    xrange = [min(allgx),max(allgx)]
    yrange = [min(allgy),max(allgy)]

    plot,xrange,yrange,background = white,color=0L,/nodata,$
         subtitle=title,xtitle=xtitle,ytitle=ytitle,xstyle=4,ystyle=4,$
         xmargin=[8,8],ymargin=[4,4];,/isotropic
xyouts,0.25,0.97,overalltitle,color=black,/normal

    x0 = !x
    y0 = !y
;    print,!x.crange

    ;I CAN SWAP THE AXES WHENEVER I WANT BY SWITCHING xaxis=0 to xaxis=1
    ;AND VICE VERSA.

    ;THESE AXES ARE THE ALTERNATE AXES IN spurion.
    ;THEY ARE IN INVERSE ANGSTROM UNITS
    AXIS, XAXIS=1, XRANGE = (!X.CRANGE), XSTYLE = 1, $
          XTITLE = xtitle,color=0L,xcharsize=0.75
    AXIS, YAXIS=1, YRANGE = (!Y.CRANGE), YSTYLE = 1, $
          YTITLE = ytitle,color=0L,ycharsize=0.65


    xtitle0 = 'X [hkl units]  || ['+ $
                strtrim(string(o1[0],format='(f4.1)'),2)+','+ $
                strtrim(string(o1[1],format='(f4.1)'),2)+','+ $
                strtrim(string(o1[2],format='(f4.1)'),2)+']';o1[0]*b1+o1[1]*b2+o1[2]*b3'
    ytitle0 = 'Y [hkl units]  || Z x X';perp o1.[b1,b2,b3],coplanar w/o1.[b1,b2,b3]&o2.[b1,b2,b3]'

    ;THESE AXES ARE THE PRIMARY AXES IN spurion
    ;THEY ARE IN UNITS OF hkl, i.e. DIMENSIONLESS!!!
    AXIS, XAXIS=0, XRANGE = (!X.CRANGE)/mago1vec, XSTYLE = 1, $
          XTITLE = xtitle0,color=0L

    ;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
    AXIS, YAXIS=0, YRANGE = (!Y.CRANGE)/yscale, YSTYLE = 1, $
          YTITLE = ytitle0,color=0L

    thick = 2

    !x = x0
    !y = y0

    ;DRAW Q WITH Ki,Kf HARMONICS
    if n_elements(Q) eq 3 then begin

        ;ki,kf HARMONICS
        for i=0,nvec-1 do begin

            if i eq 0 then begin
                xio = dotp(kiarr[*,i],xhat)
                yio = dotp(kiarr[*,i],yhat)

                ;USE -kf FOR PLOTTING
                xfo = -dotp(kfarr[*,i],xhat)
                yfo = -dotp(kfarr[*,i],yhat)

;                arrow,[0],[0],[dotp(kiarr[*,i],xhat)],[dotp(kiarr[*,i],yhat)],color=blue+red,/data
;                arrow,[0],[0],[dotp(kfarr[*,i],xhat)],[dotp(kfarr[*,i],yhat)],color=blue+green,/data

                arrow,[0],[0],xio,yio,color=blue+red,/data,thick=thick
                arrow,[xio],[yio],xio+xfo,yio+yfo,color=blue+green,/data,thick=thick
            endif;i eq 0
;            if harmonics eq 1 then begin
;                ;print,'HEY HARMONICS!!!'
;                xio = dotp(kiarr[*,i],xhat)
;                yio = dotp(kiarr[*,i],yhat)
;
;                ;USE -kf FOR PLOTTING
;                xfo = -dotp(kfarr[*,i],xhat)
;                yfo = -dotp(kfarr[*,i],yhat)
;                oplot,[xio+xfo],$
;                      [yio+yfo],$
;                      psym=2,$
;                      color=red,$
;                      symsize=5;green;red+green
;                oplot,[xfo+xio],$
;                      [yfo+yio],$
;                      psym=2,$
;                      color=red,$
;                      symsize=5;green;color=red+green
;
;;                arrow,[0],[0],[dotp(kiarr[*,i],xhat)],[dotp(kiarr[*,i],yhat)],color=blue+red,/data
;;                arrow,[0],[0],[dotp(kfarr[*,i],xhat)],[dotp(kfarr[*,i],yhat)],color=blue+green,/data
;            endif;harmonics
        endfor;i

        if harmonics eq 1 then begin
            ;print,'HARMONICS!!!!!'
            for i=0,nvec-1 do begin
                xio = dotp(kiarr[*,i],xhat)
                yio = dotp(kiarr[*,i],yhat)
                xfo = -dotp(kfarr[*,i],xhat)
                yfo = -dotp(kfarr[*,i],yhat)
                for j=0,nvec-1 do begin
                    xfoj=-dotp(kfarr[*,j],xhat)
                    yfoj=-dotp(kfarr[*,j],yhat)
                    xioj=dotp(kiarr[*,j],xhat)
                    yioj=dotp(kiarr[*,j],yhat)

                    color=long(red*double(i)/double(nvec-1))+long(blue*double(j)/double(nvec-1))

                    oplot,[xio+xfoj],$
                          [yio+yfoj],$
                          psym=2,$
                          color=color;blue;green;red+green
                    oplot,[xfo+xioj],$
                          [yfo+yioj],$
                          psym=2,$
                          color=color;blue;green;red+green
;                    arrow,[xio],[yio],[xio+dotp(kfarr[*,j],xhat)],$
;                          [yio+dotp(kfarr[*,j],yhat)],color=blue+green,/data
;                    arrow,[xfo],[yfo],[xfo+dotp(kiarr[*,j],xhat)],$
;                          [yfo+dotp(kiarr[*,j],yhat)],color=blue+red,/data
                endfor;j
            endfor;i
        endif;harmonics

        ;Q
        arrow,[0],[0],[dotp(Qarr[*,0,0],xhat)],[dotp(Qarr[*,0,0],yhat)],color=blue,/data,thick=thick

    endif

    ;PLOT THE RECIPROCAL LATTICE
    for i=0,nq-1 do begin

        h = hkl[0,i]
        k = hkl[1,i]
        l = hkl[2,i]

        g    = h*b1 + k*b2 + l*b3
;        ghkl = [h,k,l]

        gx = dotp(g,xhat)
        gy = dotp(g,yhat)

        ;THIS CONVERSION IS COVERED BY THE axis CALLS ABOVE!!!
        ;AND THIS IS NOT CORRECT HERE!!!
;        gxhkl = dotp(ghkl,o1hat)
;        gyhkl = dotp(ghkl,o2hat)

        ;PLOT ALL G's IN THE SCATTERING PLANE
        if dotp(G,crossp(o1vec,o2vec)) eq 0.0 then begin

            wset,qwid;0
            !x = x0 & !y = y0
            oplot,[gx],[gy],color=255L+(256L^2)*bytscl(alog(strf[i]+1)),psym=4

        endif;IN SCATTERING PLANE
    endfor;i

    ;IF I GET THIS WORKING RIGHT, THEN I CAN CONVERT
    ;Al POWDER RINGS AND Ki,Kf INTO Q SPACE.

    alewald,al,hklal
    ;help,al

    theta = 2.0*!PI*findgen(361)/360.0
    for i=0,n_elements(al)-1 do begin
        x = al[i]*cos(theta)
        y = al[i]*sin(theta)

        wset,qwid;0
        !x = x0 & !y = y0
        oplot,x,y,psym=0,linestyle=i mod 5,color=256L*255L

    endfor;i


    ;PLOT ENERGY v. Q
    qmag = sqrt(dotp(Q,Q))
    if ewid eq 1 then begin
        window,1
    endif else begin
        wset,ewid
    endelse

    ;plot,[qmag,qmag],[min(earr),max(earr)],color=black,background=white,/nodata
;    plot,[Q,Q],[min(earr),max(earr)],color=black,background=white,/nodata,$
;        xtitle='Q',ytitle='dE'
;    for i=0,nvec-1 do begin
;        for j=0,nvec-1 do begin
;                oplot,[qmag],[earr[i,j]],psym=4,$
;                        color=long(red*double(i)/double(nvec))+long(blue*double(i)/double(nvec))
;        endfor;j
;    endfor;i

;SIMULATE ELEMENTS OF A SCAN
    if n_elements(ki) eq 3 and n_elements(kf) eq 3 then begin
        nvec = 5
        qevalsm1 = QE(ki-[0.1*ki[0],0,0],kf,nvec)

        Qm1       = qevalsm1.Q
        kiarrm1   = qevalsm1.kiarr
        kfarrm1   = qevalsm1.kfarr
        Qarrm1    = qevalsm1.Qarr
        QMagarrm1 = qevalsm1.QMagarr
        Earrm1    = qevalsm1.Earr
        qmagm1 = sqrt(dotp(Qm1,Qm1))
    endif
;    for i=0,nvec-1 do begin
;        for j=0,nvec-1 do begin
;                oplot,[qmag],[earrm1[i,j]],psym=4,$
;                        color=long(red*double(i)/double(nvec))+long(blue*double(i)/double(nvec))
;        endfor;j
;    endfor;i
;SIMULATE ELEMENTS OF A SCAN

;USE THE FOLLOWING FUNCTION TO GET ki,kf FOR EACH POINT WHEN GIVEN THE Q
;
;getkikf_from_efixed_etran,Efixed,Etran,Qxy,kf_fixed=kf_fixed,ki_fixed=ki_fixed
;
;
    if n_elements(ki) eq 3 and n_elements(kf) eq 3 then begin
        nvec = 5
        qevalsp1 = QE(ki+[0.1*ki[0],0,0],kf,nvec)

        Qp1       = qevalsp1.Q
        kiarrp1   = qevalsp1.kiarr
        kfarrp1   = qevalsp1.kfarr
        Qarrp1    = qevalsp1.Qarr
        QMagarrp1 = qevalsp1.QMagarr
        Earrp1    = qevalsp1.Earr
        qmagp1 = sqrt(dotp(Qp1,Qp1))
    endif
;    for i=0,nvec-1 do begin
;        for j=0,nvec-1 do begin
;                oplot,[qmag],[earrp1[i,j]],psym=4,$
;                        color=long(red*double(i)/double(nvec))+long(blue*double(i)/double(nvec))
;        endfor;j
;    endfor;i
;    print,[qmagm1,qmag,qmagp1]

    edata = [earrm1,earr,earrp1]
    qdata = [qmagm1,qmag,qmagp1]
    erange = [min(edata),max(edata)]
    qrange = [min(qdata),max(qdata)]
    plot,qrange,erange,color=black,background=white,/nodata,$
        xtitle='Q',ytitle='dE',title='Hypothetical 3pt scan'

    for i=0,nvec-1 do begin
        for j=0,nvec-1 do begin
;                print,[earrm1[i,j],earr[i,j],earrp1[i,j]]
                oplot,[qmagm1,qmag,qmagp1],[earrm1[i,j],earr[i,j],earrp1[i,j]],psym=-3,$
                        color=long(red*double(i)/double(nvec-1))+long(blue*double(j)/double(nvec-1))
        endfor;j
    endfor;i


end;test_recip
;####################################################################
pro test_recipQEEfix,o1,o2,a,b,c,alpha,beta,gamma,Q,Efix,E,$
                        init=init,final=final,hkl=hkl,$
                        qwid=qwid,ewid=ewid,scanwid=scanwid,$
                        harmonics=harmonics,zhat_hkl=zhat_hkl


    ;INPUTS:
    ;a,b,c,alpha,beta,gamma
    ;o1,o2
    ;hkl,strf

    if n_elements(init)  eq 0 then init  = 0
    if n_elements(final) eq 0 then final = 0
    if n_elements(harmonics) eq 0 then harmonics = 0

    if init eq 0 and final eq 0 then final = 1

    hkl = make_hkl(31)  ;hkl, strf COULD BE INPUT

    result = recip_kikf(o1,o2,a,b,c,alpha,beta,gamma,Q,Efix,E,init=init,final=final)
    ki_hkl = result.ki
    kf_hkl = result.kf




;10) FINALLY PASS ki,kf TO ki,kf PROGRAM
    if total(finite(ki_hkl,/nan)+finite(kf_hkl,/nan)) ne 0 then begin
        return
    endif else begin
        test_recip,o1,o2,a,b,c,alpha,beta,gamma,ki=ki_hkl,kf=kf_hkl,$
                    qwid=qwid,ewid=ewid,scanwid=scanwid,harmonics=harmonics
    endelse

end;test_recipQEEfix
;###############################################################################
function Qxy_to_hkl,Qxy,a,b,c,alpha,beta,gamma,o1,o2

    ;FOR MOUSE INTERACTIONS IN THE SCATTERING PLANE:

    ;GIVEN Qxy THIS GIVES THE hkl VALUES FOR A POINT IN THE SCATTERING
    ;PLANE.  THIS WILL BE EXTREMELY USEFUL WHEN THE USER WANTS TO
    ;SELECT A POINT IN THE SCATTERING PLANE WITH THE MOUSE AND
    ;HAVE THE PROGRAM PRODUCE Q,ki,kf, AND OTHER ITEMS.


    ;ADD A zero IN THE Z COMPONENT IF ONLY x,y, COMPONENTS SPECIFIED.
    if n_elements(Qxy) eq 2 then Q = [Qxy,0.0] else Q = Qxy




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

;3) ; 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.

;4) 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))

    ;NOW GET xhat,yhat IN hkl
;5) GET THE hkl VALUES FOR xhat,yhat SO INDICES DESCRIBE THE UNIT VECTORS
    xhat_hkl = LK_cramer(b1,b2,b3,xhat)
    yhat_hkl = LK_cramer(b1,b2,b3,yhat)
    zhat_hkl = LK_cramer(b1,b2,b3,zhat)

    Qhkl = Q[0]*xhat_hkl + Q[1] * yhat_hkl + Q[2]*zhat_hkl
    return,Qhkl
end;Qxy_to_hkl

;####################################################################
function calc_Q,ki,kf,A4
;print,'RUNNING calc_Q'
    A4deg = A4
;print,'test'


    dum = clearjunk(1)
    dum = clearjunk([1,1])
    dum = clearjunk([1,1,1])

    i = calc_kpara(ki)
    f = calc_kf(kf,A4)

    A4 = !dtor*A4deg

    return, i - f
end;calc_Q


