; $Id: triangle_cam_funq.pro,v 1.3 2005/02/09 16:58:01 dimeo Exp $
;###############################################################################
;
; NAME:
;  TRIANGLE_CAM_FUNQ
;
; PURPOSE:
;  Triangle cam function implementation for the time-to-energy conversion
;  for the HFBS data.
;
; CATEGORY:
;  DAVE, HFBS, data reduction
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function triangle_cam_funQ,npts,f
; Function returns the monochromator velocity profile over
; one period of oscillation of the monochromator (half
; period of the cam rotation).
;
pi = !pi
thetalo = 0.0 & thetahi = 1.0*pi
dtheta = (thetahi-thetalo)/(npts-1)
theta = fltarr(npts)
theta = thetalo+dtheta*findgen(npts)
region1 = where(theta lt 3.0*pi/16.0,n1)
region2 = where(theta ge 3.0*pi/16.0 and theta lt 5.0*pi/16.0,n2)
region3 = where(theta ge 5.0*pi/16.0 and theta lt 0.5*pi,n3)
a = 0.09*256.0/(15.0*pi^2+8.0)
s1 = fltarr(npts) & s2 = fltarr(npts) & s3 = fltarr(npts)
a1 = fltarr(npts) & a2 = fltarr(npts) & a3 = fltarr(npts)
a1 = a & a2 = a*cos(8.*(theta-3.*!pi/16.))
a3 = -a
s1 = 0.5*a*theta^2
s2 = a*(-(1.0/64.)*cos(8.*(theta-3.0*!pi/16.))+3.0*!pi/16.*theta+1./64.-9.*!pi^2/512)
s3 = a*(-0.5*theta^2+0.5*!pi*theta+1./32.-17.*!pi^2/256.)
v1 = fltarr(npts)
v1 = a*theta
v2 = fltarr(npts)
v2 = a*((1.0/8.0)*sin(8.0*(theta-3.0*pi/16.0))+3.0*pi/16.0)
v3 = fltarr(npts)
v3 = a*(0.5*pi-theta)

newpts = n1+n2+n3
v = fltarr(newpts) & v = [v1[region1],v2[region2],v3[region3]]
s = fltarr(newpts) & s = [s1[region1],s2[region2],s3[region3]]
acc = fltarr(newpts) & acc = [a1[region1],a2[region2],a3[region3]]
finalpts = 4*newpts
vnew = fltarr(finalpts) & snew = fltarr(finalpts) & anew = fltarr(finalpts)
snew = [s[0:newpts-1],max(s[0:newpts-1])-s[0:newpts-1],s[0:newpts-1],max(s[0:newpts-1])-s[0:newpts-1]]
vnew = [v[0:newpts-1],-v[0:newpts-1],v[0:newpts-1],-v[0:newpts-1]]
anew = [acc[0:newpts-1],-acc[0:newpts-1],acc[0:newpts-1],-acc[0:newpts-1]]
snew = shift(snew,-fix(0.25*finalpts))
vnew = shift(vnew,-fix(0.25*finalpts))
anew = shift(anew,-fix(0.25*finalpts))

ango = 0.0 & anghi = 1.0*pi
dang = (anghi-ango)/(4*npts-1)
ang = fltarr(4*npts)
ang = ango+dang*findgen(4*npts)

ango1 = 0.0 & anghi1 = 2.0*pi
dang1 = (anghi1-ango1)/(finalpts-1)
ang1 = fltarr(finalpts)
ang1 = ango1+dang1*findgen(finalpts)

sout = fltarr(npts)
sout = interpol(snew,ang1,ang)
aout = fltarr(npts)
aout = interpol(anew,ang1,ang)
yout = fltarr(npts)
yout = interpol(vnew,ang1,ang)

yout = pi*f*yout
youtnew = fltarr(npts)
youtnew = congrid(yout,npts)
soutnew = congrid(sout,npts)
aoutnew = fltarr(npts)
aoutnew = (!pi*f)^2*congrid(aout,npts)

;plot,aoutnew
;oplot,youtnew,linestyle=3

data = fltarr(2,npts)
t = fltarr(npts)
tlo = 0.0 & thi = 1.0/f
dt = (thi-tlo)/(npts-1)
t = tlo+dt*findgen(npts)
k1 = 5.22707e-3 & k2 = 1261.64 & vb = k2*0.5
dsad = 4.12 & dms = 2.25
e = fltarr(npts) & efrench = fltarr(npts)
td = fltarr(npts) & tfrench = fltarr(npts)
td = t+(dsad/vb)+(dms+soutnew)/(vb+youtnew)
tfrench = t + ((dsad+dms)/vb)
efrench = k1*(k2*youtnew+youtnew^2)
e = interpol(efrench,td,tfrench)
;plot,tfrench,e,psym=0;,xrange=[0.02,0.03],yrange=[-2.,2.]
;oplot,tfrench,efrench,linestyle=3

data[0,0:npts-1] = t[0:npts-1]
data[1,0:npts-1] = shift(youtnew[0:npts-1],-fix(0.25*npts))
return,data
end