;+
;Procedure: fac_matrix_make
;
;Purpose:  generates a field aligned coordinate transformation matrix
;from an input B vector array(and sometimes a position vector array)
;then stores it in a tplot variable
;Arguments:
;   mag_var_name=the name of the tplot variable storing the magnetic field
;   vectors to be used in transformation matrix generation
;   pos_var_name(optional)=the name of the tplot variable storing the position
;   vectors to be used in transformation matrix generation
;   newname(optional)=the name of the tplot variable in which to store
;   the output
;   error(optional) = named variable that holds the error state of the 
;   computation 1 = success 0 = failure
;   other_dim(optional) = the second axis for the field aligned
;   coordinate system.
;   ************For all transformations Z = B************
;
;   valid second coord(other_dim) options:
;
;         'Xgse', translates from gse or gsm into FAC
;                    Definition(works on GSE or GSM):
;                    X Axis = on plane defined by Xgse - Z
;                    Second coordinate definition: Y = Z x X_gse
;                    Third coordinate, X completes orthogonal RHS
;                    (right hand system) triad: XYZ
;         'Rgeo',translate from geo into FAC using radial posision vector
;                    R_pos is radial position vector, radialy inwards.
;                    Second coordinate definition: Y = Z x Rpos
;                    Third coordinate, X completes orthogonal RHS XYZ.
;         'Ygeo', translate from geo into FAC using azimuthal position vector
;                    Y_pos is the azimuthal position vector, positive Westward
;                    Second coordinate definition: Y = Z x Y_pos
;                    Third coordinate, X completes orthogonal RHS XYZ
;         'Ysm', translate from geo into FAC using azimuthal sm position
;                 Y Axis = on plane defined by Phi_sm- Z
;                 Second coordinate definition: X = Phi_sm x Z
;                 Third completes orthogonal RHS XYZ
; Example:
;  fac_matrix_make,'tha_fgs',other_dim='Xgse',pos_var_name='tha_pos',out_var_name='tha_fgs_fac_mat'
;
;Notes:--> I'm a little sparse on comments here, but I think the code
;should be self documenting
;--> Should filter NaNs to supress floating point errors
;--> Contains coordinate transformation specific code, if new
;    coordinate systems are added, this code should be updated
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2007-09-07 16:45:43 -0700 (Fri, 07 Sep 2007) $
; $LastChangedRevision: 1558 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_3_01/cotrans/special/fac/fac_matrix_make.pro $
;-

pro fac_matrix_make,mag_var_name,other_dim=other_dim, pos_var_name=pos_var_name,newname = newname,error=error

;the valid_coords array is positional...so the order of coordnames
;should match their processing order in the code below

valid_coords = ['xgse','rgeo','ygeo','ysm']

error = 0

if not keyword_set(mag_var_name) then begin
    message,/continue,' fx requires mag_var_name to be set'
    return
endif

if tnames(mag_var_name) eq '' then begin
    message,/continue,' fx requires mag_var_name to be set'
    return
endif

if not keyword_set(other_dim) then begin
    other_dim='xgse'
endif

if not keyword_set(newname) then newname = mag_var_name + '_fac_mat'

other_dim_l = strlowcase(other_dim)

co_idx = where(strmatch(valid_coords, other_dim_l) ne 0)

if(co_idx[0] eq -1L) then begin
    message,/continue,' fx was passed an illegal output coordinate'+other_dim
    message,/continue,' fx requires output other_dim to be: ' + strjoin(valid_coords,',')
    return
end

get_data,mag_var_name,data=d,limits=l,dlimits=dl

;temp variable, needs to be here, cause if there are no gaps tdegap
;fails, otherwise tdegap will overwrite
store_data,mag_var_name+'_ctv_temp',data=d,limits=l,dlimits=dl

;remove any gaps in the data
tdegap,mag_var_name,newname=mag_var_name+'_ctv_temp'

tnormalize,mag_var_name+'_ctv_temp',newname='fac_mat_z_temp',error=error_state

if error_state eq 0 then begin
    message,/continue,' fx failed to normalize magnetic field vector'
    return
endif

del_data,mag_var_name+'_ctv_temp'

get_data,'fac_mat_z_temp',data=d,dlimits=dl

d_s = size(d.y,/dimensions)

if(n_elements(d_s) ne 2 || d_s[1] ne 3) then begin
       message,/continue,' fx requires mag_var_name data to be an Nx3 Array'
       return
endif

