;+
; NAME:
;       MSD_MOMENTS
;
; PURPOSE:
;       This function computes the ensemble mean and variance of individual
;       mean-squared displacement. 
;
; CALLING SEQUENCE:
;       Result = MSD_MOMENTS(t [, TimeStep = TimeStep, MicPerPix = MicPerPix,$
;                 Mydts = Mydts, N_Lag = N_Lag, Dim = Dim])
;
; INPUTS:
;       t:      The result of TRACK.PRO
;
; KEYWORD PARAMETERS:
;       TimeStep:  Time interval between successive frame/field. Default is 1.
;
;       MicPerPix: Magnification in microns per pixel. If it is single number,
;                  the magnification is assumed to be isotropic, otherwise,
;                  provide a dim-dimensional vector. Default is 1.
;
;       Mydts:     A vector of integer to calculate the output only at these
;                  lags.
;
;       N_Lag:     The number of lagtime at which the output is calculated. The
;                  lags are log-spaced points. Useless if Mydts is supplied.
;
;       Dim:       The spatial dimensionality of the data. Default is 2.
;
; OUTPUTS:
;	Result[0,*]        Lagtimes  
;	Result[1,*]        Mean displacements  
;	Result[dim+1,*]    Ensemble average of mean squared displacement  
;	Result[2*dim+1,*]  Ensemble variance of mean squared displacement
;	Result[3*dim+1,*]  Variance of Result[dim+1,*]
;	Result[4*dim+1,*]  Variance of Result[2*dim+1,*] 
;	Result[5*dim+1,*]  Transparency  
;
; EXAMPLE:
;       The array t is the output of tracking
;        m = MSD_MOMENTS(t)
;       Plot m(3,*) vs. m(0,*) gives the mean squared displacement in the x direction
;       and m(7,*)^0.5 gives a 1-sigma error bar.
;
;       Plot m(5,*) vs. m(0,*) gives the variance of individual mean squared
;       displacement in the x direction (to check heterogeneity), and m(9,*)^0.5 gives
;       the error bar.
;
; REFERENCE:
;       Savin, T. and Doyle, P.S.
;       "Statistical and sampling issues when using multiple particle tracking"
;       Phys. Rev. E, 76:021501, 2007.
;
; MODIFICATION HISTORY:
;       Written by:  Thierry Savin @ Doyle Group @ MIT - June 2006
;-

FUNCTION MSD_MOMENTS, t, TimeStep = TimeStep, MicPerPix = MicPerPix, Mydts = Mydts, N_Lag = N_Lag, Dim = Dim

IF NOT KEYWORD_SET(TimeStep) THEN TimeStep = 1d ELSE TimeStep = DOUBLE(TimeStep)
IF NOT KEYWORD_SET(MicPerPix) THEN MicPerPix = 1d ELSE MicPerPix = DOUBLE(MicPerPix)
IF NOT KEYWORD_SET(Dim) THEN Dim = 2

Dim_t = N_ELEMENTS(t[*,0])
Ind_Id = Dim_t - 1
Ind_Time = Dim_t - 2

; Characteristics of the trajectories
Id_Pos = UNIQ( t[Ind_Id,*] )    ; Subscript of the end
N_Id = N_ELEMENTS( Id_Pos )     ; Number
Id = t[Ind_Id,Id_Pos]           ; Identities
Id_Length = Id_Pos - SHIFT(Id_Pos,1) - 1 & Id_Length[0] = Id_Pos[0] ; Number of displacement at lag=1

; Create the list of Lag
Max_Lag = (N_Id GE 2) ? Id_Length[ (SORT(Id_Length))[N_Id-2] ] : MAX(Id_Length) - 1
IF KEYWORD_SET(Mydts) THEN BEGIN
    Lag = LONG(Mydts[WHERE(Mydts LE Max_Lag)])
ENDIF ELSE BEGIN
    IF NOT KEYWORD_SET(N_Lag) THEN N_Lag = 50
    IF (Max_Lag LE N_Lag) THEN Lag = LINDGEN(Max_Lag)+1 ELSE BEGIN
        N_Ind = N_Lag-1
        REPEAT BEGIN
            N_Ind++
            Index = DINDGEN(N_Ind)/(N_Ind-1)
            Lag = ROUND(Max_Lag^Index)
            Lag = Lag[UNIQ(Lag)]
        ENDREP UNTIL N_ELEMENTS(Lag) EQ N_Lag
    ENDELSE
ENDELSE
N_Lag = N_ELEMENTS(Lag)

; Calculate N_box
Time = t[Ind_Time,SORT(t[Ind_Time,*])]
Time_Pos = UNIQ(Time)
N_box = Time_Pos - SHIFT(Time_Pos,1) & N_box[0] = Time_Pos[0] + 1
N_box = TOTAL( N_box, /DOUBLE ) / N_ELEMENTS( N_Box )

; Pos, Time, Id
xy_Id = TRANSPOSE( [ t[0:Dim-1,*] * REBIN([MicPerPix],Dim,N_ELEMENTS(t[0,*]),/SAMPLE), t[ [Ind_Time,Ind_Id], * ] ] )

Result = DBLARR(5*Dim+2,N_Lag,/NOZERO) + !VALUES.D_NAN
Result[0,*] = Lag * TimeStep


