Subsystems used for the Simscape Models

function [ref] = initializeReferences(opts_param)
    %% Default values for opts
    opts = struct(   ...
        'Ts',           1e-3, ...       % Sampling Frequency [s]
        'Dy_type',      'constant', ... % Either "constant" / "triangular" / "sinusoidal"
        'Dy_amplitude', 0, ...       % Amplitude of the displacement [m]
        'Dy_period',    1, ...          % Period of the displacement [s]
        'Ry_type',      'constant', ... % Either "constant" / "triangular" / "sinusoidal"
        'Ry_amplitude', 0, ...          % Amplitude [rad]
        'Ry_period',    10, ...         % Period of the displacement [s]
        'Rz_type',      'constant', ... % Either "constant" / "rotating"
        'Rz_amplitude', 0, ...          % Initial angle [rad]
        'Rz_period',    1, ...          % Period of the rotating [s]
        'Dh_type',      'constant', ... % For now, only constant is implemented
        'Dh_pos',       [0; 0; 0; 0; 0; 0], ...  % Initial position [m,m,m,rad,rad,rad] of the top platform (Pitch-Roll-Yaw Euler angles)
        'Rm_type',      'constant', ... % For now, only constant is implemented
        'Rm_pos',       [0; pi],  ...   % Initial position of the two masses
        'Dn_type',      'constant', ... % For now, only constant is implemented
        'Dn_pos',       [0; 0; 0; 0; 0; 0]  ...  % Initial position [m,m,m,rad,rad,rad] of the top platform
    );

    %% Populate opts with input parameters
    if exist('opts_param','var')
        for opt = fieldnames(opts_param)'
            opts.(opt{1}) = opts_param.(opt{1});
        end
    end

    %% Set Sampling Time
    Ts = opts.Ts;

    %% Translation stage - Dy
    t = 0:Ts:opts.Dy_period-Ts; % Time Vector [s]
    Dy = zeros(length(t), 1);
    switch opts.Dy_type
      case 'constant'
        Dy(:) = opts.Dy_amplitude;
      case 'triangular'
        Dy(:)                     = -4*opts.Dy_amplitude + 4*opts.Dy_amplitude/opts.Dy_period*t;
        Dy(t<0.75*opts.Dy_period) =  2*opts.Dy_amplitude - 4*opts.Dy_amplitude/opts.Dy_period*t(t<0.75*opts.Dy_period);
        Dy(t<0.25*opts.Dy_period) =  4*opts.Dy_amplitude/opts.Dy_period*t(t<0.25*opts.Dy_period);
      case 'sinusoidal'
        Dy(:) = opts.Dy_amplitude*sin(2*pi/opts.Dy_period*t);
      otherwise
        warning('Dy_type is not set correctly');
    end
    Dy = struct('time', t, 'signals', struct('values', Dy));


    %% Tilt Stage - Ry
    t = 0:Ts:opts.Ry_period-Ts;
    Ry = zeros(length(t), 1);

    switch opts.Ry_type
      case 'constant'
        Ry(:) = opts.Ry_amplitude;
      case 'triangular'
        Ry(:)                     = -4*opts.Ry_amplitude + 4*opts.Ry_amplitude/opts.Ry_period*t;
        Ry(t<0.75*opts.Ry_period) =  2*opts.Ry_amplitude - 4*opts.Ry_amplitude/opts.Ry_period*t(t<0.75*opts.Ry_period);
        Ry(t<0.25*opts.Ry_period) =  4*opts.Ry_amplitude/opts.Ry_period*t(t<0.25*opts.Ry_period);
      case 'sinusoidal'

      otherwise
        warning('Ry_type is not set correctly');
    end
    Ry = struct('time', t, 'signals', struct('values', Ry));

    %% Spindle - Rz
    t = 0:Ts:opts.Rz_period-Ts;
    Rz = zeros(length(t), 1);

    switch opts.Rz_type
      case 'constant'
        Rz(:) = opts.Rz_amplitude;
      case 'rotating'
        Rz(:) = opts.Rz_amplitude+2*pi/opts.Rz_period*t;
      otherwise
        warning('Rz_type is not set correctly');
    end
    Rz = struct('time', t, 'signals', struct('values', Rz));

    %% Micro-Hexapod
    t = [0, Ts];
    Dh = zeros(length(t), 6);
    Dhl = zeros(length(t), 6);

    switch opts.Dh_type
      case 'constant'
        Dh = [opts.Dh_pos, opts.Dh_pos];

        load('./mat/stages.mat', 'micro_hexapod');

        AP = [opts.Dh_pos(1) ; opts.Dh_pos(2) ; opts.Dh_pos(3)];

        tx = opts.Dh_pos(4);
        ty = opts.Dh_pos(5);
        tz = opts.Dh_pos(6);

        ARB = [cos(tz) -sin(tz) 0;
               sin(tz)  cos(tz) 0;
               0        0       1]*...
              [ cos(ty) 0 sin(ty);
               0        1 0;
               -sin(ty) 0 cos(ty)]*...
              [1 0        0;
               0 cos(tx) -sin(tx);
               0 sin(tx)  cos(tx)];

        [Dhl] = inverseKinematicsHexapod(micro_hexapod, AP, ARB);
        Dhl = [Dhl, Dhl];
      otherwise
        warning('Dh_type is not set correctly');
    end
    Dh = struct('time', t, 'signals', struct('values', Dh));
    Dhl = struct('time', t, 'signals', struct('values', Dhl));

    %% Axis Compensation - Rm
    t = [0, Ts];
    Rm = [opts.Rm_pos, opts.Rm_pos];
    Rm = struct('time', t, 'signals', struct('values', Rm));

    %% Nano-Hexapod
    t = [0, Ts];
    Dn = zeros(length(t), 6);

    switch opts.Dn_type
      case 'constant'
        Dn = [opts.Dn_pos, opts.Dn_pos];
      otherwise
        warning('Dn_type is not set correctly');
    end
    Dn = struct('time', t, 'signals', struct('values', Dn));

    %% Save
    save('./mat/nass_references.mat', 'Dy', 'Ry', 'Rz', 'Dh', 'Dhl', 'Rm', 'Dn', 'Ts');
