FUNCTION rmrot,data, dim, rotmat = Rret, center = center
 
;+
;    
; NAME:
;   RMROT
; 
; PURPOSE:
;    Remove average rotational motion of a sample, similar to
;    rm_motion, but without the 'smoo' variable.  
;        
;    
; CALLING SEQUENCE:
;   trn = rmrot(data,dim,rotmat=RR,center=center)
;   
;    
;    
; INPUTS:
;   data: this is a track file (no pretrack yet) either 2d or 3d
;   dim:  this is the no. of dimensions you want to 'unrotate' or
;         calculate a rotation matrix for
;    
; KEYWORDS:
;   center: The procedure requires that the data at each timestamp
;           be centered about the origin.  If the data you enter
;           isn't centered, specify /center when calling.
;   rotmat: specify rotmat=RR to return the rotation matrices.
;    
; OUTPUTS:
;    trn:  a track file one-to-one with your original file, 
;          but whose coordinates have been transformed in 
;          such a way to remove cumulative bulk rotations.
;    Rret:   A set of rotation matrices that describes the rotations
;          between each timestamp.  The dimensions of RR
;          will be 3 X 3 X (number of timestamps-1)
;
;    
; MODIFICATION HISTORY:
;     Written by Gary L. Hunter - 05/21/09
;     Adapted and modified from:
;     J Challis, Journal of Biomechanics, Vol. 28, No. 6, June 1997,
;     pp.733-737
;     For more info, see 
;     Tracking rotational diffusion of colloidal clusters
;     Gary L. Hunter, Kazem V. Edmond, Mark T. Elsesser, and Eric R. Weeks
;     Optics Express, Vol. 19, Issue 18, pp. 17189-17202 (2011) 
;    
;-
    
temp = data
IF keyword_set(center) THEN BEGIN
temp = centertrack(temp, dim)
endif

tr = double(temp)

nc = n_elements(tr[*,0])

tlist = unique(tr[nc-2,sort(tr[nc-2,*])])
nts = n_elements(tlist)

Ri = identity(dim)              ; define identity matrix

Rret = dblarr(dim,dim,nts-1)


trn = tr

mxt = max(tlist)

FOR t = 0L,nts-2 DO begin

IF t MOD 100 EQ 0 THEN begin
 print,'doing time '+n2s(t)+' of '+n2s(mxt)
endif

   t1 = tlist[t]
   t2 = tlist[t+1]

   w1 = where(tr[nc-2,*] EQ t1, co1)
   w2 = where(tr[nc-2,*] EQ t2, co2)

   tr1 = tr[*,w1]
   tr2 = tr[*,w2]

   id1 = tr1[nc-1,*]
   id2 = tr2[nc-1,*]

   good = setintersection(id1,id2)
   n = n_elements(good)
   C = dblarr(dim,dim)

   FOR i = 0L,n-1 DO BEGIN 

      wid1 = where(tr1[nc-1,*] EQ good[i]) 
      wid2 = where(tr2[nc-1,*] EQ good[i]) 

      x1 = transpose(tr1[0:dim-1,wid1]) 
      x2 = transpose(tr2[0:dim-1,wid2]) 

      C = C + x2##transpose(x1) 

   ENDFOR

   C = C/n

   SVDC,C,W,U,V

   temp = identity(dim)
   temp[dim-1,dim-1] = determ(u##transpose(v))

   R = U##temp##transpose(V)

   Rret[*,*,t] = R              ;now, we can use the rotation info. 07/23/09 GLH
   
   Ri = Ri##invert(R)           ;account for past rotations

   FOR j = 0L,co2 - 1 DO BEGIN

      trn[0:dim-1,w2[j]] = transpose(Ri##transpose(trn[0:dim-1,w2[j]])) 


   ENDFOR 

ENDFOR

return,trn

END