% This example simulates a simple cross flexure
% In the post-processing step, we'll plot the effective rotation axis of
% the rigid body when it moves according to various eigenmodes of the
% system.

clear
clc
% addpath('spacar')		%Specify location of spacar folder, if not current

%% NODE POSITIONS
%           x y z
nodes = [   0 0 0;      %node 1  
            0 0.1 0;    %node 2  
            0.1 0.1 0;  %node 3
            0.1 0 0];   %node 4

        
%% ELEMENT CONNECTIVITY
%               p   q
elements = [    1   3;  %element 1
                2   3;  %element 2
                2   4]; %element 3

            
%% NODE PROPERTIES  
%node 1
nprops(1).fix               = true;         %Fix node 1

%node 3
nprops(2).mass              = 0.1;          %Mass [kg] of node 2
nprops(3).mass              = 0.1;          %Mass [kg] of node 3

%node 4
nprops(4).fix               = true;         %Fix node 4


%% ELEMENT PROPERTIES
%Property set 1
eprops(1).elems    = [1 3];            %Add this set of properties to elements 1 and 3
eprops(1).emod     = 210e9;            %E-modulus [Pa]
eprops(1).smod     = 70e9;             %G-modulus [Pa]
eprops(1).dens     = 7800;             %Density [kg/m^3]
eprops(1).cshape   = 'rect';           %Rectangular cross-section
eprops(1).dim      = [50e-3 0.2e-3];   %Width: 50 mm, thickness: 0.2 mm
eprops(1).orien    = [0 0 1];          %Orientation of the cross-section as a vector pointing along "width-direction"
eprops(1).nbeams   = 2;                %Number of beams used to model this element 
eprops(1).flex     = 1:6;        	   %Model all six deformation modes as flexible
eprops(1).color    = 'grey';           %Color
eprops(1).opacity  = 0.7;              %Opacity

%Property set 2
eprops(2).elems    = 2;                %Add this set of properties to element 2
eprops(2).dens     = 3000;             %Density [kg/m^3]
eprops(2).cshape   = 'rect';           %Rectangular cross-section
eprops(2).dim      = [50e-3 10e-3];    %Width: 50 mm, thickness: 10 mm
eprops(2).orien    = [0 0 1];          %Orientation of the cross-section as a vector pointing along "width-direction"
eprops(2).nbeams   = 1;                %Single beam for simulating this element is sufficient (since it is rigid)
eprops(2).color    = 'darkblue';       %Color

%% OPTIONAL ARGUMENTS
opt.filename    = 'ex1';     %Filename

%% CALL SPACAR_LIGHT
out = spacarlight(nodes, elements, nprops, eprops, opt);

%%%%%%%%%%%%%
%%% plotting the effective rotation axis of a rigid body, connected to node
%%% i_node, when it moves according to eigenmode with number i_mode.
%%%%%%%%%%%%

fh = out.fighandle; %handle to the Spavisual window
ah = fh.Children(1); %handle to the axis in the Spavisual window
% view(ah,70,30) %set the desired view angle

i_node = 3; %node number to which the rigid body of interest is connected
i_mode = 1; %number of the eigenmode of interest

%get rotation axis
[screwaxis,screwloc] = getPoleAxis(opt.filename,i_node,i_mode);

%mark the node that was provided, using a black dot
plot3(ah,nodes(i_node,1),nodes(i_node,2),nodes(i_node,3), 'k.', 'MarkerSize', 30);

%mark the point on the rotation axis that is *closest to the provided node*
plot3(ah,screwloc(1),screwloc(2),screwloc(3), 'r.', 'MarkerSize', 30);

%plot the rotation axis corresponding with the rigid body motion occuring in
% the eigenmode with number i_mode
ax = [screwloc screwloc+screwaxis screwloc-screwaxis];
plot3(ah,ax(1,:),ax(2,:),ax(3,:),'r')