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

%This function computes the 4th order passive loop filter componet values 
% for the 4 capacitors (C1, C2, C3 & C4) and 3 resistors (R2, R3 & R4).    
%
% 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, C4, R2, R3 & R4 : Component values of four capacitors, 
%                               C1, C2, C3 & C4, and three resistors, 
%                               R2, R3 & R4 for 4th order passive loop filter.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The transfer function of the 4th order filter is represented as:
%                             1+s.T2
%            Z(s)= -------------------------------
%                  A3.s^4 + A2.s^3 + A1.s^2 + A0.s
% where, T2= C2*R2
%        A0, A1 , A2 & A3 = filter coefficients
%
% There are parameters called 'gamma optimization factor' and 'Pole Ratios' 
% (T31 & T43)' that are used to optimize the performance of the PLL and 
% solve for loop filter equation above. Following standard default values for
% the parameters will be used:
%   Gamma = 1
%   T31   = T3/T1 = 0.5 
%   T43   = T4/T3 = 0.5

% 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, in rad/s 
OmegaC=2*pi*Fc;
% setting the optimization paramters gamma , T31 and T43 at default values
T31=0.4;
T43=0.4;
Gamma=1.115;

 % computing pole T1 using exact numerical method
%syms x;
%T1=double(vpasolve(atan(gamma/(wc*x*(1+T31+(T31*T43))))-atan(wc*x)-atan(wc*x*T31)-atan(wc*x*T31*T43)==(pi*phi)/180));

% computing for T1 using optimization toolbox

%T1=fsolve(@(x)(atan(gamma/(wc*x*(1+T31+(T31*T43))))-atan(wc*x)-atan(wc*x*T31)-atan(wc*x*T31*T43)-(pi*phi)/180),1);

Function=@(x)(atan(Gamma/(OmegaC*x*(1+T31+(T31*T43))))-atan(OmegaC*x)-atan(OmegaC*x*T31)-atan(OmegaC*x*T31*T43)-(pi*Phi)/180);
DifferentialFunction=@(x)(-OmegaC/(OmegaC^2*x^2 + 1) - (T31*OmegaC)/(T31^2*OmegaC^2*x^2 + 1) - ...
                          (T31*T43*OmegaC)/(T31^2*T43^2*OmegaC^2*x^2 + 1) - ...
                          Gamma/(OmegaC*x^2*(Gamma^2/(OmegaC^2*x^2*(T31 + T31*T43 + 1)^2) + 1)*(T31 + T31*T43 + 1)));
                      
T1=algebraicEqnSolver(Function,DifferentialFunction,1e-6,1e-6,1e3);

%computing poleS T1 & T3 using approximate method
%T1=(sec((pi/180)*phi)-tan((pi/180)*phi))/(wc*(1+T31));
T3=T1*T31;
T4=T1*T31*T43;
%computing the zero (T2) of loop filter which is independent of the loop filter order.
T2=Gamma/(OmegaC^2*(T1+T3+T4));
% 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))*(1+(OmegaC^2*T4^2))));

% Solving for components C1 and R3 using Pole T3
A1t3=A0*(T1+T3);
A2t3=A0*T1*T3;
% Now will use the filter coefficents computed above to find the respective
% Rs and Cs for the 4th order filter
C1t3=(A2t3/(T2^2))*(1+sqrt(1+((T2/A2t3)*(T2*A0-A1t3))));
C3t3=((T2*A1t3*C1t3)-(A2t3*A0)-(T2^2*C1t3^2))/(T2^2*C1t3-A2t3);
R3t3=A2t3/(C1t3*C3t3*T2);

% Similarly solving for C1 and R3 using Pole T4
A1t4=A0*(T1+T4);
A2t4=A0*T1*T4;

C1t4=(A2t4/(T2^2))*(1+sqrt(1+((T2/A2t4)*(T2*A0-A1t4))));
C3t4=((T2*A1t4*C1t4)-(A2t4*A0)-(T2^2*C1t4^2))/(T2^2*C1t4-A2t4);
R3t4=A2t4/(C1t4*C3t4*T2);

% Now we will find the actual values for C1 and R3

C1=(C1t3+C1t4)/2;
R3=(R3t3+R3t4)/2;

% Solve for C2

A1 = A0*(T1+T3+T4);
A2=A0*((T1*T3) + (T1*T4) + (T3*T4));
A3=A0*T1*T3*T4;
K0=(A2/A3)-(1/T2)-(1/(C1*R3))-((A0-C1)*T2*R3*C1)/A3;
K1=A1-(T2*A0)-(A3/(T2*R3*C1))-((A0-C1)*R3*C1);
a=A3/(T2*C1)^2;
b=T2+(R3*(C1-A0))+((A3/(T2*C1))*((1/T2)-K0));
c=K1-(K0*A3/T2);

C2=(-b+sqrt(b^2-4*a*c))/(2*a);

% Solve for other components 

C3= (T2*A3*C1)/(R3*(K0*T2*A3*C1-C2*(A3-R3*(T2*C1)^2)));
C4=A0-C1-C2-C3;
R2=T2/C2;
R4=A3/(T2*R3*C1*C3*C4);
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 || C4<0 || ...
   isreal(C4) ~= 1 || isnan(C4) == 1 || R2<0 || isreal(R2) ~= 1 || ...
   isnan(R2) == 1 || R3<0 || isreal(R3) ~= 1 || isnan(R3) == 1 || R4<0 || ...
   isreal(R4) ~= 1|| isnan(R4) == 1)

    error(message('msblks:msblksMessages:InvalidLoopFilterComponents'));

% using the transfer function we will find the pole , zeros and gain values
%Num=[T2 1];
%Den=[A3 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','FourthOrderLoopDynamicPlot');
%     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
end
end

