function u = MAX_stand_input_timemesh(t,y,alpha_conv,alpha_Boden)
%% compute the inputs u for the simulation of MAX stand

% Copyright 2023 Collaborative Research Centre Transregio 96 (CRC/TR 96)
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
% TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER
% OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%% MAX stand: description

% Due to symmetrie reasons it is sufficient to consider one stand of the 
% MAX frame. This model consists of 4 subassemblies (SA), which are assemblies
% merging several components that are not moves relative to each other.

% There are 3 kinds of contact conditions:
% IKB: Neumann boundary condition: contact condition between components of 
%      the same subassembly with q=alpha_IKB*(T1-T2) and alpha_IKB=10
% KB:  Neumann boundary condition: contact condition between two different
%      subassemblies with q=alpha_KB*(T1-T2) and alpha_KB=1000
% convection boundary conditions as well as the contact to the floor 
% modeled by Robin boundary conditions with q=alpha(T_ext-T)
% convection: alpha_conv=10, T_ext according to sensor_ambient_temperatur.m 
% (current ambient temperature - initial temperature of 20°C) (the model was 
% shifted by the initial temperature to guarantee a zero inital value
% for model order reduction)
% contact to the floor: alpha_Boden=100, T_ext=20°C
% For the Neumann boundary conditions two contact areas are needed:

% SA1: base (11 volume bodies, 25872 nodes)      
% SA2: z-stand with guide rail and bottom plate (22 volume bodies, 39527 nodes)
% SA3: motor with belt housing, bearing seat and outer bearing ring 
%      (20 volume bodies, 13551 nodes)
% SA4: spindle with inner bearing ring (6 volume bodies, 4813 nodes)

%% computation of the inputs u(t) for every subassembly

% load constants
c_pool = MAX_constants();

nsa = 4;                                  % number of SA
T_ext = mean(sensor_ambient_temperature(t));  % ambient temperature
d = axis_distance(t);                     % position z-slide on spindle
v = axis_velocity(t);                     % feed speed z-slide
alpha_IKB = 10;                           % heat transfer coefficient for IKB
alpha_KB = 1000;                          % heat transfer coefficient for KB
alpha_zahn = 100;                         % heat transfer coefficient for motor
                                          % gearwheel and toothed belt
alpha_lager = mean(alpha_bearing(v,c_pool));  % heat transfer coefficient for bearing
T_Boden = 0;                              % floor temperature (T=20°C 
                                          % -> shift -20°C -> T_Boden=0°C)                                                                                                                    
u = cell(1,nsa);
u{1} = zeros(10,1);
u{2} = zeros(30,1);
u{3} = zeros(6,1);
u{4} = zeros(23,1);

%% SA1:
u{1}(1) = alpha_IKB*(y{1}(2)-y{1}(1));        % IKB
u{1}(2) = alpha_IKB*(y{1}(1)-y{1}(2));        % IKB
u{1}(3) = alpha_IKB*(y{1}(4)-y{1}(3));        % IKB
u{1}(4) = alpha_IKB*(y{1}(3)-y{1}(4));        % IKB
u{1}(5) = alpha_IKB*(y{1}(6)-y{1}(5));        % IKB
u{1}(6) = alpha_IKB*(y{1}(5)-y{1}(6));        % IKB
u{1}(7) = alpha_KB*(y{3}(1)-y{1}(7));         % KB1 contact between SA1/u7 and SA3/u1
u{1}(8) = alpha_KB*(y{2}(9)-y{1}(8));         % KB2 contact between SA1/u8 and SA2/u9 
u{1}(9) = alpha_Boden*T_Boden;                % contact with floor
u{1}(10) = alpha_conv*T_ext;                  % convection with ambience

%% SA2:
w = 2*pi/c_pool.sp_p*abs(v);                % angular velocity
Schienen_Anteil = 0.5;

% starting position for the segments of the guide rail
character_difference_Schiene = ...
  [0;...
   0.039818503389795;...
   0.079596269812030;...
   0.119662881766638;...
   0.159553683764328;...
   0.199250223909775;...
   0.238988954254294;...
   0.278946299445552;...
   0.318794578421047;...
   0.358538954366228;...
   0.398372263599169;...
   0.438393416773297;...
   0.478205665526311;...
   0.517989935695338;...
   0.558021890596418;...
   0.597678134727693;...
   0.637686569958904;...
   0.677358944765107;...
   0.717404884338751;...
   0.757173436911123];

% length of one segment of the guide rail
character_length_Schiene     = 0.11955;

character_distance_Schiene = abs(kron(character_difference_Schiene,ones(1,length(t)))-kron(ones(length(character_difference_Schiene),1),d));
overlap_Schiene            = character_distance_Schiene<character_length_Schiene/2;

