clear all 
clc
global data
load Stationary
load data
OPz=30;
index=find(abs(z-OPz)==min(abs(z-OPz)));
P1(index)
mG_in=(6.1717e-5)*1.5;
P_feed=0.31e5+1.013e5;
P0=1.013e5;
u=[z(index)*0.01;mG_in;P_feed;P0/1e5];
%Linearizing
[A,B,C,D,E,F]=linearize(@slug_3D_lin,[],xs(1:3,index),[],u(1),u(2:4));
G=ss(A,B,C,D);
G_tf=tf(G);

GP1=G_tf(1);
GP1=zpk(GP1);
GP2=G_tf(2);
GP2=zpk(GP2);
GpT=G_tf(3);
GpT=zpk(GpT);
GW=G_tf(4);
GW=zpk(GW);
GQ=G_tf(5);
GQ=zpk(GQ);

[num_GP1,den_GP1]=tfdata(GP1,'v');
[num_GP2,den_GP2]=tfdata(GP2,'v');
[num_GpT,den_GpT]=tfdata(GpT,'v');
[num_GW,den_GW]=tfdata(GW,'v');
[num_GQ,den_GQ]=tfdata(GQ,'v');

% % %PI regulator m/prosess variabel P1
% Kc=-25;%bar^-1;
% taui=10;%s
% Gc=Kc*tf([taui 1],[taui 0]);
% S=inv(1+GP1*Gc);
% disp('Poler S')
% pole(S)
% T=1-S;
% figure(2)
% bodemag(S),hold on,bodemag(T),hold off,legend('S','T')
% W=linspace(1E-3,1E3,40000);
% [MAG,PHASE] = BODE(S,W);
% for i=1:length(MAG)
%     mag(i)=MAG(:,:,i);
% end
% loglog(W,mag)
% axis([1E-3 1E3 1E-5 10])
% hold on
% [MAGT,PHASE] = BODE(T,W);
% for i=1:length(MAGT)
%     magT(i)=MAGT(:,:,i);
% end
% loglog(W',magT,'g')
% axis([1E-3 1E3 1E-5 10])
% hold off
% xlabel('Frekvens (rad/s)')
% ylabel('Magnitude (abs) ')
% legend('S','T')
% 
% 
% g=gcf;
% set(g,'PaperUnits','centimeters','PaperType','A4')
% opts = struct('FontMode','fixed','FontSize',9,'Width',10,'Color','rgb');
% exportfig(g,'bode.eps',opts)


%PI regulator m/prosess variabel W
% Kc=-15;%bar^-1;
% taui=10;%s
% Gc=Kc*tf([taui 1],[taui 0]);
% S=inv(1+GW*Gc);
% disp('Poler S')
% pole(S)
% T=1-S;
% disp('M_S');
% norm(S,inf)
% disp('M_T');
% norm(T,inf)



% %transfer funksjon for:
% %  Y     output( P1,P2,p.T,W, Q)
% % ---    
% %  d     forstyrrelse (mG_in, mL_in)
% Gd=ss(A,E,C,F);
% Gd=Gd(:,1:2);
% Gd_tf=tf(Gd);
% %transfer funksjon for:
% %  P2    output
% % ---    
% %  z    ventilaapning

 poler_GP1=pole(GP1)
 nuller_GP1=zero(GP1)
% 
 poler_GP2=pole(GP2)
 nuller_GP2=zero(GP2)
% 
 poler_GpT=pole(GpT)
 nuller_GpT=zero(GpT)
% 
 poler_GW=pole(GW)
 nuller_GW=zero(GW)
% 
 poler_GQ=pole(GQ)
 nuller_GQ=zero(GQ)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Kaskade struktur
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% G1=GP2;
% G2=GQ;
% Kc1=-0.115;tau1=700;
% Kc2=-0.001;tau2=1200;
% 
% if tau1==0
%     K1=Kc1;
% else
% K1=Kc1*tf([tau1 1],[tau1 0]);
% end
% 
% if Kc2==0
%     K2=Kc2;
% else
% K2=Kc2*tf([tau2 1],[tau2 0]);
% end
% L1=K1*G1;
% L2=K2*G2;
% S1=1/(1+L1);
% T1=L1/(1+L1);
% S2=1/(1+L2);
% T2=L2/(1+L2);
% G=G1*S2*K2;
% 
% L=K1*G;
% S=1/(1+L);
% T=L/(1+L);