;rhs coordinates constructed on plane Xgse - B
if(other_dim_l eq valid_coords[0]) then begin

    if(dl.data_att.coord_sys ne 'gse' && $
       dl.data_att.coord_sys ne 'gsm') then begin
        message,/continue,' fx requires mag_var_name to be in GSE or GSM to generate a Xgse field aligned coordinate matrix'
        return
    endif

    get_data,'fac_mat_z_temp',data=d,dlimits=dl

    ;constructs an array of unit vectors
    ;for use in generation of series of
    ;unit vector cross products
    x_axis = transpose(rebin([1D,0D,0D],3,n_elements(d.x)))

    store_data,'fac_mat_x_temp',data={x:d.x,y:x_axis,v:d.v},dlimits=dl

    tcrossp,'fac_mat_z_temp','fac_mat_x_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
         message,/continue,' FX  failed to create zx cross product in ' + valid_coords[0]
        return
    endif

    tnormalize,'fac_mat_y_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
         message,/continue,' FX  failed to create zx cross product in ' + valid_coords[0]
        return
    endif


    tcrossp,'fac_mat_y_temp','fac_mat_z_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to create yz cross product in ' + valid_coords[0]
        return
    endif

;the code above assumes the existence of coord_sys, but it probably shouldn't
;    str_element,dl,'data_att',success=s2

;    s3 = 0

 ;   if s2 then $
 ;     str_element,dl.data_att,'coord_sys',success=s3

    
    str_element,dl,'labels',success=s1

    if s1 eq 1 then $
      dl.labels = ['Y x B','B x X_'+dl.data_att.coord_sys,'B']

endif else $

;coordinates constructed on rhs coordinate system defined by R_pos - B
if(other_dim_l eq valid_coords[1]) then begin

    if not keyword_set(pos_var_name) then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    if tnames(pos_var_name) eq '' then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    ;remove any gaps in the data
    ;tdegap,pos_var_name,newname=pos_var_name+'_ctv_temp'
    ;get_data,pos_var_name+'_ctv_temp',data=d_2,dlimits=dl_2

    get_data,pos_var_name,data=d_2,dlimits=dl_2

    ;del_data,pos_var_name+'_ctv_temp'

    d_2_s = size(d_2.y,/dimensions)

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

    if dl_2.data_att.coord_sys ne 'gei' then begin
        message,/continue,' FX  requires position data to be in gei coordinates to generate Rpos transformation'
        return
    endif

    tinterpol_mxn,pos_var_name,'fac_mat_z_temp',newname='fac_mat_pos_temp',error=error_state

    if error_state eq 0 then begin
        message,/continue,' FX  failed to interpolate position data onto mag field data'
        return
    endif

    tnormalize,'fac_mat_pos_temp',newname='fac_mat_pos_temp',error=error_state

    if error_state eq 0 then begin
        message,/continue,' FX  failed to normalize interpolated position data'
        return
    endif

    ; --vassilis begin insertion--
    case dl.data_att.coord_sys of
     'gse' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gei2gse
             end
     'gsm' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gei2gse
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gse2gsm
             end
      else : begin
               message,/continue,' FX  failed to rotate data onto mag field system'
               return
             end
    endcase
    ; --vassilis end insertion--

    tcrossp,'fac_mat_pos_temp','fac_mat_z_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to generate R_pos x B vector'
        return
    endif

    tcrossp,'fac_mat_y_temp','fac_mat_z_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to generate Y x Z vector'
        return
    endif

    del_data,'fac_mat_pos_temp'

    
    str_element,dl,'labels',success=s1

    if s1 eq 1 then $
      dl.labels = ['Y x B','B x Rpos','B']



endif else $

;begin ypos transformation
if(other_dim_l eq valid_coords[2]) then begin

    if not keyword_set(pos_var_name) then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    if tnames(pos_var_name) eq '' then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    ;remove any gaps in the data
    ;tdegap,pos_var_name,newname=pos_var_name+'_ctv_temp'

    get_data,pos_var_name,data=d_2,dlimits=dl_2

    ;del_data,pos_var_name+'_ctv_temp'

    d_2_s = size(d_2.y,/dimensions)

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

    if dl_2.data_att.coord_sys ne 'gei' then begin
        message,/continue,' FX  requires position data to be in gei coordinates to generate ypos transformation'
        return
    endif
    
    tinterpol_mxn,pos_var_name,'fac_mat_z_temp',newname='fac_mat_pos_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to interpolate position data onto mag field data'
        return
    endif

    ;obtain spherical coordinate unit vector phi in GEI
    ;first get theta,phi
    ;for some reason it truncates the names on this call so the names
    ;of the output variables are stored explicitly to make sure we can
    ;keep track of them
    xyz_to_polar,'fac_mat_pos_temp', magnitude = mag_name, theta = th_name,phi = phi_name
    get_data,phi_name,phi_t,phi_d
    get_data,th_name,theta_t,theta_d
    del_data, mag_name
    del_data, th_name
    del_data, phi_name

    ;allocate temporary storage variable
    y_2 = dindgen(n_elements(theta_t),3)

    ; next get unit vector phi coordinates in GEI system (overwrite y_2)
    y_2(*,0)=cos(theta_d*!PI/180.)*cos(phi_d*!PI/180.)
    y_2(*,1)=sin(theta_d*!PI/180.)*cos(phi_d*!PI/180.)
    y_2(*,2)=-sin(phi_d*!PI/180.)

    ; transform into mag field coordinate system
    store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v}
    case dl.data_att.coord_sys of
     'gse' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gei2gse
             end
     'gsm' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gei2gse
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gse2gsm
             end
      else : begin
               message,/continue,' FX  failed to rotate data onto mag field system'
               return
             end
    endcase

    tcrossp,'fac_mat_pos_temp','fac_mat_z_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate pos by z cross product'
        return
    endif

    tnormalize,'fac_mat_x_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to normalize x temp'
        return
    endif

    tcrossp,'fac_mat_z_temp','fac_mat_x_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate z by x cross product'
        return
    endif

    del_data,'fac_mat_pos_temp'

    
    str_element,dl,'labels',success=s1

    if s1 eq 1 then $
      dl.labels=['Y x B','B x Y_pos','B']