end

function [L] = inverseKinematicsHexapod(hexapod, AP, ARB)
% inverseKinematicsHexapod - Compute the initial position of each leg to have the wanted Hexapod's position
%
% Syntax: inverseKinematicsHexapod(hexapod, AP, ARB)
%
% Inputs:
%    - hexapod - Hexapod object containing the geometry of the hexapod
%    - AP - Position vector of point OB expressed in frame {A} in [m]
%    - ARB - Rotation Matrix expressed in frame {A}

  % Wanted Length of the hexapod's legs [m]
  L = zeros(6, 1);

  for i = 1:length(L)
    Bbi = hexapod.pos_top_tranform(i, :)' - 1e-3*[0 ; 0 ; hexapod.TP.thickness+hexapod.Leg.sphere.top+hexapod.SP.thickness.top+hexapod.jacobian]; % [m]
    Aai = hexapod.pos_base(i, :)' + 1e-3*[0 ; 0 ; hexapod.BP.thickness+hexapod.Leg.sphere.bottom+hexapod.SP.thickness.bottom-hexapod.h-hexapod.jacobian]; % [m]

    L(i) = sqrt(AP'*AP + Bbi'*Bbi + Aai'*Aai - 2*AP'*Aai + 2*AP'*(ARB*Bbi) - 2*(ARB*Bbi)'*Aai);
  end
end

function [ground] = initializeGround()
    %%
    ground = struct();

    ground.shape = [2, 2, 0.5]; % [m]
    ground.density = 2800; % [kg/m3]
    ground.color = [0.5, 0.5, 0.5];

    %% Save
    save('./mat/stages.mat', 'ground', '-append');
end

function [granite] = initializeGranite(opts_param)
    %% Default values for opts
    opts = struct('rigid', false);

    %% Populate opts with input parameters
    if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
        opts.(opt{1}) = opts_param.(opt{1});
      end
    end

    %%
    granite = struct();

    %% Static Properties
    granite.density = 2800; % [kg/m3]
    granite.volume  = 0.749; % [m3] TODO - should
    granite.mass    = granite.density*granite.volume; % [kg]
    granite.color   = [1 1 1];
    granite.STEP    = './STEPS/granite/granite.STEP';

    granite.mass_top = 4000; % [kg] TODO

    %% Dynamical Properties
    if opts.rigid
      granite.k.x = 1e12; % [N/m]
      granite.k.y = 1e12; % [N/m]
      granite.k.z = 1e12; % [N/m]
      granite.k.rx = 1e10; % [N*m/deg]
      granite.k.ry = 1e10; % [N*m/deg]
      granite.k.rz = 1e10; % [N*m/deg]
    else
      granite.k.x = 4e9; % [N/m]
      granite.k.y = 3e8; % [N/m]
      granite.k.z = 8e8; % [N/m]
      granite.k.rx = 1e4; % [N*m/deg]
      granite.k.ry = 1e4; % [N*m/deg]
      granite.k.rz = 1e6; % [N*m/deg]
    end

    granite.c.x = 0.1*sqrt(granite.mass_top*granite.k.x); % [N/(m/s)]
    granite.c.y = 0.1*sqrt(granite.mass_top*granite.k.y); % [N/(m/s)]
    granite.c.z = 0.5*sqrt(granite.mass_top*granite.k.z); % [N/(m/s)]
    granite.c.rx = 0.1*sqrt(granite.mass_top*granite.k.rx); % [N*m/(deg/s)]
    granite.c.ry = 0.1*sqrt(granite.mass_top*granite.k.ry); % [N*m/(deg/s)]
    granite.c.rz = 0.1*sqrt(granite.mass_top*granite.k.rz); % [N*m/(deg/s)]

    %% Positioning parameters
    granite.sample_pos = 0.8; % Z-offset for the initial position of the sample [m]

    %% Save
    save('./mat/stages.mat', 'granite', '-append');
end

