% LC_diode_cycle (see script to change parameters)
%
%   Simulates a DC-to-DC switcher circuit:
%
%     vin
%     |
%     L
%     |
%     ---- Diode >| --- C (to ground) --- R (to ground)
%     |
%     switch
%     |
%     ground
%
%   In particular, finds inductor current and
%   capacitor voltage trajectories. Circuit uses LC storage
%   characteristics to transform output impedance to input impedance.
%   That is, system comes to equilibrium where average input power
%   is equal to average output power but average input voltage does not
%   match average output voltage. In fact, average output voltage will
%   be higher than average input voltage.
%
%   When R is set to infinity, output voltage climbs forever. This
%   situation is a pathological singularity. Normally the device
%   "modulates" current and holds voltage constant. However, having
%   infinite output resistance forces output current to be zero, so the
%   device must modulate voltage. Input and output power *are* conserved
%   because output current is zero. In fact, there is an "equilibrium"
%   at infinite output voltage in order to maintain the power
%   equivalence.

% Copyright (c) 2007, 2008, 2009 by Theodore P. Pavlic
%
% This work is licensed under the Creative Commons
% Attribution-Noncommercial 3.0 United States License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc/3.0/us/ or
% send a letter to Creative Commons, 171 Second Street, Suite 300, San
% Francisco, California, 94105, USA.

function LC_diode_cycle

L = 10; % Inductance in Henries
C = 4; % Capacitance in Farads
R = 4; % Load resistance in Ohms ("Inf" works too)
vin = 5; % Input in Volts
f = 1; % 1 Frequency in Hz
w = 0.25; % Duty cycle (from 0 to 1)

tfinal = 250; % Simulation end time (in seconds)

control0 = 0; % Initial control (0 or 1)
i0 = 0; % initial inductor current (in A)
v0 = 0; % initial capacitor voltage (in V)

%% Simulate the system using ODE solver
%% Restrict step size so we catch every switch
wfix = abs( xor(1,control0) - w );
[tout,sysstep] = ode45( @(t,y) switcher(t,y,L,C,vin,f,wfix), ...
                        [0,tfinal], [i0;v0], ...
                        odeset( 'MaxStep', 1/(10*f) ));

% Plot capacitor and inductor waveforms from simulator
plot(tout, sysstep(:,1), '--', ...
     tout, sysstep(:,2), '-' );

% Plot control signal
t = linspace(0,tfinal,10000);
hold on
plot(t, ( xor( (t - floor(t*f)/f)*f >= wfix, control0 ) - 1.5 ), 'r-.' );
hold off

% Make the plot a little prettier
lfig = legend( '$i_L$ (A)', '$v_C$ (V)', 'Control', ...
               'Location', 'NorthWest' );
set(lfig, 'Interpreter', 'latex', 'FontSize', 14);
if( R == Inf )
    Rstr = '\infty$';
else
    Rstr = [ num2str(R) '$ $\Omega$' ];
end
xlimits = xlim;
ylimits = ylim;
text( 0.25*xlimits(2), 0.85*ylimits(2), [ '\shortstack{(' ...
         '$L = ' num2str(L) '$ H, ' ...
         '$C = ' num2str(C) '$ F, ' ...
         '$R = ' Rstr ', ' ...
         '$v_{in} = ' num2str(vin) '$ V, ' ...
         '$f = ' num2str(f) '$ Hz, ' ...
         '$w = ' num2str(w*100) '$\%, ' ...
         '$v_{in}/(1-w) = ' num2str(vin/(1-w)) '$ V)\\' ...
         '(Average($v_{out}$) = $' num2str(mean(sysstep(ceil(end/2):end,2))) '$ V, ' ...
         'Average($i_{out}$) = $' num2str(mean(sysstep(ceil(end/2):end,2))/R) '$ A, ' ...
         'Average($i_{in}$) = $' num2str(mean(sysstep(ceil(end/2):end,1))) '$ A)\\' ...
         '(Average($P_{out}$) = $' num2str(mean(sysstep(ceil(end/2):end,2))^2/R) '$ W, ' ...
         'Average($P_{in}$) = $' num2str(mean(sysstep(ceil(end/2):end,1))*vin) '$ W)}' ], ...
     'Interpreter', 'latex', 'FontSize', 14, 'BackgroundColor', 'white' );
title( 'DC-to-DC Switcher', 'Interpreter', 'latex', 'FontSize', 14 );
xlabel('Time (s)');
grid on;
set(gca, 'XTickLabel', []);
set(gca, 'XTick', 0:1/(2*f):tfinal);

    %% This function simulates the DC switcher
    function [nextstate] = switcher(t,y,L,C,vin,f,w)

        % y is [ i_L v_C ]

        if( xor( (t - floor(t*f)/f)*f >= w, control0 ) )
            % switch closed -- drive inductor current up

            nextstate = [ (1/L)*vin ; (1/C)*( 0 - y(2)/R ) ];

        else
            % switch open -- standard LC circuit

            % Only allow positive current.
            if( y(1) > 0 || ( y(1) == 0 && vin > y(2) ) )
                nextstate = [ (1/L)*(vin-y(2)) ;
                              (1/C)*( y(1) - y(2)/R ) ];
            else
                nextstate = [ 0 ; (1/C)*( y(1) - y(2)/R ) ];
            end

        end

    end

end