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