function [C1,C2,C3,R2,R3] = thirdOrderPassiveFilterDesign(Fc,Phi,Icp,Kvco,Nmax)

% This function computes the 3rd order passive loop filter componet values 
% for the 3 capacitors (C1 , C2 & C3) and 2 resistors (R2 & R3).    
%
%Input:
% Fc   - Loop bandwidth(ideally start at 1/10th of the reference frequency)
% Phi  - Desired phase margin (ideally it should be 40< phi <60)
% Icp  - Charge-pump current (follow the design spec sheet value)
% Kvco - VCO voltage sensitivity value (in Hz/V) mentioned in design spec.
% Nmax - Maximum divider ratio the PLL is supposed to work on.
%
% Ouput:
% C1, C2, C3, R2 & R3 : Component values of three capacitors, C1, C2 & C3, and  
%                       two resistors, R2 & R3, for 3rd order passive loop filter.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The transfer function of the 3rd order filter is represented as:
%                         1+s.T2
%            Z(s)= ---------------------
%                  A2.s^3 + A1.s^2 + A0.s
% where, T2 = R2*C2 
%        A0 = C1+C2+C3
%        A1 = A0*(T1+T3)
%        A2 = A0*T1*T3
%
% There are parameters called 'gamma optimization factor' and 'Pole Ratio (T31)' 
% that are used to optimize the performance of the PLL and solve for loop filter
% equation above.We will be using the following standard default values for
% these two parameters:
%   gamma = 1
%   T31   = T3/T1 = 0.6 

% Copyright 2018, The MathWorks, Inc.

if(nargin ~= 5)
    error(message('msblks:msblksMessages:InvalidNumberOfParameters'));
elseif (length(Fc) ~= 1 || Fc <= 0 || isreal(Fc) ~= 1)
    error(message('msblks:msblksMessages:InvalidLoopBandwidth'));
elseif (length(Phi) ~= 1 || Phi <= 0 || isreal(Phi) ~= 1)
    error(message('msblks:msblksMessages:InvalidPhaseMargin'));
elseif(length(Icp) ~= 1 || Icp <= 0 || isreal(Icp) ~= 1)
    error(message('msblks:msblksMessages:InvalidChargePumpCurrent','Icp'));
elseif(length(Kvco) ~= 1 || Kvco <= 0 || isreal(Kvco) ~= 1)
    error(message('msblks:msblksMessages:InvalidVoltageSensitivity','Kvco'));
elseif(length(Nmax) ~= 1 || Nmax <= 0 || isreal(Nmax) ~= 1)
    error(message('msblks:msblksMessages:InvalidDividerRatio','Nmax'));
% elseif(nargin == 5)
%     warning(message('msblks:msblksMessages:MissingPlotOption','Plot'));
%     Plot=0;
% elseif((Plot ~= 0 && Plot ~= 1))
%     warning(message('msblks:msblksMessages:InvalidPlotOption','Plot'));
%     Plot=0;
end

Kd=Icp;
N=Nmax;
% converting Fc ,loop bandwidth, to rad/s 
OmegaC=2*pi*Fc;
% setting the optimization paramters gamma and T31 to default values
T31=0.6;
Gamma=1;

% computing pole T1 by numerical method with the use of symbolic math toolbox
% syms x;
% T1=double(vpasolve(atan(gamma/(wc*x*(1+T31)))-atan(wc*x)-atan(wc*x*T31)==(pi*phi)/180));

% Another approach to solve for T1 usin optimization toolbox
%T1=fsolve(@(x)(atan(gamma/(wc*x*(1+T31)))-atan(wc*x)-atan(wc*x*T31)-(pi*phi)/180),1);

%T1=double(T1); % converting from symbolic value to double

%computing poles T1 & T3 using approximate method
T1=(sec((pi/180)*Phi)-tan((pi/180)*Phi))/(OmegaC*(1+T31));
T3=T1*T31;
%computing the zero (T2) of loop filter which is independent of the loop filter order.
T2=Gamma/(OmegaC^2*(T1+T3));
% Using the pole and zero values already calculated above we will compute the
% filter coefficents A0, A1 and A2
A0=(Kd*Kvco/(OmegaC^2*N))*sqrt((1+(OmegaC^2*T2^2))/((1+(OmegaC^2*T1^2))*(1+(OmegaC^2*T3^2))));
A1=A0*(T1+T3);
A2=A0*T1*T3;
% Now will use the filter coefficents computed above to find the respective
% Rs and Cs for the 3rd order filter
C1=(A2/(T2^2))*(1+sqrt(1+((T2/A2)*(T2*A0-A1))));
C3=((T2*A1*C1)-(A2*A0)-(T2^2*C1^2))/(T2^2*C1-A2);
C2=A0-C1-C3;
R2=T2/C2;
R3=A2/(C1*C3*T2);

if(C1<0 || isreal(C1) ~= 1 || isnan(C1) == 1 || C2<0 || isreal(C2) ~= 1 || ...
   isnan(C2) == 1 || C3<0 || isreal(C3) ~= 1 || isnan(C3) == 1 || R2<0 || ...
   isreal(R2) ~= 1 || isnan(R2) == 1 || R3<0 || isreal(R3) ~= 1 || isnan(R3) == 1)
    error(message('msblks:msblksMessages:InvalidLoopFilterComponents'));
end
%coder.extrinsic('tf','tf2zp','margin');
% using the transfer function we will find the pole , zeros and gain values
% Num=[T2 1];
% Den=[A2 A1 A0 0];
% LoopFilterTf=tf(Num,Den); % loop-filter transfer function
% NumVco=Kvco;
% DenVco=[1 0];
% VcoTf=tf(NumVco,DenVco);
% ForwardGain=(Kd*VcoTf*LoopFilterTf); % forward gain = Phase_detector * Loop-filter tf * VCO gain
% BackwardGain=1/N; % Backward gain comprises of the divider value
% OpenLoopTf=ForwardGain*BackwardGain;% open-loop transfer function
% if(Plot==1)
%     PlotTitle=message('msblks:msblksMessages:LoopDynamicPlotTitle').getString;
%     figure('Name',PlotTitle,'IntegerHandle','off','Tag','ThirdOrderLoopDynamicPlot');
%     subplot(2,2,1);
%     margin(OpenLoopTf); 
%     closed_sys=feedback(ForwardGain,BackwardGain);
%     subplot(2,2,2);
%     step(closed_sys);
%     subplot(2,2,4);
%     impulse(closed_sys);
%     subplot(2,2,3);
%     iopzmap(closed_sys);
% end
