% OscillatorI4
% This routine calculates phase diagrams for the parallel
% RLC circuit
% (a) without a driving voltage, (b) with a driving voltage

% Lin = inductance (in mH)
% Cap = capacitance (in pF)
% Res = resistance (in kOhms)

% V0a = initial voltage amplitude, part (a) (in volts)
% I0a = initial current amplitude, part (a) (in uA)

% FreqRat = ratio of driving to resonant frequency,
%           part (b)
% V0b = initial voltage amplitudes, part (b) (in volts)
% J0 = driving current amplitude, part (b) (in uA)

% Nosc = number of oscillations plotted
% Npoints = number of points in each period

% omega = natural oscillation frequency
% gamma = damping rate
% Imp = real impedance
% phoff = phase offset

% omr = resonant frequency
% gamRat = ratio of the damping and resonant frequencies
% omRat = ratio of the resonant and natural frequencies
% ADc = driving voltage amplitude;
% CDc = cosine (driving angle)
% SDc = sine (driving angle)
% VDc, VDs = in-phase and quadrature quadratures of
%            the driving voltages

% [Current,Voltage] = current and voltage, part (a)

% Part (b):
% [CurD,VolD] = driven current and voltage
% [CurT1,VolT1] = transient current and voltage 1
% [CurT2,VolT2] = transient current and voltage 2

close all

Lin = 1, Cap = 100, Res = 50, V0a = 1, 
FreqRat = 1.0, V0b = [0 1],  J0 = 10,
   % set initial values
   
L = Lin*1e-3; C = Cap*1e-12; R = Res*1e3;
   % convert to SI units
gamma = 1/(2*R*C);  omega=((1/(L*C)) - gamma^2)^(1/2);
   % compute the damping coefficient and frequency
gam = gamma/omega;  %compute the damping ratio
Imp = sqrt(L/C);  phoff = atan(gamma/omega);
   % calculate the impedance amplitude and phase offset
Ia = V0a/Imp;  I0a = 1e6*Ia; 
   % compute the initial current amplitude (part a) in uA

Jb = J0*1e-6  %convert to SI units
omr = 1/sqrt(L*C); gamRat = gamma/omr; 
omRat = omr/omega;
   % compute the resonant frequency, damping ratio,
   % and resonant to natural frequency ratio
   % (part b)
ADc = 1/((1-FreqRat^2)^2 + 4*FreqRat^2*gamRat^2)^(1/2);
CDc = 2*FreqRat*gamRat*ADc;  SDc = (1-FreqRat^2)*ADc;
VDc = FreqRat*Imp*Jb*ADc*CDc;
VDs = FreqRat*Imp*Jb*ADc*SDc;
   % in-phase and quadrature components of the driven
   % voltage
Vc1 = V0b(1) - VDc;  Vc2 = V0b(2) - VDc;
Vs1 = omRat*VDs - gam*Vc1;
Vs2 = omRat*VDs - gam*Vc2;
  % transient voltage and current amplitude coefficients

Nosc = 15, Npoints = 100, Ntot = Nosc*Npoints + 1;
   % set the plotting parameters
   
Nt = 1:Ntot; time = 2*pi*(Nt-1)/Npoints;
   % set up the time array

% Part (a):  Compute and plot voltage and current   
Voltage = V0a*cos(time).*exp(-gam*time);
Current = -I0a*sin(time + phoff).*exp(-gam*time);
   % voltage and current as functions of time
   