function [ty] = initializeTy(opts_param)
    %% Default values for opts
    opts = struct('rigid', false);

    %% Populate opts with input parameters
    if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
        opts.(opt{1}) = opts_param.(opt{1});
      end
    end

    %%
    ty = struct();

    %% Y-Translation - Static Properties
    % Ty Granite frame
    ty.granite_frame.density = 7800; % [kg/m3]
    ty.granite_frame.color   = [0.753 1 0.753];
    ty.granite_frame.STEP    = './STEPS/Ty/Ty_Granite_Frame.STEP';
    % Guide Translation Ty
    ty.guide.density         = 7800; % [kg/m3]
    ty.guide.color           = [0.792 0.820 0.933];
    ty.guide.STEP            = './STEPS/ty/Ty_Guide.STEP';
    % Ty - Guide_Translation12
    ty.guide12.density       = 7800; % [kg/m3]
    ty.guide12.color         = [0.792 0.820 0.933];
    ty.guide12.STEP          = './STEPS/Ty/Ty_Guide_12.STEP';
    % Ty - Guide_Translation11
    ty.guide11.density       = 7800; % [kg/m3]
    ty.guide11.color         = [0.792 0.820 0.933];
    ty.guide11.STEP          = './STEPS/ty/Ty_Guide_11.STEP';
    % Ty - Guide_Translation22
    ty.guide22.density       = 7800; % [kg/m3]
    ty.guide22.color         = [0.792 0.820 0.933];
    ty.guide22.STEP          = './STEPS/ty/Ty_Guide_22.STEP';
    % Ty - Guide_Translation21
    ty.guide21.density       = 7800; % [kg/m3]
    ty.guide21.color         = [0.792 0.820 0.933];
    ty.guide21.STEP          = './STEPS/Ty/Ty_Guide_21.STEP';
    % Ty - Plateau translation
    ty.frame.density         = 7800; % [kg/m3]
    ty.frame.color           = [0.792 0.820 0.933];
    ty.frame.STEP            = './STEPS/ty/Ty_Stage.STEP';
    % Ty Stator Part
    ty.stator.density        = 5400; % [kg/m3]
    ty.stator.color          = [0.792 0.820 0.933];
    ty.stator.STEP           = './STEPS/ty/Ty_Motor_Stator.STEP';
    % Ty Rotor Part
    ty.rotor.density         = 5400; % [kg/m3]
    ty.rotor.color           = [0.792 0.820 0.933];
    ty.rotor.STEP            = './STEPS/ty/Ty_Motor_Rotor.STEP';

    ty.m = 1000; % TODO [kg]

    %% Y-Translation - Dynamicals Properties
    if opts.rigid
      ty.k.ax  = 1e12; % Axial Stiffness for each of the 4 guidance (y) [N/m]
      ty.k.rad = 1e12; % Radial Stiffness for each of the 4 guidance (x-z) [N/m]
    else
      ty.k.ax  = 5e8; % Axial Stiffness for each of the 4 guidance (y) [N/m]
      ty.k.rad = 5e7; % Radial Stiffness for each of the 4 guidance (x-z) [N/m]
    end

    ty.c.ax  = 0.1*sqrt(ty.k.ax*ty.m);
    ty.c.rad = 0.1*sqrt(ty.k.rad*ty.m);

    %% Save
    save('./mat/stages.mat', 'ty', '-append');
end

function [ry] = initializeRy(opts_param)
    %% Default values for opts
    opts = struct('rigid', false);

    %% Populate opts with input parameters
    if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
        opts.(opt{1}) = opts_param.(opt{1});
      end
    end

    %%
    ry = struct();

    %% Tilt Stage - Static Properties
    % Ry - Guide for the tilt stage
    ry.guide.density = 7800; % [kg/m3]
    ry.guide.color   = [0.792 0.820 0.933];
    ry.guide.STEP    = './STEPS/ry/Tilt_Guide.STEP';
    % Ry - Rotor of the motor
    ry.rotor.density = 2400; % [kg/m3]
    ry.rotor.color   = [0.792 0.820 0.933];
    ry.rotor.STEP    = './STEPS/ry/Tilt_Motor_Axis.STEP';
    % Ry - Motor
    ry.motor.density = 3200; % [kg/m3]
    ry.motor.color   = [0.792 0.820 0.933];
    ry.motor.STEP    = './STEPS/ry/Tilt_Motor.STEP';
    % Ry - Plateau Tilt
    ry.stage.density = 7800; % [kg/m3]
    ry.stage.color   = [0.792 0.820 0.933];
    ry.stage.STEP    = './STEPS/ry/Tilt_Stage.STEP';

    ry.m = 800; % TODO [kg]

    %% Tilt Stage - Dynamical Properties
    if opts.rigid
      ry.k.tilt = 1e10; % Rotation stiffness around y [N*m/deg]
      ry.k.h    = 1e12; % Stiffness in the direction of the guidance [N/m]
      ry.k.rad  = 1e12; % Stiffness in the top direction [N/m]
      ry.k.rrad = 1e12; % Stiffness in the side direction [N/m]
    else
      ry.k.tilt = 1e4; % Rotation stiffness around y [N*m/deg]
      ry.k.h    = 1e8; % Stiffness in the direction of the guidance [N/m]
      ry.k.rad  = 1e8; % Stiffness in the top direction [N/m]
      ry.k.rrad = 1e8; % Stiffness in the side direction [N/m]
    end

    ry.c.h    = 0.1*sqrt(ry.k.h*ry.m);
    ry.c.rad  = 0.1*sqrt(ry.k.rad*ry.m);
    ry.c.rrad = 0.1*sqrt(ry.k.rrad*ry.m);
    ry.c.tilt = 0.1*sqrt(ry.k.tilt*ry.m);

    %% Positioning parameters
    ry.z_offset = 0.58178; % Z-Offset so that the center of rotation matches the sample center [m]

    %% Save
    save('./mat/stages.mat', 'ry', '-append');
end