FOR j=0, N_Lag - 1 DO BEGIN

   ; List of Displacements
    dx = xy_Id - SHIFT( xy_Id, Lag[j], 0 )
    dx = dx[ WHERE( (dx[*,Dim] EQ Lag[j]) AND (dx[*,Dim+1] EQ 0), N_dx), 0:Dim-1]

    ; Individual particle positions
    W_Good = WHERE(Id_Length GE Lag[j], N_Id)
    Norm = TOTAL(Id_Length[W_Good])
    Id_Length = Id_Length[W_Good]
    N_dx_I = Id_Length + 1 - Lag[j]
    Pos_Low = (N_Id EQ 1) ? 0 : [0, TOTAL(N_dx_I[0:N_Id-2], /CUMULATIVE, /PRESERVE_TYPE)]

    ; Weights
    Weight = DBLARR(N_dx, /NOZERO)
    FOR i = 0, N_Id - 1 DO Weight[Pos_Low[i]] = DBLARR(N_dx_I[i]) + Id_Length[i] / DOUBLE( N_dx_I[i] )
    Weight = REBIN(TEMPORARY(Weight), N_dx, Dim) / Norm

    dxb = TOTAL(dx,1) / N_dx
    dx = dx -  ( dx[*,0]*0 + 1 ) # dxb 
    dx2b = TOTAL( Weight * dx^2, 1 )
    dx4b = TOTAL( Weight * dx^4, 1 )
    dx6b = TOTAL( Weight * dx^6, 1 )
    dx8b = TOTAL( Weight * dx^8, 1 )
    
    N_Seg = Lag[j]

    ; First the one for which N_dx_I ge N_Seg
    W_Good = WHERE(N_dx_I GE N_Seg, N_Good, COMPLEMENT = W_Bad, NCOMPLEMENT = N_Bad)
    IF N_Good GE 1 THEN BEGIN
        N_dx_I_Good = N_dx_I[W_Good] - (N_dx_I[W_Good] MOD N_Seg)
        Pos_Good = [0, TOTAL(N_dx_I_Good, /CUMULATIVE, /PRESERVE_TYPE)]
        Good = LONARR(Pos_Good[N_Good], /NOZERO)
        Weight = Good * 0D
        FOR i = 0, N_Good - 1 DO BEGIN
            Good[Pos_Good[i]] = Pos_Low[ W_Good[i] ] + LINDGEN( N_dx_I_Good[i] )
            Weight[Pos_Good[i]] = DBLARR(N_dx_I_Good[i])+( Id_Length[ W_Good[i] ] + 0 )/ DOUBLE( N_dx_I_Good[i] / N_Seg )
        ENDFOR
        n = Pos_Good[N_Good] / N_seg
        Weight = REBIN(TEMPORARY(Weight), n, Dim)
        dx_2 = REBIN(dx[Good,*]^2, n, Dim)
        dx_4 = REBIN(dx[Good,*]^4, n, Dim)
        dx_6 = REBIN(dx[Good,*]^6, n, Dim)
    ENDIF

    ; Then the rest of trajectories that bring only one block displacement
    IF N_Bad GE 1 THEN BEGIN
        FOR i = 0, N_Bad - 1 DO BEGIN
            dx_I = dx[Pos_Low[W_Bad[i]]:Pos_Low[W_Bad[i]]+N_dx_I[W_Bad[i]]-1]
            Weight =( (N_Good EQ 0) && (i EQ 0) ) ? REBIN([Id_Length[W_Bad[i]]], 1, Dim) : [Weight, REBIN([Id_Length[W_Bad[i]]], 1, Dim)]
            dx_2 = ( (N_Good EQ 0) && (i EQ 0) ) ? REBIN(dx_I^2, 1, Dim) : [dx_2, REBIN(dx_I^2, 1, Dim)]
            dx_4 = ( (N_Good EQ 0) && (i EQ 0) ) ? REBIN(dx_I^4, 1, Dim) : [dx_4, REBIN(dx_I^4, 1, Dim)]
            dx_6 = ( (N_Good EQ 0) && (i EQ 0) ) ? REBIN(dx_I^6, 1, Dim) : [dx_6, REBIN(dx_I^6, 1, Dim)]
        ENDFOR
    ENDIF

    n = (N_Good EQ 0) ? N_Bad : ( n + N_Bad )
    
    Weight = TEMPORARY(Weight) / Norm
    
    dx_22b = TOTAL(Weight * dx_2^2,1)
    dx_42b = TOTAL(Weight * dx_4^2,1)
    dx_2dx_4b = TOTAL(Weight * dx_2*dx_4,1)

    Result[0*Dim+1,j] = dxb
    Result[1*Dim+1,j] = dx2b
    Result[2*Dim+1,j] = dx4b/3-dx2b^2+(dx_22b-dx2b^2) * TOTAL(Weight^2)/ TOTAL(Weight*(1-Weight))
    Result[3*Dim+1,j] = (dx_22b-dx2b^2)* TOTAL(Weight^2)/ TOTAL(Weight*(1-Weight)) + (dx4b/3-dx2b^2+(dx_22b-dx2b^2) * TOTAL(Weight^2)/ TOTAL(Weight*(1-Weight)))/N_box
    Result[4*Dim+1,j] = ((dx_42b-dx4b^2)/9+4*dx2b^2*(dx_22b-dx2b^2)+4*dx2b*(dx2b*dx4b-dx_2dx_4b)/3)* TOTAL(Weight^2)/ TOTAL(Weight) + (dx8b/105 - 4*dx2b*dx6b/15 + 2*dx2b^2*dx4b - 3*dx2b^4 - (dx4b/3-dx2b^2)^2)/N_box
    Result[5*Dim+1,j] = TOTAL(Id_Length + 1) / (N_box * ( MAX(t[Ind_Time,*])-MIN(t[Ind_Time,*]) + 1 ))

ENDFOR

RETURN, Result

END