function [Drive_model,Bearing_model,LockNut_model]=Pohon_Z(pozice,deleni)
 %% VSTUPY

 
 % KULIČKOVÝ ŠROUB
Jm=0.00154;     % [kg*m2] Moment setrvačnosti motoru
s=0.04;         % [m/ot] % Stoupani
d=0.063;        % [m] % Nominalni prumer sroubu
dW=0.010319;    % [m] % Prumer kulicek v matici
alpha=45;       % [deg] contact angle of the screw balls
mL=2000;        % [kg] Hmotnost pohyblivych casti
ro=7850;        % [kg/m3] density of the screw material        
E=2.1e+011;     % [Pa] Young's modulus of the screw material        
mi=0.3;         % [-] Poisson's constant of the screw material      
L=6.7;          % [m] Celkova delka sroubu

% KULIČKOVÁ MATICE
kM=1320000000;  % [N/m] Axiální tuhost kulickove matice
bM=100000*50;   % [Ns/m] Axialni tlumeni kulickove matice

%LOŽISKO
kAL=4900000000; % [N/m] Axiální tuhost loziska, ZARN55115-TV
bAL=100000*50;  % [Ns/m] Axiální tlumení ložiska

% ŘEMENOVÝ PŘEVOD
dw1=0.05602;    % [m] Prumer remenice na motoru (mala)
dw2=0.22918;    % [m] Prumer remenice na sroubu (velka)
p_rp=4.1;       % [-] Prevod remenoveho prevodu    
JR1=0.000161+0.0000123;   % [kg.m^2] Moment setrvacnosti (Mala remenice + pouzdro)
JR2=0.03755+0.0014;       % [kg.m^2] Moment setrvacnosti (Velka remenice + pouzdro)
krm1=562000;    % [m.N/m] Tuhost 1m řemenu  
a=0.1626;       % [m] Osova vzdalenost remenic

%% POMOCNÉ VÝPOČTY
% Kulickovy sroub
G=E/(2*(1+mi));            % [Pa] Modul pružnosti ve smyku
D=d-dW*cosd(alpha);        % [m] Efektivni prumer sroubu
A=(pi*D*D)/4;              % [m2] Efektivni plocha sroubu
J=pi*D^4/32;               % [m4] Efektivni polární kvadratický moment průřezu kuličkového šroubu
h=s/(2*pi);                % [m/rad] Prevod na SI jednotku (stoupani sroubu)

% Remenovy prevod
Ls1=sqrt(a^2-(dw2/2-dw1/2)^2);   % [m] Delka jedne vetve remene
kR=2*krm1/Ls1*dw1^2/4;           % [N.m/rad] Tuhost remene redukovana na vstupni souradnici

%% SUBMODELY
% Motor
    Mm=Jm;

% Remenovy prevod
    MR=[JR1 0; 0 JR2];
    KR=[kR -kR*p_rp;-kR*p_rp kR*p_rp*p_rp];

%Lozisko
    % in:  u_frame  u_screw     v_frame     v_screw
    % out: F_frame  F_screw
KB=[-kAL  kAL -bAL  bAL
     kAL -kAL  bAL -bAL];
    
