;+
;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.  This
;   variable must be in gse or gsm coordinates when passed to this procedure.
;   pos_var_name(optional)=the name of the tplot variable storing the position
;   vectors to be used in transformation matrix generation.  This
;   variable must be in gei coordinates when passed to this procedure.
;   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.
;   /DEGAP: Set to call TDEGAP to remove any gaps from the data. See TDEGAP for
;           for other options that can be invoked using the _extra keyword.
;           E.g. thm_fac_matrix_make, 'tha_fgs', other_dim='xgse', /degap, dt=3
;
;   ************For all transformations Z = B************
;
;   Warning about coordinate systems:
;   B field tplot variable must be in gse or gsm coordinates,
;   depending on what transformation has been selected.
;   Position tplot variable must be in gei coordinates. Gei is the default coordinate
;   system of thm_load_state.
;
;   Warning:  The resulting transformation matrices will only correctly
;   transform data from the coordinate system of the input variable to
;   the field aligned coordinate system.  So if mag_var_name is in gse
;   coordinates then you should only use the output matrices to transform
;   other data in gse coordinates.
;
;   valid second coord(other_dim) options:
;
;         'Xgse', (DEFAULT) 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
;                    Note: X_gse is a unit vector pointing in direction from
;                          earth to the sun
;         'Rgeo',translate from geo into FAC using radial position vector
;                    Rgeo is radial position vector, positive radialy outwards.
;                    Second coordinate definition: Y = Z x Rgeo (eastward)
;                    Third coordinate, X completes orthogonal RHS XYZ.
;         'mRgeo',translate into FAC using radial position vector
;                    mRgeo is radial position vector, positive radially inwards.
;                    Second coordinate definition: Y = Z x mRgeo (westward)
;                    Third coordinate, X completes orthogonal RHS XYZ.
;         'Phigeo', translate into FAC using azimuthal position vector
;                    Phigeo is the azimuthal geo position vector, positive Eastward
;                    First coordinate definition: X = Phigeo x Z (positive outwards)
;                    Second coordinate, Y ~ Phigeo (eastward) completes orthogonal RHS XYZ
;         'mPhigeo', translate into FAC using azimuthal position vector
;                    mPhigeo is minus the azimuthal geo position vector; positive Westward
;                    First coordinate definition: X = mPhigeo x Z (positive inwards)
;                    Second coordinate, Y ~ mPhigeo (Westward) completes orthogonal RHS XYZ
;         'Phism', translate into FAC using azimuthal Solar Magnetospheric vector.
;                 Phism is "phi" vector of satellite position in SM coordinates.
;                 Y Axis = on plane defined by Phism-Z, normal to Z
;                 Second coordinate definition: X = Phism x Z
;                 Third completes orthogonal RHS XYZ
;         'mPhism', translate into FAC using azimuthal Solar Magnetospheric vector.
;                 mPhism is minus "phi" vector of satellite position in SM coordinates.
;                 Y Axis = on plane defined by Phism-Z, normal to Z
;                 Second coordinate definition: X = mPhism x Z
;                 Third completes orthogonal RHS XYZ
;         'Ygsm', translate into FAC using cartesian Ygsm position as other dimension.
;                 Y Axis on plane defined by Ygsm and Z
;                 First coordinate definition: X = Ygsm 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:
;--> Should filter NaNs to supress floating point errors
;--> Contains coordinate transformation specific code, if new
;    coordinate systems are added, this code should be updated
;
;  See also:
;     cotrans,cotrans_get_coord,tvector_rotate,minvar_matrix_make,fac_crib
;
;
; $LastChangedBy: jwl $
; $LastChangedDate: 2010-02-05 15:51:58 -0800 (Fri, 05 Feb 2010) $
; $LastChangedRevision: 7220 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/trunk/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, degap=degap, _extra=ex

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

valid_coords = ['xgse','rgeo','mrgeo','phigeo','mphigeo','phism','mphism','ygsm']

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
if keyword_set(degap) then tdegap,mag_var_name,newname=mag_var_name+'_ctv_temp', _extra=ex

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

;adding position variable validation
;2008/05/06
if keyword_set(pos_var_name) then begin
  if tnames(pos_var_name) eq '' then begin
     message,/continue,'illegal position tplot variable name'
     return
  endif

  if n_elements(tnames(pos_var_name)) ne 1 then begin
     message,/continue,'can only input a single position tplot variable'
     return
  endif

  if cotrans_get_coord(pos_var_name) ne 'gei' then begin
     message,/continue,'position tplot variable must be in gei coordinates'
     return
  endif

endif

case other_dim_l of

;rhs coordinates constructed on plane Xgse - B
valid_coords[0] : 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)))

    str_element,d,'v',SUCCESS=s
    dl_temp=dl
    dl_temp.data_att.coord_sys = 'gse'
    str_element,d,'v',SUCCESS=s

    if(s) then $
      store_data,'fac_mat_x_temp',data={x:d.x,y:x_axis,v:d.v},dlimits=dl_temp $
    else $
      store_data,'fac_mat_x_temp',data={x:d.x,y:x_axis},dlimits=dl_temp

    case dl.data_att.coord_sys of
     'gse' : ; no rotation required
     'gsm' : ; no rotation required
      else : begin
               message,/continue,' FX failed to rotate Xgse vector onto input mag field system'
               return
             end
    endcase

    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 Xgse','B']