endif

;begin ypos transformation
if(other_dim_l eq valid_coords[3]) then begin

    if not keyword_set(pos_var_name) then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    if tnames(pos_var_name) eq '' then begin
        message,/continue,' FX  requires pos_var_name to be set for Rpos coordinate transformation'
        return
    endif

    ;remove any gaps in the data
    ;tdegap,pos_var_name,newname=pos_var_name+'_ctv_temp'

    get_data,pos_var_name,data=d_2,dlimits=dl_2

    ;del_data,pos_var_name+'_ctv_temp'

    d_2_s = size(d_2.y,/dimensions)

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

    if dl_2.data_att.coord_sys ne 'gei' then begin
        message,/continue,' FX  requires position data to be in gei coordinates to generate ypos transformation'
        return
    endif

    cotrans,pos_var_name,pos_var_name+'_t',/gei2gse
    cotrans,pos_var_name+'_t',pos_var_name+'_t',/gse2gsm
    cotrans,pos_var_name+'_t',pos_var_name+'_t',/gsm2sm
    
    tinterpol_mxn,pos_var_name+'_t','fac_mat_z_temp',newname='fac_mat_pos_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to interpolate position data onto mag field data'
        return
    endif

    ;obtain spherical coordinate unit vector phi in GEI
    ;first get theta,phi
    ;for some reason it truncates the names on this call so the names
    ;of the output variables are stored explicitly to make sure we can
    ;keep track of them
    xyz_to_polar,'fac_mat_pos_temp', magnitude = mag_name, theta = th_name,phi = phi_name
    get_data,phi_name,phi_t,phi_d
    get_data,th_name,theta_t,theta_d
    del_data, mag_name
    del_data, th_name
    del_data, phi_name

    ;allocate temporary storage variable
    y_2 = dindgen(n_elements(theta_t),3)

    ; next get unit vector phi coordinates in GEI system (overwrite y_2)
    y_2(*,0)=cos(theta_d*!PI/180.)*cos(phi_d*!PI/180.)
    y_2(*,1)=sin(theta_d*!PI/180.)*cos(phi_d*!PI/180.)
    y_2(*,2)=-sin(phi_d*!PI/180.)

    ;transform into mag field coordinate system
    ;how/should I perform this transform?
    store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v}
    case dl.data_att.coord_sys of
     'gse' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/sm2gsm
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/gsm2gse
             end
     'gsm' : begin
               cotrans,'fac_mat_pos_temp','fac_mat_pos_temp',/sm2gsm
             end
      else : begin
               message,/continue,' FX  failed to rotate data onto mag field system'
               return
             end
    endcase

    tcrossp,'fac_mat_pos_temp','fac_mat_z_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate pos by z cross product'
        return
    endif

    tnormalize,'fac_mat_x_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to normalize x temp'
        return
    endif

    tcrossp,'fac_mat_z_temp','fac_mat_x_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate z by x cross product'
        return
    endif

    del_data,'fac_mat_pos_temp'
    
    str_element,dl,'labels',success=s1

    ;what are the labels?
    if s1 eq 1 then $
      dl.labels=['Y x B','B x Ysm','B']

endif


get_data,'fac_mat_x_temp',data=d_x

del_data,'fac_mat_x_temp'

get_data,'fac_mat_y_temp',data=d_y

del_data,'fac_mat_y_temp'

get_data,'fac_mat_z_temp',data=d_z

del_data,'fac_mat_z_temp'

;generate output variable and store it
out = dindgen(d_s[0],3,3)

out[*,0,*] = d_x.y
out[*,1,*] = d_y.y
out[*,2,*] = d_z.y

;out = transpose(out, [0, 2, 1])

d = {x:d.x,y:out,v:d.v}

dl.data_att.coord_sys=valid_coords[0]

;tvector rotate uses this flag to determine if it should take axes
;labels from the transformation matrix

str_element,dl,'labflag',success=s

if s eq 1 then $
  dl.labflag = 1 $
else $
  str_element,dl,'labflag',value=1,/add_replace


store_data,newname,data=d,dlimits=dl

error = 1

return

end