%% Application: State and Output Trajectories of LTI Systems using Laplace Transform - Example 1
%
% Table of Contents
% Introduction
% LTI System of 1st Order: a Series RC Circuit
%     Step 0: Declare the Variables
%     Step 1: Declare the LTI Equations
%     Step 2: Solve the ODE Using the Laplace Transform
%     Step 3: Substitute Numerical Values
%     Step 4: Plot the results.
%
%
%% Introduction
% Symbolic workflows keep calculations in their natural symbolic form rather than numerical form. This approach helps to understand the solution's properties and uses exact symbolic values. Numbers are substituted for symbolic variables only when a numerical result is required or cannot be continued in a symbolic manner.
%
% Typically, the steps are:
% Declare the variables and the (state and output) equations.
% Solve the equations using the Laplace transform.
% Substitute numerical values.
% Plot the results.
% Analyse the results.
%
%% LTI System of 1st Order: a Series RC Circuit
%
% We want to determine state and output trajectories for the LTI system:
%
%
% \left\{ \begin{array}{rcl}
% \dot{x}(t) &=& \frac{1}{C} u(t)\\
% y(t) &=& x(t) + Ru(t)
% \end{array} \right.
% with initial condition x(0) = x_)\, and input u(t) = A \cos(\omega t) \cdot 1(t)\,
%
%
%% Step 0: Declare the Variables

clear

syms t s x(t) y(t) u(t) X Y C R A omega % the variables 
                                      % in the time and frequency domain

assume(t >=0)
assume([A C R omega] >0) % a compact way to declare A, C, R and omega as positive real
                         % symbolic variables


u(t) = A * cos(omega * t) % Note: you DO NOT NEED to explicity write the Heaviside function!
x0VAL = -1.45; % the initial condition [V]
Aval = 2.50;  % V
Cval = 1.5;   % F
Rval = 10.0;  % Ohm
OMval = 2*pi*0.1; % rad/s

%%
%
%% Step 1: Declare the LTI Equations
%
% Firstly, declare the symbolic derivative \dot x of the state x(t)

dot_x = diff(x,t);

% Then, declare the ODE state equation and the algebraic output equation. Pay attention to the == symbols

eqnX = dot_x  == (1/C) * u; % state eq.
eqnY = y == x+R*u;              % output eq.

% What does the equations look like?

disp(eqnX)
disp(eqnY)

%
%% Step 2: Solve the ODE Using the Laplace Transform
%
% Let's compute the Laplace transform of the ODE:

eqnLT = laplace(eqnX, t, s);

% What does the equation in the Laplace transform domain look like?

disp(eqnLT)

%
% Note: In order to continue and solve the algebraic equation with the Laplace transform, we must replace all terms indicating Laplace transforms of the function x(t) with the symbolic variable X corresponding to the Laplace transform X(s) .
%

eqnLT1 = subs(eqnLT,laplace(x,t,s),X)

%
% What about the output equation? Let's compute the Laplace transform of the output equation:

eqnYLT = laplace(eqnY, t, s);
disp(eqnYLT)

% Note: In order to continue and solve the algebraic equation with the Laplace transform, we must replace all terms indicating Laplace transforms of the function x(t) and y(t) with the symbolic variables X, Y, corresponding to the Laplace transforms X(s) and Y(s).

eqnY_LT1 = subs(eqnYLT,[laplace(x,t,s) laplace(y,t,s)],[X Y]);
disp(eqnY_LT1)

%
% Finally, solve the ODE

Xsol = solve(eqnLT1, X)
Xsol = simplify(expand(Xsol))

%
% Let's apply the inverse Laplace transform:

xsol = ilaplace(Xsol)

% That's it!

eqnY_LT1 = subs(eqnY_LT1,X,Xsol)

Ysol = solve(eqnY_LT1,Y)
ysol = ilaplace(Ysol)

%
%% Step 3: Substitute Numerical Values
%

xSol_IC = subs(xsol, x(0), x0VAL)
xSol_IC = subs(xSol_IC, [A C R omega], [Aval Cval Rval OMval])

ySol_IC = subs(ysol, x(0), x0VAL)
ySol_IC = subs(ySol_IC, [A C R omega], [Aval Cval Rval OMval])

%
%
%% Step 4: Plot the results.
%

figure('Units','normalized','Position',[0.1,0.1,0.9,0.9])

subplot(2,1,1)
fplot(xSol_IC,[0 15], 'LineWidth',2)                      
title('State trajectory')
ylabel('x(t)')
xlabel('t')
grid on
subplot(2,1,2);
fplot(ySol_IC,[0 15], 'LineWidth',2)                        
title('Output trajectory')
ylabel('y(t)')
xlabel('t')
grid on

%%
%
%
%
%
%
%
%