% script for analyzing the RC circuit by varying the R value
clear all
close all
clc

% parameter values for the resistance (two cases: R1 and R2) and the capacitance
% (Cap)
R1=100;
R2=1000;
Cap=1e-6;

% define a (a1 and a2), b, c, d
a1=-1/(R1*Cap);
a2=-1/(R2*Cap);
b=1/Cap;
c=1;
d=0;

x0=10; % initial condition

% compute the time constant
tau1=R1*Cap;
tau2=R2*Cap;

% define the ss representation for the two circuits (R1-C and R2-C)
% xdot=a*x+b*u;
% y=c*x+d*u;
sysR1=ss(a1,b,c,d);
sysR2=ss(a2,b,c,d);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compute the free evolution %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
% simulate the free evolution
initial(sysR1,sysR2, x0);
hold on
grid on
legend('R1','R2')

% define time-vector
time=0:tau1/20:50*tau1;
%time_v=linespace(0,5e-3,500);

% store the free evolution ...
y_free1=initial(sysR1,x0,time);
y_free2=initial(sysR2,x0,time);

% ... and then plot
figure(2)
hold on
grid on
plot(time,y_free1,'-b','linewidth',2)
plot(time,y_free2,'-r','linewidth',2)

% another way to use the plot command  
% plot(time,y_free1,'-b',time,y_free2,'-r','linewidth',2)

xlabel('time')
title('free evolution')
ylabel('y_{free1},y_{free2}')
legend('y_{free1}','y_{free2}')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compute the step response %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
% drawing the step response
step(sysR1,sysR2)
hold on
grid on
legend('R1','R2')


figure(4)
hold on
grid on
% step response with amplitude 5e-3
% change the step options (by default the amplitude is set to 1)
opt = stepDataOptions('StepAmplitude',5e-3);
step(sysR1,sysR2,opt)
hold on
grid on
legend('R1','R2')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% compute the response to a square wave %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the response to a square wave with period equal to 5ms
tau_sq=5e-3;
[u_sq,time_sq] = gensig('square',tau_sq);

% plot the response to the generated square wave using the comand lsim
figure(5) 
lsim(sysR1,u_sq,time_sq)
hold on
grid on


% store the response to the square input...
y_sq=lsim(sysR1,u_sq,time_sq);

% ... and then plot the normalized output (y/R1 - with R1 the static gain of the circuit) 
figure(6)
hold on 
grid on
plot(time_sq,u_sq,'-k','linewidth',2);
plot(time_sq,y_sq/R1,'-m','linewidth',2);
xlabel('time')
ylabel('u,y')
title('response to a square wave')
legend('u','y')

% another way to use the plot command
figure(7)
hold on
grid on
plot(time_sq,u_sq,'-k',time_sq,y_sq/R1,'-m','linewidth',2);
xlabel('time')
ylabel('u,y')
title('response to a square wave')
legend('input u','output y')