% Elementy KŠ
% meni se delka elemnetu 
switch deleni
    case 1
        % 2 elementy před maticí
            Me1ax=ro*A*(pozice/2)/6*[2 1;1 2];
            Ke1ax=E*A/(pozice/2)*[1 -1;-1 1];
            Me1tor=ro*J*(pozice/2)/6*[2 1;1 2];
            Ke1tor=G*J/(pozice/2)*[1 -1;-1 1];
        % 6 elementů za maticí
            Me2ax=ro*A*((L-pozice)/6)/6*[2 1;1 2];
            Ke2ax=E*A/((L-pozice)/6)*[1 -1;-1 1];
            Me2tor=ro*J*((L-pozice)/6)/6*[2 1;1 2];
            Ke2tor=G*J/((L-pozice)/6)*[1 -1;-1 1];

    case 2        
        % 4 elementy před matici
            Me1ax=ro*A*(pozice/4)/6*[2 1;1 2];
            Ke1ax=E*A/(pozice/4)*[1 -1;-1 1];
            Me1tor=ro*J*(pozice/4)/6*[2 1;1 2];
            Ke1tor=G*J/(pozice/4)*[1 -1;-1 1];
        % 4 elementy za matici
            Me2ax=ro*A*((L-pozice)/4)/6*[2 1;1 2];
            Ke2ax=E*A/((L-pozice)/4)*[1 -1;-1 1];
            Me2tor=ro*J*((L-pozice)/4)/6*[2 1;1 2];
            Ke2tor=G*J/((L-pozice)/4)*[1 -1;-1 1];
    case 3
        % 6 elementů před maticí
            Me1ax=ro*A*(pozice/6)/6*[2 1;1 2];
            Ke1ax=E*A/(pozice/6)*[1 -1;-1 1];
            Me1tor=ro*J*(pozice/6)/6*[2 1;1 2];
            Ke1tor=G*J/(pozice/6)*[1 -1;-1 1];
        % 2 elementy za maticí
            Me2ax=ro*A*((L-pozice)/2)/6*[2 1;1 2];
            Ke2ax=E*A/((L-pozice)/2)*[1 -1;-1 1];
            Me2tor=ro*J*((L-pozice)/2)/6*[2 1;1 2];
            Ke2tor=G*J/((L-pozice)/2)*[1 -1;-1 1];
end;

% Kulickova matice
    % in:   u_frame     u_screw     fi_screw    v_frame     v_screw     om_screw
    % out:  F_frame     F_screw     T_screw
KKM=[-kM       kM       kM*h     -bM       bM       bM*h
      kM      -kM      -kM*h      bM      -bM      -bM*h
      kM*h    -kM*h    -kM*(h^2)  bM*h    -bM*h    -bM*(h^2)];
 
 
 %% Tvorba celkove matice K a M 
 %JEDNOSTRANNÝ NÁHON
 %coordinate indexes
    fiM1=1; % 1st global coordinate - 1st (left) motor rotation
    fi1=2;  % 2nd global coordinate - screw rotation - motor 1 end
    fi2=3;  % 3rd global coordinate - screw rotation
    fi3=4;  % 4th global coordinate - screw rotation - middle of the first section
    fi4=5;  % 5th global coordinate - screw rotation
    fi5=6;  % 6th global coordinate - screw rotation - under the lock nut of the ball screw
    fi6=7;  % 7th global coordinate - screw rotation
    fi7=8;  % 8th global coordinate - screw rotation - middle of the second section
    fi8=9;  % 9th global coordinate - screw rotation
    fi9=10; % 10th global coordinate - screw rotation - motor 2 end
    x1=11;  % 11th global coordinate - screw translation - motor 1 end with axial bearing
    x2=12;  % 12th global coordinate - screw translation
    x3=13;  % 13th global coordinate - screw translation - middle of the first section
    x4=14;  % 14th global coordinate - screw translation
    x5=15;  % 15th global coordinate - screw translation - under the lock nut of the ball screw
    x6=16;	% 16th global coordinate - screw translation
    x7=17;  % 17th global coordinate - screw translation - middle of the second section
    x8=18;  % 18th global coordinate - screw translation
    x9=19;  % 19th global coordinate - screw translation
    nCoord=19; % number of coordinates
     
%Vysledne matice
    M=zeros(nCoord);    % iniciace
    K=M;                % iniciace
    
%JEDNOSTRANNÝ NÁHON
    % 1 motor
    M=Add_submatrix(M,Mm,fiM1);
    % Remenovy prevod pro 1. motor (LEVÝ)
    M=Add_submatrix(M,MR,[fiM1, fi1]);    
    K=Add_submatrix(K,KR,[fiM1, fi1]);
    
