;+
;Procedure: tvector_rotate
;
;Purpose: rotates a set of vectors by a set of coordinate
; transformation matrices inputs and outputs are tplot
; variables. This is designed mainly for use with
; fac_matrix_make.pro and minvar_matrix_make, but can be used
; for more general purposes
;
;Warning:
; The transformation matrices generated by
;fac_matrix_make,thm_fac_matrix_make, and minvar_matrix_make are
;defined relative to a specific coordinate system. This means that if
;you use vectors in one coordinate system to generate the rotation
;matrices they will only correctly transform data from that coordinate
;system into the functional coordinate system. For example if you use
;magnetic field data in gsm to generate Rgeo transformation matrices
;using fac_matrix_make then the vectors being provided to tvector
;rotate to be transformed by those matrices should only be in gsm coordinates.
;
;Arguments:
;
; mat_var_in: the name of the tplot variable storing input matrices
; The y component of the tplot variable's data struct should be
; an Mx3x3 array, storing a list of transformation matrices.
;
;
; vec_var_in: the name of the tplot variable storing input
;vectors. You can use globbing to rotate several tplot variables
;storing vectors with a single matrix. The y component of the tplot variable's
; data struct should be an Nx3 array.
;
;
;
; 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:
; mat_var_in should store rotation matrices. (ie the columns of any
; matrix in mat_var_in should form an orthonormal basis)
;
; If M!=N M must be >= 2 (where M,N refer to the dimensions of the
; input variables.
;
; Also If the x components of mat_var_in and vec_var_in data structs
; do not match then. The matrices will be interpolated to match the
; cadence of the vector data. Interpolation is done by turning the
; matrices into quaternions and performing spherical linear
; interpolation. As noted above this interpolation behavior will
; not be predictable if the matrices do anything other than rotate.
;
; The user should note that if the timestamps of mat_var_in are not
; monotonic(ascending or identical) or if the timestamps vec_var_in are not strictly
; monotonic(always ascending) the program will not work correctly. It is the user's
; responsibility to remove duplicate or non sequential timestamps from the data
;
; PROGRAMMERS: you might want to consider adding a check to
; verify that input matrices are actually rotation matrices
; this would just verify that M dot transpose(M) = Identity matrix
;
; SEE ALSO: mva_matrix_make.pro, fac_matrix_make.pro
;
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2009-03-30 14:17:08 -0700 (Mon, 30 Mar 2009) $
; $LastChangedRevision: 5454 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_5_00/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'
;loop over possibly many vectors
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
;load the matrices and validate dimensions
get_data,mat_var_in,data=m_d,dlimits=m_dl ;rotation matrices and input abcissa values
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
;load the vectors and validate dimensions
get_data,tn[i],data=v_d,dlimits=v_dl ;output abcissa values
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
idx = where(m_d.x ne v_d.x)
if n_elements(m_d.x) eq n_elements(v_d.x) && idx[0] eq -1 then begin
m_d_y = m_d.y
endif else begin
if n_elements(m_d.x) gt 1 then begin
idx = where((m_d.x[1:n_elements(m_d.x)-1]-m_d.x[0:n_elements(m_d.x)-2]) lt 0)
if(idx[0] ne -1) then begin
message,'matrix timestamps not monotonic please remove any non-ascending timestamps from your data'
endif
endif
if n_elements(v_d.x) gt 1 then begin
idx = where((v_d.x[1:n_elements(v_d.x)-1]-v_d.x[0:n_elements(v_d.x)-2]) le 0)
if(idx[0] ne -1) then begin
message, 'vector timestamps not monotonic please remove any non-ascending or repeated timestamps from your data'
endif
endif
;using quaternion library transform into quanternions
q_in = mtoq(m_d.y)
;interpolate quaternions
q_out = qslerp(q_in,m_d.x,v_d.x)
;turn quaternions back into matrices
m_d_y = qtom(q_out)
endelse
if(size(m_d_y,/n_dim) eq 0 && m_d_y[0] eq -1L) then begin
message,/continue,'failed to interpolate rotation matrix array'
return
endif
;rotate vectors using matrices
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
;update meta data
;and store output
;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,/continue,'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