;+
;Procedure: tvector_rotate
;
;Purpose:  rotates a set of vectors by a set of coordinate
;          transformation matrices inputs and outputs are tplot variables
;Arguments:
;
; mat_var_in: the name of the tplot variable storing input matrices
;
; vec_var_in: the name of the tplot variable storing input vectors(can
; use normal globbing)
;
;
; newname(optional): the name of the output variable, defaults to 
;                    vec_var_in + '_rot'
;                    If you use type globbing in the vector variable
;                    This option will be disabled
;
; error(optional): named variable in which to return the error state
; of the computation.  1 = success 0 = failure
;
;
;CALLING SEQUENCE:
;
; tvector_rotate,'matrix_var','vector_var',newname = 'out_var',error=error
;
; tvector_rotate,'matrix_var','vector_var'
;
;Notes: 
;--->I'm a little sparse on comments here, but I think the code
;should be self documenting
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2007-08-26 20:06:03 -0700 (Sun, 26 Aug 2007) $
; $LastChangedRevision: 1493 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_3_01/cotrans/special/tvector_rotate.pro $
;-


;helper function
;calculates the norm of a bunch of vectors simultaneously
function ctv_norm_vec_rot,v

COMPILE_OPT HIDDEN

if not keyword_set(v) then return, -1L

if size(v,/n_dim) ne 2 then return,-1L

return, sqrt(total(v^2,2))

end

;helper function
;normalizes a bunch of vectors simultaneously
function ctv_normalize_vec_rot,v

COMPILE_OPT HIDDEN

if not keyword_set(v) then return, -1L

if size(v,/n_dim) ne 2 then return,-1L

n_a = ctv_norm_vec_rot(v)

if(size(n_a,/n_dim) eq 0 && n_a eq -1L) then return,-1L

v_s = size(v,/dimension)

;calculation is pretty straight forward
;we turn n_a into an N x D so computation can be done element by element
n_b = rebin(n_a,v_s[0],v_s[1])

return,v/n_b

end

;helper function
;vectorized fx to multiply n matrices by n vectors
function ctv_mx_vec_rot,m,x

COMPILE_OPT HIDDEN

;input checks
if(not keyword_set(m)) then return, -1L

if(not keyword_set(x)) then return, -1L

m_s = size(m,/dimension)

x_s = size(x,/dimension)

;make sure number of dimensions in input arrays is correct
if(n_elements(m_s) ne 3) then return, -1L

if(n_elements(x_s) ne 2) then return, -1L

;make sure dimensions match
if(not array_equal(x_s,[m_s[0],m_s[1]])) then return,-1L

if(not array_equal(m_s,[x_s[0],x_s[1],x_s[1]])) then return,-1L

;calculation is pretty straight forward
;we turn x into an N x 3 x 3 so computation can be done element by element
y_t = rebin(x,x_s[0],x_s[1],x_s[1])

;custom multiplication requires rebin to stack vector across rows,
;not columns
y_t = transpose(y_t, [0, 2, 1])

;9 multiplications and 3 additions per matrix
y = total(y_t*m,3)

return, y

end

pro tvector_rotate,mat_var_in,vec_var_in,newname = newname,error=error

if keyword_set(error) then error = 0

if not keyword_set(mat_var_in) then begin
    message,/continue,'FX requires the matrix variable to be set'
    return
endif

if tnames(mat_var_in) eq '' then begin
    message,/continue,'FX requires the matrix variable to be set'
    return
endif

if not keyword_set(vec_var_in) then begin
    message,/continue,'FX requires vector variable to be set'
    return
endif

tn = tnames(vec_var_in)

if size(tn,/n_dim) eq 0 && tn eq '' then begin
    message,/continue,'FX requires vector variable to be set'
    return
endif

if not keyword_set(newname) then $
  newname = [vec_var_in + '_rot'] $
else $
  newname = [newname]


if n_elements(tn) gt 1 then newname = tn + '_rot'

