clc
clear all
close all

% analysis of RC circuit in parallel config

R1=100;
R2=1000;
C=1e-6;

% define a, b, c, d

a1=-1/(R1*C);
a2=-1/(R2*C);
tau1=-1/a1;
tau2=-1/a2;

b=1/C;
c=1;
d=0;

% define ss rep
sys1=ss(a1,b,c,d);
sys2=ss(a2,b,c,d);

% initial condition
x0=1;

% free evolution
figure (1)
initial(sys1,x0)
hold on
grid on
initial(sys2,x0)
legend('ss1','ss2')

figure(2)
initial(sys1,sys2,x0)
hold on
grid on
legend

% simulate the step resonse

figure(3)
step(sys1,sys2)
grid on
hold

% define transfer function of RC 

s=tf('s');
W1=b*c/(s-a1)+d;
W2=b*c/(s-a2)+d;

num_1=1/C;
den_1=[1 -a1];

num_2=1/C;
den_2=[1 -a2];

W1a = tf(num_1,den_1);
W2a=tf(num_2,den_2);

[num_1b,den_1b] = ss2tf(a1,b,c,d);  % =W1a; =W1
 

% plot the step response
figure(4)
step(W1,W2)
hold on
grid on

% plot the response to a step signal with amplitude U0
figure(5)
U0=2;
opt = RespConfig;
opt.Amplitude = U0;
step(sys1,sys2,opt)

% a different way
opt1 = stepDataOptions('StepAmplitude',2);
figure(6)
step(sys1,sys2,opt1);

%  compute the pole= a defined by ss
pole_1=pole(W1); % pole(sys1)
pole_2=pole(W2); % pole(sy2)

% compute the static gain equal to c*b/(-a)
W0_1=dcgain(W1); 
W0_2=dcgain(W2);