% Control3A_2.m
% Solution to problem 2 of control 3. May 14, 2021.
% JMS, May.2021

% Set data
clear
R = [0.5e-3;1e-3];   % circuit resistances (Ohm)
C = [120e-6;80e-6];  % circuit capacities (F)
L = [50e-9;100e-9];  % circuit self-inductances (H)
M = 30e-9;           % mutual inductance (H)
q0 = 1e-6;           % initial charge of C1 (C)
tf = 100e-6;         % simulation time (s)
dt = 0.02e-6;        % simulation time interval (s)

% Solve circuit's differential equations
Current_Variation = @(Q,I,t) ...
  [L(2)*(Q(1)/C(1)+R(1)*I(1)) - M*(Q(2)/C(2)+R(2)*I(2)); ...
   L(1)*(Q(2)/C(2)+R(2)*I(2)) - M*(Q(1)/C(1)+R(1)*I(1))] / (M^2-L(1)*L(2));
Q0 = [q0;0];
I0 = [0;0];
[Q,I] = newtonRK( Current_Variation, 1, Q0, I0, dt, tf );

% Plot charges as a function of time
figure(11)
t = 0:dt:tf;
plot(t,Q,'LineWidth',1.5)
xlim([0,tf])
legend('Q_1','Q_2')
xlabel('t (s)')
ylabel('Q (C)')
set(gca,'FontSize',15)
grid on

% Plot intensities
figure(12)
plot(t,I,'LineWidth',1.5)
xlim([0,tf])
legend('I_1','I_2')
xlabel('t (s)')
ylabel('I (A)')
set(gca,'FontSize',15)
grid on

% Plot energies dissipated in each resistance
figure(13)
P = R.*I.^2;
W(1,:) = integral(t,P(1,:));
W(2,:) = integral(t,P(2,:));
plot(t,W,'LineWidth',1.5)
xlim([0,tf])
legend('W_1','W_2','Location','SouthEast')
xlabel('t (s)')
ylabel('W (J)')
set(gca,'FontSize',15)
grid on
