% LC_diode_sim (see script to change parameters)
%
%   Simulates a circuit:
%
%     vin --- L --- Diode >| --- C --- ground
%
%   In particular, finds inductor current and
%   capacitor voltage trajectories. The circuit can be visualized as
%   a bell that vibrates when struck with a hammer but can only expand in
%   size. Compare to an LC resonator, which can expand and contract (i.e.,
%   "ring").

% 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_sim

L = 1; % 1 Henry
C = 1; % 1 Farad
vin = 5; % 5 Volts

tfinal = 10; % 10 seconds

i0 = 0; % initial inductor current
v0 = 0; % initial capacitor voltage

[tout,sysstep] = ode45( @(t,y) switchoff(t,y,L,C,vin), [0 tfinal], [i0;v0] );

plot(tout, sysstep(:,1), '--', tout, sysstep(:,2), '-');
lfig = legend( '$i_L$ (A)', '$v_C$ (V)', 'Location', 'NorthWest' );
set(lfig, 'Interpreter', 'latex', 'FontSize', 14);
title( [ '\shortstack{LC + diode summary\\(' ...
         '$L = ' num2str(L) '$ H, ' ...
         '$C = ' num2str(C) '$ F, ' ...
         '$v_{in} = ' num2str(vin) '$ V)}' ], ...
      'Interpreter', 'latex', 'FontSize', 14 );
xlabel('Time (s)');
ylim(ylim.*[0 1] + [-0.5 0.5]);

    function [nextstate] = switchoff(t,y,L,C,vin)

        % y is [ i_L v_C ]

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

    end

end