% Imie Nazwisko Indeks
% Grupa nr. 1, Wtorek, 7:30
% LAB 14
% MATLAB version: 23.2.0.2436196 (R2023b) Update 6

% Wersja 3 25I2024
close all; clear;

%%  ======================= Zadeklarowanie zmiennych z miniprojektu =======================
indeks = [1,2];                % indeks
zm_wym_powi = 0.6;             % Wymiana Powietrza [%]
zm_str_dach = 0.4;             % Straty przez Dach [%]
% Zmienna skalująca wartości Nominalne
zm_Qg   = 1.0;                 % zmienna mocy  % 1.0, 0.7, 0.3
zm_Tzew = 0;                   % zmienna temp  % 0,   10,  20
zm_fp   = 1.0;                 % zmienna fp    % 1.0, 0.7, 0.3
% Zmienne skoku
zm_skok_dTzew = 5;             % 0 lub np. 1
zm_skok_dQg   = 0;             % 0 lub np. 0.05*QgN (5% zakresu)
zm_skok_dfp   = 0;             % 0 lub np. 0.2*fpN (20% zakresu)
%Wartości Nominalne
TzewN = -20;                    % °C
TwewN =  20;                    % °C
TpN =    10;                    % °C
QgN = 20000;                    % W
Vw = indeks(1)*indeks(2)*2.5;   % m3    objętość wewnętrzna
Vp = Vw*0.5;                    % m3    objętość podasza
cp = 1000;                      % kg/m3
ro = 1.2;                       % kg/m3
%Identyfikacja parametrów statycznych
fpN = Vw / (12*60*60);  % m^3
Kp = (  (cp*ro*fpN)*(TpN-TzewN)     )/(  zm_str_dach*(TwewN-TpN)  );
K2 = (  zm_wym_powi*Kp*(TwewN-TpN)  )/(  TpN-TzewN    );
K1 = (  QgN-Kp*(TwewN-TpN)          )/(  TwewN-TzewN  );
% parametry "dynamiczne"
Cvw = cp * ro * Vw;
Cvp = cp * ro * Vp;
fp0   = fpN   * zm_fp  ;% - zm_skok_dfp;

%% =======================      V Transmitancje      =======================
  M=[Cvw*Cvp, (Kp+K2+cp*ro*fp0)*Cvw+(K1+Kp)*Cvp, K1*(Kp+K2+cp*ro*fp0)+Kp*(K2+cp*ro*fp0)];
  Qg_1   = [Cvp, K2+Kp+cp*ro*fp0];
% Tzew_1 = [Cvp*K1, K1*K2+K1*Kp+K2*Kp+K1*cp*ro*fp0+Kp*cp*ro*fp0];
  Qg_2   = Kp;
% Tzew_2 = [Cvw*K2+Cvw*cp*ro*fp0, K1*K2+K1*Kp+K2*Kp+K1*cp*ro*fp0+Kp*cp*ro*fp0];

%% ======================= LAB 14 =======================
s = tf('s');
A = M(1);
B = M(2);
C = M(3);

S1 = ( -B-sqrt(B^2-4*A*C) )/( 2*A );
S2 = ( -B+sqrt(B^2-4*A*C) )/( 2*A );

%T1 = ( -1 )/( S1 );
%T2 = ( -1 )/( S2 );
T1 = ( -1 )/( S1 ) - 73.0552563151182;
T2 = ( -1 )/( S2 ) + 97815.4302563151;

D = A* S1* S2;
E = Qg_1(1);
F = Qg_1(2);
k = Qg_2;
%T3 = E/F;
T3 = E/F - 1160;

G1 = ( (s*T3+1)*F )/( D*(s*T1+1)*(s*T2+1) );
G2 = ( k          )/( D*(s*T1+1)*(s*T2+1) );
% G1 = tf( (s*T3+1)*F/D, (s*T1+1)*(s*T2+1) );
% G2 = tf( k/D,          (s*T1+1)*(s*T2+1) );
% G1 = tf( [Qg_1(1),Qg_1(2)], M);
% G2 = tf(  Qg_2,           M);

%% ======================= Rysowanie wykresów =======================
% 1. Rysujemy poziomą linię, jest na wysokości 20log(k)
% 2. Rysujemy linię:
%      Tak by miała 20 wysokości (np. od -60 do - 80)
%      Tak by miała od 10^-4 do 10^-3.
%      Co daje 20 db/dek
% 3. Przystawiamy ją do wykresy.

figure(1);
%step(G1); grid on;
bode(G1); grid on;
title("Bode Diagram ( (s*T3+1)*F )/( D*(s*T1+1)*(s*T2+1) )");
figure(2);
nyquist(G1); grid on;

figure(3);
%step(G2); grid on;
bode(G2); grid on;
title("Bode Diagram ( k )/( D*(s*T1+1)*(s*T2+1) )");
figure(4);
nyquist(G2); grid on;

% Wartości z wykresu zgadzają się z wartościami z matlaba.
disp("k1 = "+F/D+"        20*log10(k1) = "+20*log10(F/D) );
disp("k2 = "+k/D+"       20*log10(k2) = "+20*log10(k/D) );
disp("T1 = "+T1+"      Frequency f1 = "+1/T1);
disp("T2 = "+T2+"    Frequency f2 = "+1/T2);
disp("T3 = "+T3+"         Frequency f3 = "+1/T3);
disp("G1 = ( (s*T3+1)*F )/( D*(s*T1+1)*(s*T2+1) )");
disp("G2 = ( k          )/( D*(s*T1+1)*(s*T2+1) )");