switch deleni
    case 1
        %TORZNI ELEMENTY
        M=Add_submatrix(M,Me1tor,[fi1, fi2]);  K=Add_submatrix(K,Ke1tor,[fi1, fi2]); % 1st torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi2, fi3]);  K=Add_submatrix(K,Ke1tor,[fi2, fi3]); % 2nd torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi3, fi4]);  K=Add_submatrix(K,Ke2tor,[fi3, fi4]); % 3rd torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi4, fi5]);  K=Add_submatrix(K,Ke2tor,[fi4, fi5]); % 4th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi5, fi6]);  K=Add_submatrix(K,Ke2tor,[fi5, fi6]); % 5th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi6, fi7]);  K=Add_submatrix(K,Ke2tor,[fi6, fi7]); % 6th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi7, fi8]);  K=Add_submatrix(K,Ke2tor,[fi7, fi8]); % 7th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi8, fi9]);  K=Add_submatrix(K,Ke2tor,[fi8, fi9]); % 8th torsional element of the screw

        %AXIALNI ELEMENTY
        M=Add_submatrix(M,Me1ax,[x1, x2]);  K=Add_submatrix(K,Ke1ax,[x1, x2]); % 1st axial element of the screw
        M=Add_submatrix(M,Me1ax,[x2, x3]);  K=Add_submatrix(K,Ke1ax,[x2, x3]); % 2nd axial element of the screw
        M=Add_submatrix(M,Me2ax,[x3, x4]);  K=Add_submatrix(K,Ke2ax,[x3, x4]); % 3rd axial element of the screw
        M=Add_submatrix(M,Me2ax,[x4, x5]);  K=Add_submatrix(K,Ke2ax,[x4, x5]); % 4th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x5, x6]);  K=Add_submatrix(K,Ke2ax,[x5, x6]); % 5th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x6, x7]);  K=Add_submatrix(K,Ke2ax,[x6, x7]); % 6th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x7, x8]);  K=Add_submatrix(K,Ke2ax,[x7, x8]); % 7th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x8, x9]);  K=Add_submatrix(K,Ke2ax,[x8, x9]); % 8th axial element of the screw
    case 2
        %TORZNI ELEMENTY
        M=Add_submatrix(M,Me1tor,[fi1, fi2]);  K=Add_submatrix(K,Ke1tor,[fi1, fi2]); % 1st torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi2, fi3]);  K=Add_submatrix(K,Ke1tor,[fi2, fi3]); % 2nd torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi3, fi4]);  K=Add_submatrix(K,Ke1tor,[fi3, fi4]); % 3rd torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi4, fi5]);  K=Add_submatrix(K,Ke1tor,[fi4, fi5]); % 4th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi5, fi6]);  K=Add_submatrix(K,Ke2tor,[fi5, fi6]); % 5th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi6, fi7]);  K=Add_submatrix(K,Ke2tor,[fi6, fi7]); % 6th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi7, fi8]);  K=Add_submatrix(K,Ke2tor,[fi7, fi8]); % 7th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi8, fi9]);  K=Add_submatrix(K,Ke2tor,[fi8, fi9]); % 8th torsional element of the screw

        %AXIALNI ELEMENTY
        M=Add_submatrix(M,Me1ax,[x1, x2]);  K=Add_submatrix(K,Ke1ax,[x1, x2]); % 1st axial element of the screw
        M=Add_submatrix(M,Me1ax,[x2, x3]);  K=Add_submatrix(K,Ke1ax,[x2, x3]); % 2nd axial element of the screw
        M=Add_submatrix(M,Me1ax,[x3, x4]);  K=Add_submatrix(K,Ke1ax,[x3, x4]); % 3rd axial element of the screw
        M=Add_submatrix(M,Me1ax,[x4, x5]);  K=Add_submatrix(K,Ke1ax,[x4, x5]); % 4th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x5, x6]);  K=Add_submatrix(K,Ke2ax,[x5, x6]); % 5th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x6, x7]);  K=Add_submatrix(K,Ke2ax,[x6, x7]); % 6th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x7, x8]);  K=Add_submatrix(K,Ke2ax,[x7, x8]); % 7th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x8, x9]);  K=Add_submatrix(K,Ke2ax,[x8, x9]); % 8th axial element of the screw
    case 3
        %TORZNI ELEMENTY
        M=Add_submatrix(M,Me1tor,[fi1, fi2]);  K=Add_submatrix(K,Ke1tor,[fi1, fi2]); % 1st torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi2, fi3]);  K=Add_submatrix(K,Ke1tor,[fi2, fi3]); % 2nd torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi3, fi4]);  K=Add_submatrix(K,Ke1tor,[fi3, fi4]); % 3rd torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi4, fi5]);  K=Add_submatrix(K,Ke1tor,[fi4, fi5]); % 4th torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi5, fi6]);  K=Add_submatrix(K,Ke1tor,[fi5, fi6]); % 5th torsional element of the screw
        M=Add_submatrix(M,Me1tor,[fi6, fi7]);  K=Add_submatrix(K,Ke1tor,[fi6, fi7]); % 6th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi7, fi8]);  K=Add_submatrix(K,Ke2tor,[fi7, fi8]); % 7th torsional element of the screw
        M=Add_submatrix(M,Me2tor,[fi8, fi9]);  K=Add_submatrix(K,Ke2tor,[fi8, fi9]); % 8th torsional element of the screw

        %AXIALNI ELEMENTY
        M=Add_submatrix(M,Me1ax,[x1, x2]);  K=Add_submatrix(K,Ke1ax,[x1, x2]); % 1st axial element of the screw
        M=Add_submatrix(M,Me1ax,[x2, x3]);  K=Add_submatrix(K,Ke1ax,[x2, x3]); % 2nd axial element of the screw
        M=Add_submatrix(M,Me1ax,[x3, x4]);  K=Add_submatrix(K,Ke1ax,[x3, x4]); % 3rd axial element of the screw
        M=Add_submatrix(M,Me1ax,[x4, x5]);  K=Add_submatrix(K,Ke1ax,[x4, x5]); % 4th axial element of the screw
        M=Add_submatrix(M,Me1ax,[x5, x6]);  K=Add_submatrix(K,Ke1ax,[x5, x6]); % 5th axial element of the screw
        M=Add_submatrix(M,Me1ax,[x6, x7]);  K=Add_submatrix(K,Ke1ax,[x6, x7]); % 6th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x7, x8]);  K=Add_submatrix(K,Ke2ax,[x7, x8]); % 7th axial element of the screw
        M=Add_submatrix(M,Me2ax,[x8, x9]);  K=Add_submatrix(K,Ke2ax,[x8, x9]); % 8th axial element of the screw
