%% Assignment
%
% G(s) = \frac{10}{(1+10\,s)\,(1+2\,s) }\quad \Longrightarrow \quad
% \left\{ \begin{array}{rcl}
% g&=&0 \\
% \mu &=& 10 \; \Longrightarrow \; \mu_{dB} = 20\, \text{dB} \\
% \tau_1 &=& 10 \\
% \tau_2 &=& 2 \\
% \omega_1 &=& \frac{1}{10} \\
% \omega_2 &=& \frac{1}{2}
% \end{array} \right.
%
% G\left(j\,\frac{1}{10} \right) = ? \qquad G\left(j\,\frac{1}{2} \right) = ?
% \angle G(j\, \bar \omega) = -90^{\circ} \;\longrightarrow \; \bar \omega =?
%
% We would like to answer the questions using the Symbolic Toolbox or the Control Toolbox functions.
%

clear variables

%
%% Solving by Using the Symbolic Math Toolbox

syms s 
syms omega real
assume(omega >=0)

G(s) = 10/((1+10*s)*(1+2*s))

G_om(omega ) = subs(G, s, 1i*omega)

%
%% First Question: Magnitude and Phase of the Frequency Response at a Given Angular Frequency
% G\left(j\,\frac{1}{10} \right) = ? \qquad G\left(j\,\frac{1}{2} \right) = ?

G10 = G_om(1/10)

G2 = G_om(1/2)

%
% The magnitude of the frequency response at \omega = \frac{1}{10} rad/sec :

magG10sym = abs(G10)

% and the computation of the phase \angle G\left(j\, \frac{1}{10} \right)

angle(G10)

%
% The magnitude of the frequency response at \omega = \frac{1}{2} rad/sec :

magG2sym = abs(G2)

% and the computation of the phase \angle G\left(j\, \frac{1}{2} \right)

angle(G2)

%
%% Second Question: Solving a Trigonometric Equation
% \angle G(j\, \bar \omega) = -90^{\circ} \;\longrightarrow \; \bar \omega =?

phaseG(omega) = simplify(angle(G_om))

%%
%% Foolish Approach
% Let's rely on the MATLAB function solve:


EQ1 = phaseG ==-pi/2;

sol = solve(EQ1,omega,"ReturnConditions",true)
sol.omega
sol.conditions

% Note: the provided solution is nonsense. We must rewrite  the equation to be solved more clearly and efficiently.
%
%% Yet Another Foolish Approach
% Let rewrite the equation describing the problem, by explicitly  determining real and imaginary part of G(j \omega):

imagG_om = simplify(imag(G_om))
realG_om = simplify(real(G_om))

% Now, let's write the computation of \arg G(\j\,\omega) by explicitly computing the \arctan function, according to Part 8, pages 36-42 of the course material (in particular, refer to page 37)

phaseG2 = atan(simplify(imagG_om/realG_om))-pi

% Now, let's write the equation:

EQ2 = phaseG2 == -pi/2
solve(EQ2, omega)

% Note: now it should be clear why the previously presented approach was nonsense, and why there is no solution with the current approach:
% \tan \alpha is undefined when \alpha = (2k+1)\,\frac{\pi}{2} \quad , \; k \in \mathbb{Z}
% we need a different approach to solve the equation.
%
%% An Effective Approach
%
% Recall that (Part 8, page 40 of the course material)
%
%
% Thus, the straightforward way to solve the equation is to put the real part of the frequency response equal to zero
%
% \left\{ \begin{array}{rcl}
% G(j\,\omega ) &=& R(\omega) + j\, I(\omega) \\
%  &=& \left[ -\, \frac{10\,(20 \omega^2-1)}{(4 \, \omega^2+1)\,(100\,\omega^2+1)} \right] + j\, \left[ -\, \frac{120\,\omega}{(4 \, \omega^2+1)\,(100\,\omega^2+1)} \right] \\ \\
% \angle G(j\,\omega) &=& -\, \frac{\pi}{2}
% \end{array} \right.  \qquad \Longrightarrow \quad
% \left[ -\, \frac{10\,(20 \omega^2-1)}{(4 \, \omega^2+1)\,(100\,\omega^2+1)} \right]  = 0
%

OMsol = solve(realG_om==0)
barOM = sym(1)/sym(sqrt(20))
double(barOM)

% Refer to Part 8, page 59 of the course material.
%
%% Solving by Using the Control Toolbox

clear s

s = tf('s')
Gs = 10/((1+10*s)*(1+2*s))

%
%% First Question: Magnitude and Phase of the Frequency Response at a Given Angular Frequency

s10 = 1i*0.1; % <<-- NB
s2 = 1i*0.5; 

FreqResp_Gs10 = evalfr(Gs, s10)
FreqResp_Gs2 = evalfr(Gs, s2)

%


magG10 = abs(FreqResp_Gs10) 
double(magG10sym)

magG2 = abs(FreqResp_Gs2) 
double(magG2sym)

%

rad2deg(angle(FreqResp_Gs10))

rad2deg(angle(FreqResp_Gs2))

%
%% Second Question: Solving a Trigonometric Equation
%
% \angle G(j\, \bar \omega) = -90^{\circ} \;\longrightarrow \; \bar \omega =?
%
% There is no direct way to solve this equation. You can only find an approximate solution, by "exploring" the Bode diagram of the phase of the frequency response.
%
% Put the mouse cursor on the phase diagram, and move it, still staying on the curve, until you find a point close enough to the desired phase

figure;
bodeplot(Gs)

%