classdef Solver < handle
    
 properties
Werkstoffmodell;
f_start;        % frequency of the excitation/ vibration [Hz] at the beginning
f_rate;         %[Hz/s] Steigerung der Frequenz
f_end;          % frequency of the excitation/ vibration [Hz] at the end
f;
sig_v;          % [MPa] Vorspannung
eps_n_r;        % Dehnung rechts
eps_n_l;        % Dehnung links
Fvr;
Fvl;
F;
dt;             % [s] Zeitschritt, time increase per step [s]
t_end;
t;
N;
F_max;          % [N] Kraftamplitude
Area;           % [m^2] Probenquerschnitt
m;              % Masse [kg]
c;              % [N*s/m] Dämfpung
k;              % [N/m] Steifigkeit
l_0;            % [m] Länge der Proben im entspannten Zustand
l_v;            % Laenge der  Feder vorgespannt
l_l;            % Laenge der  linken Feder zum Zeitpunkt t
l_r;            % Laenge der  rechten Feder zum Zeitpunkt t
Fl;             % Kraft aus Probe links
Fr;             % Kraft aus Probe rechts

%ODE45
t1;
y1;
x;

end

methods
        
function obj = Solver(argDuration,argCe,argK,argl_0,argArea,argM,argTimeIncrease,argForceAmplitude,argPreStress,argE_A,argE_M,argE_T,argE_n_K,argKprob,arg_eps_Tr_gr_r,arg_eps_Tr_gr_l,arg_sig_fh,arg_sig_fr,argH,argC,arg_delta,argT_0,argT_R,argA_s,argA_f,argM_s,argM_f,argFrequencyStart,argFrequencyIncrease,argFrequencyEnd)
obj.Werkstoffmodell = Werkstoffmodell(argE_A,argE_M,argE_T,argE_n_K,argKprob,arg_eps_Tr_gr_r,arg_eps_Tr_gr_l,arg_sig_fh,arg_sig_fr,argH,argC,arg_delta,argT_0,argT_R,argA_s,argA_f,argM_s,argM_f);      % Werkstoffmodell

% SOLVE  d2x/dt2*m + c*dx/dt + k*x + Fr(x) - Fl(x) = F(t)
% initial conditions: x(0) = 0, x'(0)=0 

obj.l_0 = argl_0; % [m] Länge der Proben im entspannten Zustand
obj.m = argM; %Masse [kg] 
obj.k = argK; %Federsteifigkeit [N/m] 
obj.c = argCe; %Dämpfung [N*s/m] 

obj.F_max = argForceAmplitude; %Anregungsamplitude [N]
obj.sig_v = argPreStress; %Vorspannung in [MPa]
obj.Area = argArea; %Balkenfläche [m²]

obj.dt = argTimeIncrease; %Zeitschritt [s] 
obj.f_start = argFrequencyStart; %Startfrequenz [Hz] 
obj.f_end = argFrequencyEnd; %Endfrequenz [Hz] 
obj.f_rate = argFrequencyIncrease; 
obj.f = obj.f_start:obj.f_rate*obj.dt:obj.f_end; 

obj.t_end = argDuration;
t=0:obj.dt:obj.t_end;   % time scale 


obj.F =  @(t) obj.F_max*chirp(t,obj.f_start,obj.t_end,obj.f_end,'lin',-90); %Anregungskraft in [N] 

figure
plot (t,obj.F(t))
title ('Force [N]')

obj.N = zeros (1,length(obj.F)+1); 

[obj.t1,obj.y1] = ode45(@(t,x) [x(2); (-x(2)*obj.c-obj.k*x(1)+obj.F(t))/obj.m] ,t ,obj.N); %Runge-Kutta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% obj.Fr =  @(x) Frfun(x); % Wie lässt sie diese Funktion mit dem Werkstoffmodell verknüpfen?
% obj.Fl =  @(x) Flfun(x); % Wie lässt sie diese Funktion mit dem Werkstoffmodell verknüpfen?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [obj.t1,obj.y1] = ode45(@(t,x) [x(2); (-x(2)*obj.c-obj.k*x(1)-obj.Fr(x)+obj.Fl(x)+obj.F(t))/obj.m] ,t ,obj.N); %Neue DGL

    end
    
end
end