zeller_Schiene = sum(overlap_Schiene);

% friction loss for the segments of the guide rails
f_fric_Schiene = kron(ones(length(character_difference_Schiene),1),(c_pool.mu_factor1*Schienen_Anteil*c_pool.mu_M)*w./(c_pool.sp_ASegment*zeller_Schiene)).*overlap_Schiene;
%--------------------------------------------------------------------------

u{2}(1) = alpha_IKB*(y{2}(2)-y{2}(1));        % IKB
u{2}(2) = alpha_IKB*(y{2}(1)-y{2}(2));        % IKB
u{2}(3) = alpha_IKB*(y{2}(4)-y{2}(3));        % IKB
u{2}(4) = alpha_IKB*(y{2}(3)-y{2}(4));        % IKB
u{2}(5) = alpha_IKB*(y{2}(6)-y{2}(5));        % IKB
u{2}(6) = alpha_IKB*(y{2}(5)-y{2}(6));        % IKB
u{2}(7) = alpha_IKB*(y{2}(8)-y{2}(7));        % IKB
u{2}(8) = alpha_IKB*(y{2}(7)-y{2}(8));        % IKB
u{2}(9) = alpha_KB*(y{1}(8)-y{2}(9));         % KB2 contact between SA1/u8 and SA2/u9
u{2}(10) = alpha_conv*T_ext;                  % convection
u{2}(11:30) = mean(f_fric_Schiene,2);         % 20 segments of guide rail
                                              % contact with z-slide (only
                                              % friction)
                                            
%% SA3:
W_sensor = 0.402378257201514; % for angular velocity w = 1;
                              % motor power loss (=W_sensor*w)
Aussenring_Anteil = 0.5;
Innenring_Anteil = 0.5;
M_0  = (4501*c_pool.l_f_0*c_pool.l_d_m^3)*(c_pool.l_ny_40*w).^(2/3);
M    = M_0 + c_pool.l_M_1 + c_pool.l_M_d;
% 
f_fric_aussen = mean((Aussenring_Anteil/c_pool.l_A)*M.*w);
f_fric_innen  = mean((Innenring_Anteil/c_pool.l_A)*M.*w);

%--------------------------------------------------------------------------

u{3}(1) = alpha_KB*(y{1}(7)-y{3}(1));         % KB1 contact between SA1/u7 and SA3/u1
u{3}(2) = alpha_conv*T_ext;                   % convection
u{3}(3) = mean((c_pool.Anteil_rotor/c_pool.mo_ARotor)*W_sensor.*w);   % motor: rotor
u{3}(4) = mean((c_pool.Anteil_stator/c_pool.mo_AStator)*W_sensor.*w); % motor: stator
u{3}(5) = alpha_zahn*(y{4}(23)-y{3}(5));      % motor: gearwheel
                                              % KB contact between SA3/u5 and SA4/u23
u{3}(6) = f_fric_aussen+alpha_lager*(y{4}(2)-y{3}(6));  % friction loss outer bearing ring
                                            % KB contact between SA3/u6 and SA4/u2
                                            
%% SA4:
Spindel_Anteil = 0.5*(1/4);

% starting position for the segments of the spindle
character_difference_Spindel = ...
  [0;...
   0.037219042222223;...
   0.073600000972222;...
   0.110399999861112;...
   0.147200000000001;...
   0.184000001111111;...
   0.220800000972222;...
   0.257599999861111;...
   0.294400000000000;...
   0.331200001111112;...
   0.368000000972224;...
   0.404799999861111;...
   0.441600000000000;...
   0.478400001111112;...
   0.515200000972222;...
   0.551999999861110;...
   0.589429351250001;...
   0.626229351250001;...
   0.663029352222222;...
   0.699791259861111];

% length of one segment of the spindle
character_length_Spindel     = 0.1472;

character_distance_Spindel = abs(kron(character_difference_Spindel,ones(1,length(t)))-kron(ones(length(character_difference_Spindel),1),d));%character_distance = abs(character_difference-d);
overlap_Spindel            = character_distance_Spindel<character_length_Spindel/2;

zeller_Spindel = sum(overlap_Spindel);

% friction loss for the segments of the spindle
f_fric_Spindel = (c_pool.mu_factor1*Spindel_Anteil*c_pool.mu_M)*w./(c_pool.sp_ASegment*(zeller_Spindel+eps)).*overlap_Spindel;% Reibverluste

%--------------------------------------------------------------------------

u{4}(1) = alpha_conv*T_ext;                   % convection
u{4}(2) = f_fric_innen+alpha_lager*(y{3}(6)-y{4}(2)); % friction loss inner bearing ring
                                              % KB contact between SA3/u6 and SA4/u2                                                                                        