plot(.001*Current,Voltage,'linewidth',1,'color','r')
% axis([-300,300,-1.0,1.0])
% title('Phase plot:  Voltage vs. Current')
% xlabel('Current (\muA)')
% ylabel ('Voltage (V)')
annotation('arrow','Position',[0.691,0.5,0,0.005])
set(gca,'fontname','Times New Roman','fontsize',20);
set(gca,'LineWidth',1.5,'TickLength',[0.025 0.025]);
title('Phase plot:  Voltage vs. Current','fontname','Times New Roman','fontsize',20)
xlabel('Current (mA)','fontname','Times New Roman','fontsize',20)
ylabel ('Voltage (V)','fontname','Times New Roman','fontsize',20)
xlim([-.500 .500]) %define limits for x axis
ylim([-1 1]) %define limits for y axis
set(gcf, 'WindowStyle', 'normal');
set(gca, 'Unit', 'inches');
set(gca, 'TickLabelInterpreter','LaTeX',...
    'YTick',[-1:.5:1],'YTickLabel',{'$-1$','$-0.5$','$0$','$0.5$','$1$'})

set(gca, 'TickLabelInterpreter','LaTeX',...
    'XTick',[-.5:.25:.5],'XTickLabel',{'$-0.5$','$-0.25$','$0$','$0.25$','$0.5$'}) %define the ticks and use LaTeX interpreter to have real minus signs
set(gcf, 'PaperPosition',1.9*[0 0 4 3],'PaperSize',1.9*[4 3]) %define the paper size and position
saveas(gcf,'3.pdf','pdf') %save file as a .pdf file


% Part (b):  Compute and plot the voltage and current
% evolution, starting from the given initial conditions.
% The initial current is assumed to equal zero.

FreqD = FreqRat*omRat;   % driving frequency
VolD = VDc*cos(FreqD*time) - VDs*sin(FreqD*time);
CurD = -(1/(omRat*Imp))*(VDs*cos(FreqD*time)...
    + VDc*sin(FreqD*time)); CurD = CurD*1e6
VolT1 = (Vc1*cos(time) - Vs1*sin(time)).*exp(-gam*time);
CurT1 = -(1/(omRat*Imp))*((Vc1-gamRat*Vs1)*sin(time)...
    +(Vs1+gamRat*Vc1)*cos(time)).*exp(-gam*time);
CurT1 = CurT1*1e6
VolT2 = (Vc2*cos(time) - Vs2*sin(time)).*exp(-gam*time);
CurT2 = -(1/(omRat*Imp))*((Vc2-gamRat*Vs2)*sin(time)...
    +(Vs2+gamRat*Vc2)*cos(time)).*exp(-gam*time);
CurT2 = CurT2*1e6
VolT1 = VolT1 + VolD;  CurT1 = CurT1 + CurD;
VolT2 = VolT2 + VolD;  CurT2 = CurT2 + CurD;
close all
plot(.001*CurD,VolD,'LineWidth',4)
hold on
plot(.001*CurT1,VolT1,'k',.001*CurT2,VolT2)
set(gca,'fontname','Times New Roman','fontsize',20);
set(gca,'LineWidth',1.5,'TickLength',[0.025 0.025]);
title('Phase plot:  Voltage vs. Current','fontname','Times New Roman','fontsize',20)
xlabel('Current (mA)','fontname','Times New Roman','fontsize',20)
ylabel ('Voltage (V)','fontname','Times New Roman','fontsize',20)
xlim([-.500 .500]) %define limits for x axis
ylim([-1 1]) %define limits for y axis
set(gcf, 'WindowStyle', 'normal');
set(gca, 'Unit', 'inches');
set(gca, 'TickLabelInterpreter','LaTeX',...
    'YTick',[-1:.5:1],'YTickLabel',{'$-1$','$-0.5$','$0$','$0.5$','$1$'})

set(gca, 'TickLabelInterpreter','LaTeX',...
    'XTick',[-.5:.25:.5],'XTickLabel',{'$-0.5$','$-0.25$','$0$','$0.25$','$0.5$'}) %define the ticks and use LaTeX interpreter to have real minus signs
set(gcf, 'PaperPosition',1.9*[0 0 4 3],'PaperSize',1.9*[4 3]) %define the paper size and position
saveas(gcf,'4.pdf','pdf') %save file as a .pdf file