for i = 0,n_elements(tn) -1L do begin

    tinterpol_mxn,mat_var_in,tn[i],newname = mat_var_in+'_interp_temp',error=error_state

    if error_state eq 0 then begin
        message,/continue,'TVECTOR_ROTATE failed to interpolate rotation matrix arrat'

        return

    endif

    get_data,mat_var_in+'_interp_temp',data = m_d,dlimits=m_dl

    del_data,mat_var_in+'_interp_temp'

    m_d_s = size(m_d.y,/dimensions)

    if(n_elements(m_d_s) ne 3 || m_d_s[1] ne 3 || m_d_s[2] ne 3) then begin

        message,/continue,'FX requires matrix data to contain an Nx3x3 array'
        return
    endif

    get_data,tn[i],data=v_d,dlimits=v_dl

    v_d_s = size(v_d.y,/dimensions)

    if(n_elements(v_d_s) ne 2 || v_d_s[1] ne 3) then begin
        message,/continue,'FX requires vector data to contain an Nx3 Array'
        return
    endif

    ;renormalize matrix columns
    xt = ctv_normalize_vec_rot(m_d.y[*,*,0])

    if(size(xt,/n_dim) eq 0 && xt eq -1L) then begin
        message,/continue,'FX failed to renormalize matrix x column'
        return
    endif

    yt = ctv_normalize_vec_rot(m_d.y[*,*,1])

    if(size(yt,/n_dim) eq 0 && yt eq -1L) then begin
        message,/continue,'FX failed to renormalize matrix y column'
        return
    endif

    zt = ctv_normalize_vec_rot(m_d.y[*,*,2])

    if(size(zt,/n_dim) eq 0 && zt eq -1L) then begin
        message,/continue,'FX failed to renormalize matrix z column'
        return
    endif

    m_d.y[*,*,0] = xt
    m_d.y[*,*,1] = yt
    m_d.y[*,*,2] = zt

    v_t = ctv_mx_vec_rot(m_d.y,v_d.y)
   
    if(size(v_t,/n_dim) eq 0 && v_t eq -1L) then begin
        message,/continue,'FX failed to perform mx operation'
        return
    endif

    ;if the vector dlimits struct is unavailable create one
    if(size(v_dl,/n_dim) eq 0 && v_dl eq 0) then begin
        message,/continue,'dlimits for input vector should be set'
        message,/continue,'Creating default dlimits structure'

        v_data_att = {coord_sys:'undefined',units:''}
        v_dl = {data_att:v_data_att,labflag:0,labels:make_array(v_d_s[1],/string,value='')}

    endif

    ;if the matrix dlimits struct is unavailable create one
    if(size(m_dl,/n_dim) eq 0 && m_dl eq 0) then begin
        message,/contiune,'m_dl coord_sys unset, unable to make coordinate system determination'
        
        m_data_att = {coord_sys:'undefined',units:''}
        m_dl = {data_att:m_data_att,labflag:0,labels:make_array(v_d_s[1],/string,value='')}
    endif


    ;make sure it has data_att
    str_element,v_dl,'data_att',SUCCESS=s

    if(s eq 0) then begin

        v_data_att = {coord_sys:'undefined',units:''}
        
        str_element,/add,v_dl,'data_att',v_data_att

    endif

    ;make sure it has ytitle

    str_element,v_dl,'ytitle',SUCCESS=s

    if(s eq 0) then begin

        v_ytitle = ''

        str_element,/add,v_dl,'ytitle',v_ytitle

    endif

    ;make sure it has appropriate data_att elements
    str_element,v_dl.data_att,'coord_sys',SUCCESS=s

    if(s eq 0) then begin

        v_data_att = v_dl.data_att 

        str_element,/add,v_data_att,'coord_sys','undefined'

        str_element,/add,v_dl,'data_att',v_data_att

        
    endif

    str_element,v_dl.data_att,'units',SUCCESS=s

    if(s eq 0) then begin

        v_data_att = v_dl.data_att 

        str_element,/add,v_data_att,'units',''

        str_element,/add,v_dl,'data_att',v_data_att

    endif

    v_dl.data_att.coord_sys = m_dl.data_att.coord_sys
    v_dl.data_att.units = v_dl.data_att.units+'_rot'
    v_dl.ytitle = newname[i] + '!C' + v_dl.data_att.units 

    ;no longer take labels from transformation matrix
    ;if m_dl.labflag eq 1 then v_dl.labels = m_dl.labels

    str_element,v_d,'v',SUCCESS=s

    if(s) then $
      v_d = {x:v_d.x,y:v_t,v:v_d.v} $
    else $
      v_d = {x:v_d.x,y:v_t} 

    store_data,newname[i],data=v_d,dlimits = v_dl

endfor

error = 1

return

end