end;    
    
    %% VÝSTUP - STATE-SPACE MODEL
    [V,Lam] = eig(K,M);             % modal transformation, matrix V of eigen vektors is normed to unit mass
    [L,I]=sort(diag(Lam));          % sorting of eigen values
    L=round(L*1e6)*1e-6;            % rounding of eigen values
    Lam=diag(L);                    % back to matrix form
    V=V(:,I);                       % sorting of eigen vectors according to eigen values
    
    dzeta=0.1;                      % default value of the modal damping
    Dzeta=dzeta*diag(double(L~=0)); % matrix form of modal damping, rigid modes with zero damping
    
    om=sqrt(L);                     % eigen frequencies
    Cq=2*Dzeta*diag(om);            % matrix of damping
    
    % state-space
    A=[zeros(size(Lam)) eye(size(Lam));
           -Lam            -Cq];
    B=[zeros(size(V'));V'];
    C=[    V             zeros(size(V))
              zeros(size(V))         V];
    D=[zeros(size(V));zeros(size(V))];
    
    
    switch deleni
        case 1
            in=[fiM1 fi3 x1 x3];    
            out=[fiM1 fi3 x1 x3,... 
                fiM1+nCoord fi3+nCoord x1+nCoord x3+nCoord];    
        case 2
            in=[fiM1 fi5 x1 x5];    
            out=[fiM1 fi5 x1 x5,... 
                fiM1+nCoord fi5+nCoord x1+nCoord x5+nCoord];    
        case 3
            in=[fiM1 fi7 x1 x7];    
            out=[fiM1 fi7 x1 x7,... 
                fiM1+nCoord fi7+nCoord x1+nCoord x7+nCoord];   
    end;
    
    B=B(:,in);
    C=C(out,:);
    D=D(out,in);                    
    
    Drive_model=ss(A,B,C,D);        % continuous state-space model
    Bearing_model=ss(KB);
    LockNut_model=ss(KKM);
end

%% Přiřazení submatic do celkové matice

function out=Add_submatrix(in,sub,ind)
    out=in;
    out(ind,ind)=out(ind,ind)+sub;
end 
   