end

;coordinates constructed on rhs coordinate system defined by B - Rgeo
valid_coords[1] : begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 Rgeo 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_z_temp','fac_mat_pos_temp',newname='fac_mat_y_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,' FX  failed to generate B x Rgeo vector'
        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 normalize y-vector after cross-product'
      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'
    del_data,pos_var_name_temp


    str_element,dl,'labels',success=s1

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

end

;coordinates constructed on rhs coordinate system defined by B - mRgeo
valid_coords[2] : begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 mRgeo 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 Rgeo x B = B x mRgeo vector'
        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 normalize y-vector after cross-product'
      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'
    del_data,pos_var_name_temp


    str_element,dl,'labels',success=s1

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

end

;begin Phigeo transformation
valid_coords[3] : begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 Phigeo 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 (it is identical to phi GEO)
    ;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


    str_element,d,'v',SUCCESS=s

    if(s) then $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v} $
    else $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2}

    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 Phigeo 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'
    del_data,pos_var_name_temp

    str_element,dl,'labels',success=s1

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

end

;begin mPhigeo transformation
valid_coords[4] : begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 mPhigeo 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 (it is identical to phi GEO)
    ;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


    str_element,d,'v',SUCCESS=s

    if(s) then $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v} $
    else $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2}

    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_z_temp','fac_mat_pos_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate mPhigeo 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'
    del_data,pos_var_name_temp

    str_element,dl,'labels',success=s1

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

end

;begin Phism transformation
valid_coords[5]: begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 Phism transformation'
        return
    endif

    cotrans,pos_var_name_temp,pos_var_name_temp,/gei2gse
    cotrans,pos_var_name_temp,pos_var_name_temp,/gse2gsm
    cotrans,pos_var_name_temp,pos_var_name_temp,/gsm2sm

    tinterpol_mxn,pos_var_name_temp,'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 SM
    ;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 SM system (overwrite y_2)
    y_2(*,0)=-sin(phi_d*!PI/180.)
    y_2(*,1)=cos(phi_d*!PI/180.)
    y_2(*,2)=0.

    str_element,d,'v',SUCCESS=s

    if(s) then $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v} $
    else $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2}

    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 Phism 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'
    del_data,pos_var_name_temp

    str_element,dl,'labels',success=s1

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

end


;begin mPhism transformation
valid_coords[6]: begin

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

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

    ;create temp variable so position vector is not overwritten
    get_data, pos_var_name, data=d_2,limits=l_2,dlimits=dl_2
    pos_var_name_temp = pos_var_name+'_temp'
    store_data, pos_var_name_temp, data=d_2,limits=l_2,dlimits=dl_2

    ;remove any gaps in the data
    if keyword_set(degap) then tdegap,pos_var_name_temp,newname=pos_var_name_temp, _extra=ex

    get_data,pos_var_name_temp,data=d_2,dlimits=dl_2

    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 mPhism transformation'
        return
    endif

    cotrans,pos_var_name_temp,pos_var_name_temp,/gei2gse
    cotrans,pos_var_name_temp,pos_var_name_temp,/gse2gsm
    cotrans,pos_var_name_temp,pos_var_name_temp,/gsm2sm

    tinterpol_mxn,pos_var_name_temp,'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 SM
    ;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 SM system (overwrite y_2)
    y_2(*,0)=-sin(phi_d*!PI/180.)
    y_2(*,1)=cos(phi_d*!PI/180.)
    y_2(*,2)=0.

    str_element,d,'v',SUCCESS=s

    if(s) then $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2, v:d.v} $
    else $
      store_data,'fac_mat_pos_temp',data={x:d.x,y:y_2}

    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_z_temp','fac_mat_pos_temp',newname='fac_mat_x_temp',error=error_state

    if(error_state eq 0) then begin
        message,/continue,'FX failed to calculate mPhism 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'
    del_data,pos_var_name_temp

    str_element,dl,'labels',success=s1

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

end



;begin Ygsm transformation (similar to other_dim=Xgse)
valid_coords[7]: 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 Ygsm 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
    y_axis = transpose(rebin([0D,1D,0D],3,n_elements(d.x)))

    str_element,d,'v',SUCCESS=s
    dl_temp=dl
    dl_temp.data_att.coord_sys = 'gsm'

    if(s) then $
      store_data,'fac_mat_y_temp',data={x:d.x,y:y_axis,v:d.v},dlimits=dl_temp $
    else $
      store_data,'fac_mat_y_temp',data={x:d.x,y:y_axis},dlimits=dl_temp

    case dl.data_att.coord_sys of
     'gse' : begin
               cotrans,'fac_mat_y_temp','fac_mat_y_temp',/gsm2gse
             end
     'gsm' : ; no rotation required
      else : begin
               message,/continue,' FX failed to rotate Ygsm vector onto input mag field system'
               return
             end
    endcase

    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 Ygsm'
        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 Ygsm in Ygsm'
        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 create zx cross product in Ygsm'
        return
    endif

    str_element,dl,'labels',success=s1

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

end

else: begin
      message,/continue,'no valid other_dimension'
               return
      end
endcase


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])


str_element,d,'v',SUCCESS=s

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

;dl.data_att.coord_sys=valid_coords[0]

str_element,dl,'data_att.source_sys',dl.data_att.coord_sys,/add
dl.data_att.coord_sys=other_dim_l

;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