u{4}(3:22) = mean(f_fric_Spindel,2);          % 20 spindle segments
                                              % contact with z-Schlitten
                                              % (only friction)

u{4}(23) = alpha_zahn*(y{3}(5)-y{4}(23));     % toothed belt
                                              % KB contact between SA3/u5 and SA4/u23 

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% local functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% compute heat transfer coefficient for the bearing depending on velocity v

function alpha = alpha_bearing(v,c_pool)

%--------------------- relevant constants ----------------------
%W_Grund = 5000;

%eta_k = 220e-6;    % [m^2/s] kinematic viscosity
%rho   = 900;       % [kg/m^3] density of fluid
%eta   = eta_k*rho; % [Ns/m^2] dynamic viscosity
%s     = 10e-6;     % [m] gap

n_min = 10;
w_min = n_min/60*2*pi;
%-----------------------------------------------------------------

w = 2*pi/c_pool.sp_p*abs(v);

w_cor = w_min + (w>w_min).*(w-w_min);
L = (c_pool.l_i_wk*590*c_pool.l_d_wk^2)*...
    sqrt(log(w_cor*(c_pool.l_d_i+c_pool.l_d_wk)/2000)-log(c_pool.l_d_wk)+7);
alpha = (2/c_pool.l_A)*L;

end

%%%%%%%%%%%%%%%%%%%%%%%%

%% compute the feed speed (velocity) v of the z-slide at time t

function v = axis_velocity(t)

% velocity of the z-slide for every hour of the process to be simulated
% (overall 16 hours)

modulo = 3600; % [s] (=1h)
velocity = [0.0, 2.5, 2.5, 2.5, 5.0, 5.0, 5.0, 20.0, 20.0, 20.0, 50.0, ...
            50.0, 50.0, 0.0, 0.0, 0.0]/60; % [m/s]
max_index = length(velocity);
index = fix(t/modulo)+1;
index = max_index + (index<max_index).*(index-max_index);

v = velocity(index);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%

%% compute the position d of the z-slide on the spindle and guide rail at
%  time t

function d = axis_distance(t)

modulo = 3600;  %[s] (=1h)
L = 0.75715;    %length of the axis [m] 
% velocity, start and stop position of the z-slide for every hour of the 
% process to be simulated (overall 16 hours)
velocity  = [0.0, 2.5, 2.5, 2.5, 5.0, 5.0, 5.0, 20.0, 20.0, 20.0, 50.0,...
             50.0, 50.0, 0.0, 0.0, 0.0]/60; % [m/s]
start     = [0, 0, 0, L/2, 0, 0, L/2, 0, 0, L/2, 0, 0, L/2, 0, 0, 0];
stop      = [0, L/2, L, L/2, L/2, L, L/2, L/2, L, L/2, L/2, L, L/2, 0, 0, 0];
max_index = length(start);
position  = zeros(1,max_index);
pos = 0;

for i = 2:13
    v = velocity(i);
    rest = abs(start(i)-pos);
    
    %--------------
    time = 3600 - rest/v;
    
    s = time*v;
    richtung = mod(fix(s/stop(i)),2);
    distance = mod(s,stop(i));
    pos = start(i) + richtung*stop(i) + distance*(1 - 2*richtung);
    position(i+1) = pos;
end

rest_time = mod(t,modulo);
index = fix(t/modulo)+1;
index = max_index + (index<max_index).*(index-max_index);

v = velocity(index);

s = rest_time.*v - abs(start(index)-position(index));

richtung = mod(fix(s./(stop(index)-eps)),2);
distance = mod(s,stop(index));
pos = start(index) + distance.*(1 - 2*richtung) + richtung.*stop(index);

d = (s<0).*(position(index) + round((position(index)<start(index))-1/2).*(rest_time.*v))...
    + (s>=0).*pos;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%
%% compute ambient temperature T at time t

function T = sensor_ambient_temperature(t)

%--------------------- relevant constants ----------------------
modulo = 3600; % s (=1h)

% the system is shifted by the initial temperature of 20°C to guarantee a
% zero initial condition for model order reduction
T_air = [3.405, 1.719, 2.228, 3.402, 3.764, 4.499, 4.876, 4.688, 3.935,...
         3.049, 2.43, 2.07, 1.832, 1.451, 1.652, 1.53, 1.185, 1.048,... 
         1.096, 0.497]; %[°C]
max_index = length(T_air);
%-----------------------------------------------------------------

index = fix(t/modulo)+1;
index = max_index + (index<max_index).*(index-max_index);

T = T_air(index);

end
