FUNCTION rotmsd,rotmat,rotvec = rotvec, maxtime = maxtime, gettheta = gettheta, getvec = getvec, quiet = quiet, noplot = noplot, mydts = mydts

;+
; PURPOSE:
;   Calculate rotational msds (MSADS) of a vector given a set of
;   rotation matrices.
;
; REQUIRED INPUTS:
;   'rotmat' is the rotation matrix produced by rmrot.pro
;
; OPTIONAL INPUTS:
;   'rotvec' is the initial orientation vector you wish to rotate
;   according to rotmat.  Default = [1,0,0]
;   
;   'maxtime' is the maximum dt to use in msad calculations.
;
;   See msd.pro for information on keywords: quiet, noplot, mydts.
;
; OPTIONAL OUTPUTS:
;   
;   'gettheta' optional keyword to output the displacements
;   in rotation space.  See REF [1] for more info.
;
;   'getvec' optional keyword to output the orientation
;   vector as a function of time
;
;   See msd.pro for information on keywords: quiet, noplot, mydts.
;
; REFERENCES:
;  [1] Gary L. Hunter, Kazem V. Edmond, Mark T. Elsesser, and Eric
;  R. Weeks Optics Express, Vol. 19, Issue 18, pp. 17189-17202 (2011)
;
; HISTORY:
;  Written by Gary L. Hunter, 05/2009
;-

dim = n_elements(rotmat[0,*])

nts =  float(n_elements(rotmat[0,0,*]))

IF ~keyword_set(maxtime) THEN maxtime = nts


;see if 'rotvec' is specified
IF ~keyword_set(rotvec) THEN begin

x0 = double([1,0,0])

ENDIF ELSE BEGIN

x0 = double(rotvec)

ENDELSE 
;
;make sure orientation vector is properly normalized
x0 = x0/norm(x0)

;set up orientation vector array
x = dblarr(4,nts+1)
x[0:2,0] = x0
x[3,*] = lindgen(nts+1)

;apply rotations
FOR i = 1L,nts DO x[0:2,i] =  transpose(rotmat[*,*,i-1] ## transpose(x[0:2,i-1]))
print,'rotations completed'

;theta is to be an array almost exactly like a typical 3d track file.
;only columns [0:2] and [7:8] have meaning
;columns [0:2] are positions in rotation space corresponding to
;cumulative rotations about axes [x,y,z] respectively
theta = fltarr(9, nts+1)

;populate theta array
FOR i = 1L, nts DO BEGIN 
temp = crosspn(x[0:2, i], x[0:2, i-1])
theta[0:2, i] = asin(norm(temp))*temp/norm(temp)+theta[0:2, i -1]
endfor

;create timestamp column
theta[7, *] = lindgen(nts+1)

IF keyword_set(gettheta) THEN BEGIN
res = theta
ENDIF ELSE begin
 IF keyword_set(getvec) THEN BEGIN
 res = x
 ENDIF ELSE BEGIN
 res = msd(theta, dim=3, maxtime=maxtime, /keepdrift, quiet=quiet, noplot=noplot,mydts=mydts)
 ENDELSE
ENDELSE

return, res


end