function [rz] = initializeRz(opts_param)
    %% Default values for opts
    opts = struct('rigid', false);

    %% Populate opts with input parameters
    if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
        opts.(opt{1}) = opts_param.(opt{1});
      end
    end

    %%
    rz = struct();

    %% Spindle - Static Properties
    % Spindle - Slip Ring
    rz.slipring.density = 7800; % [kg/m3]
    rz.slipring.color   = [0.792 0.820 0.933];
    rz.slipring.STEP    = './STEPS/rz/Spindle_Slip_Ring.STEP';
    % Spindle - Rotor
    rz.rotor.density    = 7800; % [kg/m3]
    rz.rotor.color      = [0.792 0.820 0.933];
    rz.rotor.STEP       = './STEPS/rz/Spindle_Rotor.STEP';
    % Spindle - Stator
    rz.stator.density   = 7800; % [kg/m3]
    rz.stator.color     = [0.792 0.820 0.933];
    rz.stator.STEP      = './STEPS/rz/Spindle_Stator.STEP';

    % Estimated mass of the mooving part
    rz.m = 250; % [kg]

    %% Spindle - Dynamical Properties

    if opts.rigid
      rz.k.rot  = 1e10; % Rotational Stiffness (Rz) [N*m/deg]
      rz.k.tilt = 1e10; % Rotational Stiffness (Rx, Ry) [N*m/deg]
      rz.k.ax   = 1e12; % Axial Stiffness (Z) [N/m]
      rz.k.rad  = 1e12; % Radial Stiffness (X, Y) [N/m]
    else
      rz.k.rot  = 1e6; % TODO - Rotational Stiffness (Rz) [N*m/deg]
      rz.k.tilt = 1e6; % Rotational Stiffness (Rx, Ry) [N*m/deg]
      rz.k.ax   = 2e9; % Axial Stiffness (Z) [N/m]
      rz.k.rad  = 7e8; % Radial Stiffness (X, Y) [N/m]
    end

    % Damping
    rz.c.ax   = 0.1*sqrt(rz.k.ax*rz.m);
    rz.c.rad  = 0.1*sqrt(rz.k.rad*rz.m);
    rz.c.tilt = 0.1*sqrt(rz.k.tilt*rz.m);
    rz.c.rot  = 0.1*sqrt(rz.k.rot*rz.m);

    %% Save
    save('./mat/stages.mat', 'rz', '-append');
end

