% EExtra_3.m
% Solution to problem 3 of ConvExtra. JMS, Jun.2024

clear

% Set parameters
V0 = 1e-3;    % source voltace amplitude (V)
L1 = 1e-7;    % 1st self inductance (H)
L2 = 3e-8;    % 2nd self inductance (H)
M = 2e-8;     % mutual inductance (H)
w = 2*pi*1e6; % source voltage frequency (rad/s)
dt = 5e-9;    % time step (s)
tf = 1.5e-5;  % final time (s)

% Integrate differential equation
Q0 = [0;0];     % initial condenser charges (C)
dQdt0 = [0;0];  % initial condenser intensities (C/s)
m = 1;          % fictitious mass
t = 0:dt:tf;
myforce = @(Q,dQdt,t)EExtra_VarCorr( Q, V0*sin(w*t) );
[Q,dQdt,d2Qdt2] = newtonRK( myforce, m, Q0, dQdt0, dt, tf );

% Plot charges
figure(11)
plot(t,Q,'LineWidth',2)
legend('Q_1','Q_2')
xlabel('time (s)')
ylabel('charge (C)')
set(gca,'FontSize',15)
grid on

% Plot intensities
figure(12)
plot(t,dQdt,'LineWidth',2)
legend('I_1','I_2')
xlabel('time (s)')
ylabel('dQ/dt (A)')
set(gca,'FontSize',15)
grid on

% Find and plot powers
P1 = L1*dQdt(1,:).*d2Qdt2(1,:);
P2 = L2*dQdt(2,:).*d2Qdt2(2,:);
PM = 0.5*M*( dQdt(1,:).*d2Qdt2(2,:) + dQdt(2,:).*d2Qdt2(1,:) );
figure(13)
plot( t,P1, t,P2, t,PM, 'LineWidth',2 )
legend('P_1','P_2','P_M')
xlabel('time (s)')
ylabel('power (W)')
set(gca,'FontSize',15)
grid on
