%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   slug_3D_linearize.m
%
%   Script for linearizing the simplified slug model made as 
%   an aid in tuning the model
%
%   Written by :    Espen Storkaas (espen.storkaas@chemeng.ntnu.no)
%   Date:           08.10.02
%   
%   For questions, comments or documentation, contact author
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
global data 

% Inputs for initialization
%=========================================================================\
%Definerer parametere
%=========================================================================/
mG_in=0.362; %l/min denne blir satt også i s_slug!!!!!!! bør kanskje prøve å sette den bare ett sted
mL_in=8.64; %l/min denne blir satt også i s_slug!!!!!!! bør kanskje prøve å sette den bare ett sted
z_setpkt=0.3;%  %ventilåpning det skal lineariseres rundt 
z_arbpkt=70;%Arbeidspunkt for pådraget i ytre sløyfe. Se simulinkdiagram stasjonaer.sim
%=========================================================================\
%Simulerer for å finne P1, P2 og h1 stasjonært
%=========================================================================/
simtid=3*3600;
sim('stasjonaer');
P1_stationary=1e5*P1_simulert.signals.values(length(P1_simulert.signals.values)); %henter ut siste verdi
P2_stationary=1e5*P2_simulert.signals.values(length(P2_simulert.signals.values)); %henter ut siste verdi
h1_stationary=h1_simulert.signals.values(length(h1_simulert.signals.values)); %henter ut siste verdi
%=========================================================================\
%Finner den stasjonære løsningen x0 rundt z=0.3
%=========================================================================/
[x0,y_stasj,data]=initialize2(z_setpkt,mG_in,mL_in,P1_stationary,h1_stationary,P2_stationary); 

%========================================================================
% sys=slug_3D_lin(0,x0,[z_setpkt;mG_in;mL_in]); %bare en test for å se om
% x0 faktik er en stasjonærløsning returverdien bør være ca. lik 0.
% disp(sys)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Løsning på hvorfor jeg får "trying to linearize around a point which is   
% not in steadystate" : x0 i steady-state er funnet ved hjelp av funksjonen
% initialize2. Denne funksjonen tar inn en rekke parametere, deriblant
% mL_in og mG_in (De andre verdiene er funnet vha simulering i simulink, med en kaskaderegulator). mL_in og mG_in har her verdier forskjellig fra 0. Når jeg
% skal linearisere ser jeg imidlertid på et system som kun har EN input,
% nemlig u=z. Problemet oppstår når jeg vil linearisere rundt et system med
% kun EN input. Jeg tenkte at jeg kunne sette u=z_setpkt i funksjonskallet
% til linearize. Det er imidlertid dette som forårsaker problemene. Jeg er
% imidlertid usikker på om dette faktisk har noe å si for resultatet av
% transferfunksjonene. Skal jeg finne en transferfunksjon, settes jo som
% regel alle andre innganger til 0. Kan prøve å finne transferfunksjoner
% fra hele u i stedet for bare fra u=z_setpkt. Dette gjøres ved å plassere
% alle de alternative målingene i målevektoren Y, og sette IU=1 (velge u=z)
% når jeg benytter meg av funksjonen ss2tf i stedet. 
%
% Konklusjon:
% Forskjellen er altså
% at jeg først ønsket å lage en tilstandsrom-modell med en input u=z for så
% å finne en transferfunksjon. Nå finner jeg en tilstandsrom-modell med 3
% input (u=[z mG_in mL_in]), og velger hvilken inngang jeg ser på ved bruk
% av argumentet IU når jeg benytter meg av ss2tf. (I dette tilfellet er IU=1 pga u=z)
% Til slutt kan det nevnes at jeg foretok en sammenligning av de to
% fremgangsmåtene. Resultatet ved bruk av de to metodene var identisk til tross for at jeg ved bruk av den
% første metoden fikk feilmeldingen "trying to linearize around a ..etc". 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=========================================================================\
%Lineariserer rundt den stasjonære løsningen x0 for z=0.3
%=========================================================================/
if ~isreal(y_stasj);
    disp('something is wrong')
end
u=[z_setpkt;mG_in;mL_in]; %u0 = u da u er konstant.
%Linearizing
[A,B,C,D]=linearize(@slug_3D_lin,[],x0(1:3),[],u);

disp('================================================================')
disp('Polene til systemet')
disp(eig(A))
%=========================================================================\
%Finner transferfunksjonen mellom u=z og y=[DP rho_T w Q alpha_L alpha_Lm]
%=========================================================================/
iu=1; % skal ha transferfunksjonen fra iu`th input - dvs. z.
[tf_teller,tf_nevner]=ss2tf(A,B,C,D,iu); %
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=DP')
h1=tf(tf_teller(1,:),tf_nevner) %DP
disp('Nullpunkt')
roots(tf_teller(1,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=rho_T')
h2=tf(tf_teller(2,:),tf_nevner) %rho_T  
disp('Nullpunkt')
roots(tf_teller(2,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=w')
h3=tf(tf_teller(3,:),tf_nevner) %w   
disp('Nullpunkt')
roots(tf_teller(3,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=Q')
h4=tf(tf_teller(4,:),tf_nevner) %Q
disp('Nullpunkt')
roots(tf_teller(4,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=alpha_L')
h5=tf(tf_teller(5,:),tf_nevner) %alpha_L
disp('Nullpunkt')
roots(tf_teller(5,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=alpha_Lm')
h6=tf(tf_teller(6,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(6,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=rho_T*DP')
h7=tf(tf_teller(7,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(7,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=z*DP/rho_T')
h8=tf(tf_teller(8,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(8,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=DP/rho_T')
h9=tf(tf_teller(9,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(9,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=z*rho_T*DP')
h10=tf(tf_teller(10,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(10,:))
disp('================================================================')
disp('Transferfunksjon fra u=z -> y=rho_T*DP')
h11=tf(tf_teller(11,:),tf_nevner) %alpha_Lm
disp('Nullpunkt')
roots(tf_teller(11,:))
disp('================================================================')
%=========================================================================\
% Måling y=DP/rho
%=========================================================================/
disp('Lineærkombinasjon av DP - rho')
num_h1=tfdata(h1,'v');
num_h1=num_h1*1/data.lin_rho;
num_h2=tfdata(h2,'v');
num_h2=(num_h2*data.lin_DP)/(data.lin_rho^2)
H_num_dp_minus_rho=num_h1-num_h2;
H_dp_minus_rho=tf(H_num_dp_minus_rho,tf_nevner)
disp('Nullpunkt');
roots([H_num_dp_minus_rho])
disp('================================================================')

%=========================================================================\
% Måling av y=DP*rho
%=========================================================================/
disp('Lineærkombinasjon av DP+rho')
num_h1=tfdata(h1,'v');
num_h1=num_h1*1/data.lin_rho;
num_h2=tfdata(h2,'v');
num_h2=(num_h2*data.lin_DP)/(data.lin_rho^2)
H_num_dp_pluss_rho=num_h1+num_h2;
H_dp_pluss_rho=tf(H_num_dp_pluss_rho,tf_nevner)
disp('Nullpunkt');
roots([H_num_dp_pluss_rho])
disp('================================================================')