function [u,kappa_ext] = MiniHex_input(t,y)

% Compute input for the four subassemblies (SA)

% 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.


%% MiniHex: Explanation of the contact areas (CA) of the subassemblies (SA)

% SA1: belt housing with motor
       % CA1 - load area of the stator                        (heat flux)
       % CA2 - load area of the rotor                         (heat flux)
       % CA3 - friction area: outer ring of the bearing       (heat flux)
       %       (contact between SA1 and SA2)
       % CA4 - outer surface: upper part of the motor         (convection)
       % CA5 - outer surface: middle part of the motor        (convection)
       % CA6 - outer surface: lower part of the motor         (convection)
       % CA7 - outer surface: (convection)
       % CA8 - outer surface: (convection)
       % CA9 - outer surface: (convection)

% SA2: spindle
       % CA1 - contact area between motor and spindle (heat flux)
       %       (heat transfered from motor to spindle by the belt)
       % CA2 - friction area: inner ring of the bearing       (heat flux)
       %       (contact area between SA2 and SA1)
       % CA3:CA30 - contact areas: spindle segments           (heat flux)
       %            (contact areas between SA2 and SA3 
       %            as well as SA2 and SA4, depends on position of the nut)
       % CA31 - end part of spindle (to SA1)                  (convection)
       % CA32 - end part of spindle (to SA4)                  (convection)
       
% SA3: spindle nut
       % CA1 - friction area: spindle nut                     (heat flux)
       %       (contact between SA3 and SA2)
       % CA2 - contact area between spindle nut and feed tube (heat flux)
       %       (contact between SA3 and SA4)
       % CA3 - outer surface: spindle nut                     (convection)

% SA4: feed tube
       % CA1 - contact area between feed tube and spindle nut (heat flux)
       %       (contact between SA4 and SA3)
       % CA2:CA8 - contact area between inner feed tube 
       %           segments and spindle segments              (heat flux)
       %           (contact between SA4 and SA2)
       % CA9 - outer surface: beginning of tube               (convection)
       % CA10 - outer surface: end part of tube with joint    (convection)
       % CA11:CA17 - outer feed tube segments                 (convection)
       %             (contact with the ambience)

%% Compute the thermally dependent heat transfer coefficients kappa_ext and
%% the inputs u(t) for the contact areas of each subassembly       

% sensor1: ambient temperature
% sensor2: heat induced by the motor (motor power loss)

[c_pool, VARIANTE1, Messreihe1, Messreihe2,...
    Messreihe1_zyklisch, T_Grund] =  MiniHex_constants();

load('character_difference.mat','character_difference')

length_nut=57/1000;
length_tube=75/1000;
% sensor 1 (ambient temperature)
T_ext = sensor_ambientTemperature(t,T_Grund,Messreihe2);
% sensor 2 (motor power loss)
sensor2 = sensor_motorPowerLoss(t,c_pool,VARIANTE1,Messreihe1,Messreihe1_zyklisch);
% compute velocity and position
[v, d] = achse1_velocity_distance(t,VARIANTE1,Messreihe1,Messreihe1_zyklisch);

kappa_ext=cell(1,4); % temperature dependent heat transfer coefficient with ambience
u=cell(1,4); % input
u{1}=zeros(9,1);
u{2}=zeros(32,1);
u{3}=zeros(3,1);
u{4}=zeros(17,1);

% contact 1 (SA1/CA1 and sensor2)
Anteil_stator = 0.5;
u{1}(1)=(Anteil_stator/c_pool.mo_AStartor)*sensor2;

% contact 2 (SA1/CA2 and sensor2)
Anteil_rotor = 0.0; 
u{1}(2)=(Anteil_rotor/c_pool.mo_ARotor)*sensor2;

%%

w = 2*pi/c_pool.sp_p*abs(v);
n_min = 10;
w_min = n_min/60*2*pi;
w_cor = w_min + (w>w_min).*(w-w_min);

% contact 3 (SA1/CA3 und SA2/CA2)
L1 = (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);
alphaKB3=(2/c_pool.l_A)*L1;
Aussenring_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;
frictionKB3=(Aussenring_Anteil/c_pool.l_A)*M.*w;
u{1}(3)=alphaKB3*(y{2}(2)-y{1}(3))+frictionKB3;
u{2}(2)=alphaKB3*(y{1}(3)-y{2}(2))+frictionKB3;

