function [ C_chi ] = pid_chi(lambda)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

disp('CHIDAMBARAM TUNING CALCULATION:');
disp('  ')
disp('Controller structure:')
disp('  ')
disp('C=Kc(1 +( 1/(tau_i* s))+ tau_d* s)*((alpha*s +1)/(beta*s +1))')


%%
s=tf('s');
theta=1; %time delay

%%
G=model;
[num ,den]=tfdata(G);
den=den{1};
num=num{1};

b1=num(2);
b0=num(3);
a1=den(2);
a2=den(3);


p=-b1/b0;
Kp=b0/a2;
a1=a1/a2;
a2=1/a2;




%% Here the tuning parameter is calculated
% lambda=130;
c0=-1;
c1=-1;
c2=-1;
c3=-1;
c4=-1;
d1=-1;
d2=-1;



y1=(-a2^2)+(0.5*theta*p*a2); 
y2=(-0.5*theta*p*a1)-(a2*p)-(0.5*theta*a2);
y3=(-a1*a2)-(a2*p)-(0.5*theta*a2);
y4=(-0.5*theta*p)+a2;

z1=(a2*(lambda^3))+(1.5*theta*a2*(lambda^2))-(0.5*theta*a1*(lambda^3))-(3*lambda*(a2^2))-(theta*(a2^2))-(p*(a2^2));
z2=(3*a2*(lambda^2))+(1.5*theta*a2*lambda)-(0.5*p*theta*a2)-(0.5*theta*(lambda^3))-(3*a1*a2*lambda)-(a1*a2*theta)-(a1*a2*p);

%% Defining h, x, n and beta
n2=((y3*z1)-(y1*z2))/((y2*y3)-(y1*y4));
n1=((y2*z2)-(z1*y4))/((y2*y3)-(y1*y4));

h=(3*lambda)+theta-n1+p;

x1=((0.5*theta*(lambda^3))-(0.5*theta*n2*p))/h;
x2=((lambda^3)+(1.5*theta*(lambda^2))+(n2*p)+(0.5*theta*n2)-(0.5*theta*n1*p))/h;
x3=((3*(lambda^2))+(1.5*theta*lambda)-(n2)+(n1*p)+(0.5*theta*n1)-(0.5*p*theta))/h;

alpha=0.5*theta;
beta=x1/a2;


%% Assembling the controller and plotting responses

Kc=n1/(Kp*h);
taui=n1;
taud=n2/n1;




%% Stability stuff

%Cs
c4=(a2*beta*taui)+(0.5*theta*taui*taud*Kp*Kc*p);
c3=(a1*beta*taui)+(a2*taui)-(0.5*theta*taui*taud*Kp*Kc)+(0.5*theta*taui*Kp*Kc*p)-(taui*taud*Kp*Kc*p);
c2=(taui*beta)+(a1*taui)-(0.5*theta*taui*Kp*Kc)+(taui*taud*Kp*Kc)+(0.5*theta*Kp*Kc*p)-(taui*Kp*Kc*p);
c1=taui-(0.5*theta*Kp*Kc)+(taui*Kp*Kc)-(Kp*Kc*p);
c0=Kp*Kc;

%Ds
d1=((c1*c4)-(c2*c3))/c3;
d2=((c0*(c3^2))-(c4*(c1^2))+(c1*c2*c3))/((c1*c4)-(c2*c3));

C_chi=Kc*(1+(1/(taui*s))+(taud*s/(4*s+1)))*((alpha*s+1)/(beta*s+1)); %defining controller

end