function [micro_hexapod] = initializeMicroHexapod(opts_param)
    %% Default values for opts
    opts = struct(...
      'rigid', false, ...
      'AP',    zeros(3, 1), ... % Wanted position in [m] of OB with respect to frame {A}
      'ARB',   eye(3) ...       % Rotation Matrix that represent the wanted orientation of frame {B} with respect to frame {A}
    );

    %% Populate opts with input parameters
    if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
        opts.(opt{1}) = opts_param.(opt{1});
      end
    end

    %% Stewart Object
    micro_hexapod = struct();
    micro_hexapod.h        = 350; % Total height of the platform [mm]
    micro_hexapod.jacobian = 270; % Distance from the top of the mobile platform to the Jacobian point [mm]

    %% Bottom Plate - Mechanical Design
    BP = struct();

    BP.rad.int   = 110;   % Internal Radius [mm]
    BP.rad.ext   = 207.5; % External Radius [mm]
    BP.thickness = 26;    % Thickness [mm]
    BP.leg.rad   = 175.5; % Radius where the legs articulations are positionned [mm]
    BP.leg.ang   = 9.5;   % Angle Offset [deg]
    BP.density   = 8000;  % Density of the material [kg/m^3]
    BP.color     = [0.6 0.6 0.6]; % Color [rgb]
    BP.shape     = [BP.rad.int BP.thickness; BP.rad.int 0; BP.rad.ext 0; BP.rad.ext BP.thickness];

    %% Top Plate - Mechanical Design
    TP = struct();

    TP.rad.int   = 82;   % Internal Radius [mm]
    TP.rad.ext   = 150;  % Internal Radius [mm]
    TP.thickness = 26;   % Thickness [mm]
    TP.leg.rad   = 118;  % Radius where the legs articulations are positionned [mm]
    TP.leg.ang   = 12.1; % Angle Offset [deg]
    TP.density   = 8000; % Density of the material [kg/m^3]
    TP.color     = [0.6 0.6 0.6]; % Color [rgb]
    TP.shape     = [TP.rad.int TP.thickness; TP.rad.int 0; TP.rad.ext 0; TP.rad.ext TP.thickness];

    %% Struts
    Leg = struct();

    Leg.stroke     = 10e-3; % Maximum Stroke of each leg [m]
    if opts.rigid
      Leg.k.ax     = 1e12; % Stiffness of each leg [N/m]
    else
      Leg.k.ax     = 2e7; % Stiffness of each leg [N/m]
    end
    Leg.ksi.ax     = 0.1;   % Modal damping ksi = 1/2*c/sqrt(km) []
    Leg.rad.bottom = 25;    % Radius of the cylinder of the bottom part [mm]
    Leg.rad.top    = 17;    % Radius of the cylinder of the top part [mm]
    Leg.density    = 8000;  % Density of the material [kg/m^3]
    Leg.color.bottom  = [0.5 0.5 0.5]; % Color [rgb]
    Leg.color.top     = [0.5 0.5 0.5]; % Color [rgb]

    Leg.sphere.bottom = Leg.rad.bottom; % Size of the sphere at the end of the leg [mm]
    Leg.sphere.top    = Leg.rad.top; % Size of the sphere at the end of the leg [mm]
    Leg.m             = TP.density*((pi*(TP.rad.ext/1000)^2)*(TP.thickness/1000)-(pi*(TP.rad.int/1000^2))*(TP.thickness/1000))/6; % TODO [kg]
    Leg = updateDamping(Leg);


    %% Sphere
    SP = struct();

    SP.height.bottom  = 27; % [mm]
    SP.height.top     = 27; % [mm]
    SP.density.bottom = 8000; % [kg/m^3]
    SP.density.top    = 8000; % [kg/m^3]
    SP.color.bottom   = [0.6 0.6 0.6]; % [rgb]
    SP.color.top      = [0.6 0.6 0.6]; % [rgb]
    SP.k.ax           = 0; % [N*m/deg]
    SP.ksi.ax         = 10;

    SP.thickness.bottom = SP.height.bottom-Leg.sphere.bottom; % [mm]
    SP.thickness.top    = SP.height.top-Leg.sphere.top; % [mm]
    SP.rad.bottom       = Leg.sphere.bottom; % [mm]
    SP.rad.top          = Leg.sphere.top; % [mm]
    SP.m                = SP.density.bottom*2*pi*((SP.rad.bottom*1e-3)^2)*(SP.height.bottom*1e-3); % TODO [kg]

    SP = updateDamping(SP);

    %%
    Leg.support.bottom = [0 SP.thickness.bottom; 0 0; SP.rad.bottom 0; SP.rad.bottom SP.height.bottom];
    Leg.support.top    = [0 SP.thickness.top; 0 0; SP.rad.top 0; SP.rad.top SP.height.top];

    %%
    micro_hexapod.BP = BP;
    micro_hexapod.TP = TP;
    micro_hexapod.Leg = Leg;
    micro_hexapod.SP = SP;

    %%
    micro_hexapod = initializeParameters(micro_hexapod);

    %% Setup equilibrium position of each leg
    micro_hexapod.L0 = inverseKinematicsHexapod(micro_hexapod, opts.AP, opts.ARB);

    %% Save
    save('./mat/stages.mat', 'micro_hexapod', '-append');

    %%
    function [element] = updateDamping(element)
      field = fieldnames(element.k);
      for i = 1:length(field)
        element.c.(field{i}) = 2*element.ksi.(field{i})*sqrt(element.k.(field{i})*element.m);
      end
    end

    %%
    function [stewart] = initializeParameters(stewart)
        %% Connection points on base and top plate w.r.t. World frame at the center of the base plate
        stewart.pos_base = zeros(6, 3);
        stewart.pos_top = zeros(6, 3);

        alpha_b = stewart.BP.leg.ang*pi/180; % angle de décalage par rapport à 120 deg (pour positionner les supports bases)
        alpha_t = stewart.TP.leg.ang*pi/180; % +- offset angle from 120 degree spacing on top

        height = (stewart.h-stewart.BP.thickness-stewart.TP.thickness-stewart.Leg.sphere.bottom-stewart.Leg.sphere.top-stewart.SP.thickness.bottom-stewart.SP.thickness.top)*0.001; % TODO

        radius_b = stewart.BP.leg.rad*0.001; % rayon emplacement support base
        radius_t = stewart.TP.leg.rad*0.001; % top radius in meters

        for i = 1:3
          % base points
          angle_m_b = (2*pi/3)* (i-1) - alpha_b;
          angle_p_b = (2*pi/3)* (i-1) + alpha_b;
          stewart.pos_base(2*i-1,:) =  [radius_b*cos(angle_m_b), radius_b*sin(angle_m_b), 0.0];
          stewart.pos_base(2*i,:) = [radius_b*cos(angle_p_b), radius_b*sin(angle_p_b), 0.0];

          % top points
          % Top points are 60 degrees offset
          angle_m_t = (2*pi/3)* (i-1) - alpha_t + 2*pi/6;
          angle_p_t = (2*pi/3)* (i-1) + alpha_t + 2*pi/6;
          stewart.pos_top(2*i-1,:) = [radius_t*cos(angle_m_t), radius_t*sin(angle_m_t), height];
          stewart.pos_top(2*i,:) = [radius_t*cos(angle_p_t), radius_t*sin(angle_p_t), height];
        end

        % permute pos_top points so that legs are end points of base and top points
        stewart.pos_top = [stewart.pos_top(6,:); stewart.pos_top(1:5,:)]; %6th point on top connects to 1st on bottom
        stewart.pos_top_tranform = stewart.pos_top - height*[zeros(6, 2),ones(6, 1)];

        %% leg vectors
        legs = stewart.pos_top - stewart.pos_base;
        leg_length = zeros(6, 1);
        leg_vectors = zeros(6, 3);
        for i = 1:6
          leg_length(i) = norm(legs(i,:));
          leg_vectors(i,:)  = legs(i,:) / leg_length(i);
        end

        stewart.Leg.lenght = 1000*leg_length(1)/1.5;
        stewart.Leg.shape.bot = [0 0; ...
                                stewart.Leg.rad.bottom 0; ...
                                stewart.Leg.rad.bottom stewart.Leg.lenght; ...
                                stewart.Leg.rad.top stewart.Leg.lenght; ...
                                stewart.Leg.rad.top 0.2*stewart.Leg.lenght; ...
                                0 0.2*stewart.Leg.lenght];

        %% Calculate revolute and cylindrical axes
        rev1 = zeros(6, 3);
        rev2 = zeros(6, 3);
        cyl1 = zeros(6, 3);
        for i = 1:6
          rev1(i,:) = cross(leg_vectors(i,:), [0 0 1]);
          rev1(i,:) = rev1(i,:) / norm(rev1(i,:));

          rev2(i,:) = - cross(rev1(i,:), leg_vectors(i,:));
          rev2(i,:) = rev2(i,:) / norm(rev2(i,:));

          cyl1(i,:) = leg_vectors(i,:);
        end


        %% Coordinate systems
        stewart.lower_leg = struct('rotation', eye(3));
        stewart.upper_leg = struct('rotation', eye(3));

        for i = 1:6
          stewart.lower_leg(i).rotation = [rev1(i,:)', rev2(i,:)', cyl1(i,:)'];
          stewart.upper_leg(i).rotation = [rev1(i,:)', rev2(i,:)', cyl1(i,:)'];
        end

        %% Position Matrix
        stewart.M_pos_base = stewart.pos_base + (height+(stewart.TP.thickness+stewart.Leg.sphere.top+stewart.SP.thickness.top+stewart.jacobian)*1e-3)*[zeros(6, 2),ones(6, 1)];

        %% Compute Jacobian Matrix
        aa = stewart.pos_top_tranform + (stewart.jacobian - stewart.TP.thickness - stewart.SP.height.top)*1e-3*[zeros(6, 2),ones(6, 1)];
        stewart.J  = getJacobianMatrix(leg_vectors', aa');
    end

    %%
    function J  = getJacobianMatrix(RM, M_pos_base)
        % RM: [3x6] unit vector of each leg in the fixed frame
        % M_pos_base: [3x6] vector of the leg connection at the top platform location in the fixed frame
        J = zeros(6);
        J(:, 1:3) = RM';
        J(:, 4:6) = cross(M_pos_base, RM)';
    end
end

function [axisc] = initializeAxisc()
    %%
    axisc = struct();

    %% Axis Compensator - Static Properties
    % Structure
    axisc.structure.density = 3400; % [kg/m3]
    axisc.structure.color   = [0.792 0.820 0.933];
    axisc.structure.STEP    = './STEPS/axisc/axisc_structure.STEP';
    % Wheel
    axisc.wheel.density     = 2700; % [kg/m3]
    axisc.wheel.color       = [0.753 0.753 0.753];
    axisc.wheel.STEP        = './STEPS/axisc/axisc_wheel.STEP';
    % Mass
    axisc.mass.density      = 7800; % [kg/m3]
    axisc.mass.color        = [0.792 0.820 0.933];
    axisc.mass.STEP         = './STEPS/axisc/axisc_mass.STEP';
    % Gear
    axisc.gear.density      = 7800; % [kg/m3]
    axisc.gear.color        = [0.792 0.820 0.933];
    axisc.gear.STEP         = './STEPS/axisc/axisc_gear.STEP';

    axisc.m      = 40; % TODO [kg]

    %% Axis Compensator - Dynamical Properties
    % axisc.k.ax   = 1; % TODO [N*m/deg)]
    % axisc.c.ax   = (1/1)*sqrt(axisc.k.ax/axisc.m);

    %% Save
    save('./mat/stages.mat', 'axisc', '-append');
end

function [] = initializeMirror(opts_param)
    %% Default values for opts
    opts = struct(...
        'shape', 'spherical', ... % spherical or conical
        'angle', 45 ...
    );

    %% Populate opts with input parameters
    if exist('opts_param','var')
        for opt = fieldnames(opts_param)'
            opts.(opt{1}) = opts_param.(opt{1});
        end
    end

    %%
    mirror = struct();
    mirror.h = 50; % height of the mirror [mm]
    mirror.thickness = 25; % Thickness of the plate supporting the sample [mm]
    mirror.hole_rad = 120; % radius of the hole in the mirror [mm]
    mirror.support_rad = 100; % radius of the support plate [mm]
    mirror.jacobian = 150; % point of interest offset in z (above the top surfave) [mm]
    mirror.rad = 180; % radius of the mirror (at the bottom surface) [mm]

    mirror.density = 2400; % Density of the mirror [kg/m3]
    mirror.color = [0.4 1.0 1.0]; % Color of the mirror

    mirror.cone_length = mirror.rad*tand(opts.angle)+mirror.h+mirror.jacobian; % Distance from Apex point of the cone to jacobian point

    %% Shape
    mirror.shape = [...
        0 mirror.h-mirror.thickness
        mirror.hole_rad mirror.h-mirror.thickness; ...
        mirror.hole_rad 0; ...
        mirror.rad 0 ...
    ];

    if strcmp(opts.shape, 'spherical')
        mirror.sphere_radius = sqrt((mirror.jacobian+mirror.h)^2+mirror.rad^2); % Radius of the sphere [mm]

        for z = linspace(0, mirror.h, 101)
            mirror.shape = [mirror.shape; sqrt(mirror.sphere_radius^2-(z-mirror.jacobian-mirror.h)^2) z];
        end
    elseif strcmp(opts.shape, 'conical')
        mirror.shape = [mirror.shape; mirror.rad+mirror.h/tand(opts.angle) mirror.h];
    else
        error('Shape should be either conical or spherical');
    end

    mirror.shape = [mirror.shape; 0 mirror.h];

    %% Save
    save('./mat/stages.mat', 'mirror', '-append');
end

function [nano_hexapod] = initializeNanoHexapod(opts_param)
    %% Default values for opts
    opts = struct('actuator', 'piezo');

    %% Populate opts with input parameters
    if exist('opts_param','var')
        for opt = fieldnames(opts_param)'
            opts.(opt{1}) = opts_param.(opt{1});
        end
    end

    %% Stewart Object
    nano_hexapod = struct();
    nano_hexapod.h        = 90;  % Total height of the platform [mm]
    nano_hexapod.jacobian = 175; % Point where the Jacobian is computed => Center of rotation [mm]
%     nano_hexapod.jacobian = 174.26; % Point where the Jacobian is computed => Center of rotation [mm]

    %% Bottom Plate
    BP = struct();

    BP.rad.int   = 0;   % Internal Radius [mm]
    BP.rad.ext   = 150; % External Radius [mm]
    BP.thickness = 10;  % Thickness [mm]
    BP.leg.rad   = 100; % Radius where the legs articulations are positionned [mm]
    BP.leg.ang   = 5;   % Angle Offset [deg]
    BP.density   = 8000;% Density of the material [kg/m^3]
    BP.color     = [0.7 0.7 0.7]; % Color [rgb]
    BP.shape     = [BP.rad.int BP.thickness; BP.rad.int 0; BP.rad.ext 0; BP.rad.ext BP.thickness];

    %% Top Plate
    TP = struct();

    TP.rad.int   = 0;   % Internal Radius [mm]
    TP.rad.ext   = 100; % Internal Radius [mm]
    TP.thickness = 10;  % Thickness [mm]
    TP.leg.rad   = 90;  % Radius where the legs articulations are positionned [mm]
    TP.leg.ang   = 5;   % Angle Offset [deg]
    TP.density   = 8000;% Density of the material [kg/m^3]
    TP.color     = [0.7 0.7 0.7]; % Color [rgb]
    TP.shape     = [TP.rad.int TP.thickness; TP.rad.int 0; TP.rad.ext 0; TP.rad.ext TP.thickness];

    %% Leg
    Leg = struct();

    Leg.stroke     = 80e-6; % Maximum Stroke of each leg [m]
    if strcmp(opts.actuator, 'piezo')
        Leg.k.ax   = 1e7;   % Stiffness of each leg [N/m]
    elseif strcmp(opts.actuator, 'lorentz')
        Leg.k.ax   = 1e4;   % Stiffness of each leg [N/m]
    else
        error('opts.actuator should be piezo or lorentz');
    end
    Leg.ksi.ax     = 10;    % Maximum amplification at resonance []
    Leg.rad.bottom = 12;    % Radius of the cylinder of the bottom part [mm]
    Leg.rad.top    = 10;    % Radius of the cylinder of the top part [mm]
    Leg.density    = 8000;  % Density of the material [kg/m^3]
    Leg.color.bottom  = [0.5 0.5 0.5]; % Color [rgb]
    Leg.color.top     = [0.5 0.5 0.5]; % Color [rgb]

    Leg.sphere.bottom = Leg.rad.bottom; % Size of the sphere at the end of the leg [mm]
    Leg.sphere.top    = Leg.rad.top; % Size of the sphere at the end of the leg [mm]
    Leg.m             = TP.density*((pi*(TP.rad.ext/1000)^2)*(TP.thickness/1000)-(pi*(TP.rad.int/1000^2))*(TP.thickness/1000))/6; % TODO [kg]

    Leg = updateDamping(Leg);


    %% Sphere
    SP = struct();

    SP.height.bottom  = 15; % [mm]
    SP.height.top     = 15; % [mm]
    SP.density.bottom = 8000; % [kg/m^3]
    SP.density.top    = 8000; % [kg/m^3]
    SP.color.bottom   = [0.7 0.7 0.7]; % [rgb]
    SP.color.top      = [0.7 0.7 0.7]; % [rgb]
    SP.k.ax           = 0; % [N*m/deg]
    SP.ksi.ax         = 0;

    SP.thickness.bottom = SP.height.bottom-Leg.sphere.bottom; % [mm]
    SP.thickness.top    = SP.height.top-Leg.sphere.top; % [mm]
    SP.rad.bottom       = Leg.sphere.bottom; % [mm]
    SP.rad.top          = Leg.sphere.top; % [mm]
    SP.m                = SP.density.bottom*2*pi*((SP.rad.bottom*1e-3)^2)*(SP.height.bottom*1e-3); % TODO [kg]

    SP = updateDamping(SP);

    %%
    Leg.support.bottom = [0 SP.thickness.bottom; 0 0; SP.rad.bottom 0; SP.rad.bottom SP.height.bottom];
    Leg.support.top    = [0 SP.thickness.top; 0 0; SP.rad.top 0; SP.rad.top SP.height.top];

    %%
    nano_hexapod.BP = BP;
    nano_hexapod.TP = TP;
    nano_hexapod.Leg = Leg;
    nano_hexapod.SP = SP;

    %%
    nano_hexapod = initializeParameters(nano_hexapod);

    %% Save
    save('./mat/stages.mat', 'nano_hexapod', '-append');

    %%
    function [element] = updateDamping(element)
        field = fieldnames(element.k);
        for i = 1:length(field)
            if element.ksi.(field{i}) > 0
              element.c.(field{i}) = 1/element.ksi.(field{i})*sqrt(element.k.(field{i})/element.m);
            else
              element.c.(field{i}) = 0;
            end
        end
    end

    %%
    function [stewart] = initializeParameters(stewart)
        %% Connection points on base and top plate w.r.t. World frame at the center of the base plate
        stewart.pos_base = zeros(6, 3);
        stewart.pos_top = zeros(6, 3);

        alpha_b = stewart.BP.leg.ang*pi/180; % angle de décalage par rapport à 120 deg (pour positionner les supports bases)
        alpha_t = stewart.TP.leg.ang*pi/180; % +- offset angle from 120 degree spacing on top

        height = (stewart.h-stewart.BP.thickness-stewart.TP.thickness-stewart.Leg.sphere.bottom-stewart.Leg.sphere.top-stewart.SP.thickness.bottom-stewart.SP.thickness.top)*0.001; % TODO

        radius_b = stewart.BP.leg.rad*0.001; % rayon emplacement support base
        radius_t = stewart.TP.leg.rad*0.001; % top radius in meters

        for i = 1:3
          % base points
          angle_m_b = (2*pi/3)* (i-1) - alpha_b;
          angle_p_b = (2*pi/3)* (i-1) + alpha_b;
          stewart.pos_base(2*i-1,:) =  [radius_b*cos(angle_m_b), radius_b*sin(angle_m_b), 0.0];
          stewart.pos_base(2*i,:) = [radius_b*cos(angle_p_b), radius_b*sin(angle_p_b), 0.0];

          % top points
          % Top points are 60 degrees offset
          angle_m_t = (2*pi/3)* (i-1) - alpha_t + 2*pi/6;
          angle_p_t = (2*pi/3)* (i-1) + alpha_t + 2*pi/6;
          stewart.pos_top(2*i-1,:) = [radius_t*cos(angle_m_t), radius_t*sin(angle_m_t), height];
          stewart.pos_top(2*i,:) = [radius_t*cos(angle_p_t), radius_t*sin(angle_p_t), height];
        end

        % permute pos_top points so that legs are end points of base and top points
        stewart.pos_top = [stewart.pos_top(6,:); stewart.pos_top(1:5,:)]; %6th point on top connects to 1st on bottom
        stewart.pos_top_tranform = stewart.pos_top - height*[zeros(6, 2),ones(6, 1)];

        %% leg vectors
        legs = stewart.pos_top - stewart.pos_base;
        leg_length = zeros(6, 1);
        leg_vectors = zeros(6, 3);
        for i = 1:6
          leg_length(i) = norm(legs(i,:));
          leg_vectors(i,:)  = legs(i,:) / leg_length(i);
        end

        stewart.Leg.lenght = 1000*leg_length(1)/1.5;
        stewart.Leg.shape.bot = [0 0; ...
                                stewart.Leg.rad.bottom 0; ...
                                stewart.Leg.rad.bottom stewart.Leg.lenght; ...
                                stewart.Leg.rad.top stewart.Leg.lenght; ...
                                stewart.Leg.rad.top 0.2*stewart.Leg.lenght; ...
                                0 0.2*stewart.Leg.lenght];

        %% Calculate revolute and cylindrical axes
        rev1 = zeros(6, 3);
        rev2 = zeros(6, 3);
        cyl1 = zeros(6, 3);
        for i = 1:6
          rev1(i,:) = cross(leg_vectors(i,:), [0 0 1]);
          rev1(i,:) = rev1(i,:) / norm(rev1(i,:));

          rev2(i,:) = - cross(rev1(i,:), leg_vectors(i,:));
          rev2(i,:) = rev2(i,:) / norm(rev2(i,:));

          cyl1(i,:) = leg_vectors(i,:);
        end


        %% Coordinate systems
        stewart.lower_leg = struct('rotation', eye(3));
        stewart.upper_leg = struct('rotation', eye(3));

        for i = 1:6
          stewart.lower_leg(i).rotation = [rev1(i,:)', rev2(i,:)', cyl1(i,:)'];
          stewart.upper_leg(i).rotation = [rev1(i,:)', rev2(i,:)', cyl1(i,:)'];
        end

        %% Position Matrix
        stewart.M_pos_base = stewart.pos_base + (height+(stewart.TP.thickness+stewart.Leg.sphere.top+stewart.SP.thickness.top+stewart.jacobian)*1e-3)*[zeros(6, 2),ones(6, 1)];

        %% Compute Jacobian Matrix
        aa = stewart.pos_top_tranform + (stewart.jacobian - stewart.TP.thickness - stewart.SP.height.top)*1e-3*[zeros(6, 2),ones(6, 1)];
        stewart.J  = getJacobianMatrix(leg_vectors', aa');
    end

    function J  = getJacobianMatrix(RM,M_pos_base)
        % RM: [3x6] unit vector of each leg in the fixed frame
        % M_pos_base: [3x6] vector of the leg connection at the top platform location in the fixed frame
        J = zeros(6);
        J(:, 1:3) = RM';
        J(:, 4:6) = cross(M_pos_base, RM)';
    end
end

function [cedrat] = initializeCedratPiezo(opts_param)
  %% Default values for opts
  opts = struct();

  %% Populate opts with input parameters
  if exist('opts_param','var')
      for opt = fieldnames(opts_param)'
          opts.(opt{1}) = opts_param.(opt{1});
      end
  end

  %% Stewart Object
  cedrat = struct();
  cedrat.k = 10e7; % Linear Stiffness of each "blade" [N/m]
  cedrat.ka = 10e7; % Linear Stiffness of the stack [N/m]

  cedrat.c = 0.1*sqrt(1*cedrat.k); % [N/(m/s)]
  cedrat.ca = 0.1*sqrt(1*cedrat.ka); % [N/(m/s)]

  cedrat.L = 80; % Total Width of the Actuator[mm]
  cedrat.H = 45; % Total Height of the Actuator [mm]
  cedrat.L2 = sqrt((cedrat.L/2)^2 + (cedrat.H/2)^2); % Length of the elipsoidal sections [mm]
  cedrat.alpha = 180/pi*atan2(cedrat.L/2, cedrat.H/2); % [deg]

  %% Save
  save('./mat/stages.mat', 'cedrat', '-append');
end

function [sample] = initializeSample(opts_param)
    %% Default values for opts
    sample = struct('radius', 100, ...
                    'height', 300, ...
                    'mass',   50,  ...
                    'offset', 0,   ...
                    'color',  [0.45, 0.45, 0.45] ...
    );

    %% Populate opts with input parameters
    if exist('opts_param','var')
        for opt = fieldnames(opts_param)'
            sample.(opt{1}) = opts_param.(opt{1});
        end
    end

    %%
    sample.k.x = 1e8;
    sample.c.x = 0.1*sqrt(sample.k.x*sample.mass);

    sample.k.y = 1e8;
    sample.c.y = 0.1*sqrt(sample.k.y*sample.mass);

    sample.k.z = 1e8;
    sample.c.z = 0.1*sqrt(sample.k.z*sample.mass);

    %% Save
    save('./mat/stages.mat', 'sample', '-append');
end