% contact 4 (SA1/CA4:9 and sensor1)
delta_T1(1:6,1) = abs(y{1}(4:9,1)-T_ext);
alphaKB4=4 + 2.2*delta_T1.^0.33;
u{1}(4:9)=alphaKB4*T_ext;

% contact 5 (SA2/CA1 and sensor2)
Anteil_Riehmenleitung = 0.1;
u{2}(1)=(Anteil_Riehmenleitung/c_pool.sp_ARiehmenaufnahme)*sensor2;

% contact 6 (SA3/CA1 und SA2/CA3:30) 
character_distance_spindel = abs(kron(character_difference{6},ones(1,length(t)))-kron(ones(length(character_difference{6}),1),d));%character_distance = abs(character_difference-d);
overlap_spindel            = character_distance_spindel<length_nut/2;

L2 = (c_pool.mu_i_wk*590*c_pool.mu_d_wk^2)*sqrt(log(w_cor*(c_pool.sp_d+c_pool.mu_d_wk)/2000)-log(c_pool.mu_d_wk)+7);
alphaspindel=(c_pool.mu_i_u*L2./(c_pool.sp_ASegment*sum(overlap_spindel)))*overlap_spindel; %input spindle
alphamutter=(c_pool.mu_i_u/c_pool.mu_A)*L2;                                                 %input nut
Spindel_Anteil = 0.5;
Mutter_Anteil = 0.5;

for i=1:28
    if overlap_spindel(i)==1
       u{2}(2+i)=alphaspindel(i)*(y{3}(1)-y{2}(2+i))+...
                 +(c_pool.mu_factor1*Spindel_Anteil*c_pool.mu_M)*w./(c_pool.sp_ASegment*sum(overlap_spindel));               
    end
end

