% Circuit constants
C1 = 50 * 10^(-9);
C2 = 1000 * 10^(-9);
R = 2 * 10 ^ 3;
w = 2 * pi * 50;
VMAX = 220 * sqrt(2);
% Simulation time limit
TMAX = 0.1;
% Simulation time resolution
TRES = 0.0001;
for ti = 0 : 1 : (TMAX / TRES) - 1
% tt: integer time index to use in vectors
% t : actual time value
t = ti * TRES;
% Intermediate variables that simplifies the expression
k = R * (C1 + C2) * w;
k21 = (k * k) + 1;
phi = atan(-1 / k);
% Calculate voltage sample
v_trans = VMAX * (C1 / (C1 + C2)) * (1 / k21) * exp(-w * t / k);
v_stead = VMAX * (C1 / (C1 + C2)) * (k * sqrt(k21) / k21) * cos((w * t) + phi);
v(ti + 1) = v_trans + v_stead;
% Calculate current sample
i(ti + 1) = v(ti + 1) / R;
end
% Index for plotting
tt = 0 : TRES : (TMAX / TRES) - 1;
% Truncate workspace variables
tt = tt(1 : TMAX / TRES);
i = i(1 : TMAX / TRES);
v = v(1 : TMAX / TRES);
% Plot current
hold off;
plot(tt, v, 'r');
hold on;
plot(tt, 1000 * i, 'b');
xlabel('Time (s)');
ylabel('Current and Voltage (mA and V)');