u{3}(1)=alphamutter*(overlap_spindel'*y{2}(3:30)/sum(overlap_spindel)-y{3}(1))+...
        + c_pool.mu_factor1*Mutter_Anteil*c_pool.mu_M/c_pool.mu_A*w;

% contact 7 (SA2/CA31:32 and sensor1)
delta_T2(1:2,1) = abs(y{2}(31:32,1)-T_ext);
alphaKB7=[(1.8+2.2*delta_T2(1,1).^0.33+3.75*(c_pool.sp_d/2)^0.4*w.^0.7) 0]';
u{2}(31:32)=alphaKB7*T_ext;

% contact 8 (SA2/CA3:30 and sensor1)
delta_T5(1:28,1) = abs(y{2}(3:30,1)-T_ext);
character_distance_pos = -(kron(character_difference{8},ones(1,length(t)))-kron(ones(length(character_difference{8}),1),d));%character_distance_pos = (character_difference-d);
overlap_frei           = (character_distance_pos>length_nut/2);
alphaKB8=(1.8 + 2.2*delta_T5.^0.33 + 3.75*(c_pool.sp_d/2)^0.4*w.^0.7);

for i=1:28
    if overlap_frei(i)==1
       u{2}(2+i)=alphaKB8(i)*(T_ext-y{2}(2+i));
    end
end

% contact 9 (SA3/CA2 und SA4/CA1)
alphaKB9=700;
u{3}(2)=alphaKB9*(y{4}(1)-y{3}(2));
u{4}(1)=alphaKB9*(y{3}(2)-y{4}(1));
 
% contact 10 (SA3/CA3 and sensor1)
delta_T3 = abs(y{3}(3)-T_ext);
alphaKB10=1.8 + 2.2*delta_T3.^0.33 + 4.7*abs(v)^0.6/c_pool.vr_d_a^0.4;
u{3}(3)=alphaKB10*T_ext;

% contacts 11-17 (SA4/CA2:8 und SA2/CA3:30)
% contact 11 (SA4/CA2 und SA2/CA3:30)
% contact 12 (SA4/CA3 und SA2/CA3:30)
% contact 13 (SA4/CA4 und SA2/CA3:30)
% contact 14 (SA4/CA5 und SA2/CA3:30)
% contact 15 (SA4/CA6 und SA2/CA3:30)
% contact 16 (SA4/CA7 und SA2/CA3:30)
% contact 17 (SA4/CA8 und SA2/CA3:30)

character_distance_rohr=cell(1,7);
overlap_rohr=cell(1,7);
alpha_rohr=(1.3 + 0.026/c_pool.sp_s + (1.24/c_pool.sp_s^0.22)*(((c_pool.sp_d+c_pool.vr_d_i)/4)*w).^0.78 + c_pool.vr_factor1);

for k=1:7
    character_distance_rohr{k} = abs(kron(character_difference{10+k},ones(1,length(t)))-kron(ones(length(character_difference{10+k}),1),d));%character_distance = abs(character_difference-d);
    overlap_rohr{k}            = character_distance_rohr{k}<length_tube/2;
    
    if sum(overlap_rohr{k}~=0)
       u{4}(1+k)=alpha_rohr*(((overlap_rohr{k})'*y{2}(3:30))/(sum(overlap_rohr{k})+eps)-y{4}(1+k));
    end
    
    for i=1:28
        if overlap_rohr{k}(i,1)==1
           u{2}(2+i)=alpha_rohr*(y{4}(1+k)-overlap_rohr{k}(i,1)*(y{2}(2+i)));  
        end
    end
end

% contact 18 (SA4/CA9:17 and sensor1)
delta_T4(1:9,1) = abs(y{4}(9:17)-T_ext);
alphaKB18=1.8 + 2.2*delta_T4.^0.33 + 4.7*abs(v)^0.6/c_pool.vr_d_a^0.4 +4;
u{4}(9:17)=alphaKB18*T_ext;

kappa_ext{1}=alphaKB4;
kappa_ext{2}=alphaKB7;
kappa_ext{3}=alphaKB10;
kappa_ext{4}=alphaKB18;     

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% local functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute velocity and position

function [v,d] = achse1_velocity_distance(t, VARIANTE1, Messreihe1, ...
                 Messreihe1_zyklisch)

if VARIANTE1

    %--------------------- constants ----------------------
    kanal_num_v = 3;
    kanal_num_d = 9;
    modulo = 1/4; % [s]
    max_index = size(Messreihe1, 1);
    min_index = 1;
    %-----------------------------------------------------------------
    index = fix((t-39.25)/modulo)+1;
    index = min_index + (index>min_index).*(index-min_index);
    index = max_index + (index<max_index).*(index-max_index);

    v = 2*Messreihe1(index,kanal_num_v)';
    d = Messreihe1(index,kanal_num_d)';

else

    %--------------------- constants ----------------------
    kanal_num_v = 3;
    kanal_num_d = 9;
    modulo = 1/4; % [s]
    t_min = 410;
    t_max = 6310;
    t_cor = t; t_cor(t<t_min) = t_min; t_cor(t>t_max) = t_max;
    max_index = size(Messreihe1_zyklisch, 1);
    %-----------------------------------------------------------------
    index = fix((t_cor-t_min)/modulo)+1;
    index = mod(index-1,max_index)+1;
    
    v = 2*Messreihe1_zyklisch(index,kanal_num_v)';
    v((t<t_min) | (t>t_max)) = 0;

    d = Messreihe1_zyklisch(index,kanal_num_d)';
    d((t>t_max)) = 0;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute ambient temperature 
% (shifted to guarantee a zero initial condition for reduced order models)

function T = sensor_ambientTemperature(t, T_Grund, Messreihe2)

%--------------------- constants ----------------------
kanal_num = 22;
modulo = 10; % [s]
max_index = size(Messreihe2, 1);
%-----------------------------------------------------------------
index = fix(t/modulo)+1;
index = max_index + (index<max_index).*(index-max_index);

T = Messreihe2(index,kanal_num)'-T_Grund;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute motor power loss
function g = sensor_motorPowerLoss(t, c_pool, VARIANTE1, Messreihe1,...
                                   Messreihe1_zyklisch)

Factor1 = 6;
    
[v,~] = achse1_velocity_distance(t,VARIANTE1,Messreihe1,Messreihe1_zyklisch);
w = 2*pi/c_pool.sp_p*abs(v);

M_Mutter = c_pool.mu_M;
M_0      = (4501*c_pool.l_f_0*c_pool.l_d_m^3)*(c_pool.l_ny_40*w).^(2/3);
M_Lager  = M_0 + c_pool.l_M_1 + c_pool.l_M_d;

g = Factor1*0.24*(M_Mutter+M_Lager).*w;

end
