clear all
close all
%Plot the waveforms for a real choke

alpha = 2*pi;
Vcc = 5;
P = 1;
L = 100e-6;
Imax = 4*P/Vcc;
Om = 2*pi*27e6;
N = 300000;
t = linspace(0,50000,N)*2e-9;

I1 = (Imax/2/pi)*(alpha - sin(alpha))/(1 - cos(alpha/2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%THE CRUCIAL BIT
%Vs = 2*P./I1;
Vs = Vcc; %% why limit ourselves to this??
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Vs_ratio = Vs/2/Vcc
Vq = 1 - Vs/2/Vcc
Rload = Vs/I1

RclassA = 2*Vcc/Imax;
RL = Rload / RclassA

gamma = Rload/L;
%Pdis = Vcc*Idc - Vs*I1/2;

Co = -gamma*Imax*(gamma/(Om^2 + gamma^2) - 1/gamma*cos(alpha/2))/(1 - cos(alpha/2));
C1 = 0;
Io = zeros(1,N);
Id = zeros(1,N);

% Solve for the supply current from the DE...given Id
for i=1:N

Id(i) = Imax*(cos(Om*t(i)) - cos(alpha/2))/(1 - cos(alpha/2));
if(Id(i) > 0)

     Io(i) = Co*exp(-gamma*t(i));

     Io(i) = Io(i) + Imax*gamma*((gamma*cos(Om*t(i))+Om*sin(Om*t(i)))/(Om^2 + gamma^2) - cos(alpha/2)/gamma)/(1 - cos(alpha/2));

     C1 = exp(gamma*t(i))*Io(i);
     else
     Id(i) = 0;
     Io(i) = C1*exp(-gamma*t(i));
     Co = Io(i) - Imax*gamma*((gamma*cos(Om*t(i))+Om*sin(Om*t(i)))/(Om^2 + gamma^2) - cos(alpha/2)/gamma)/(1 - cos(alpha/2));

     Co = Co*exp(gamma*t(i));

end
end

Iload = Io - Id;

figure(2)
subplot(2,1,1)
plot(t,Id,t,Io,'r',t,Iload,'k','LineWidth',3)
title('Currents')
axis([t(N/2) t(N/2*1.02) -1 1 ])

subplot(2,1,2)

%Need to filter out all but the fundamental for Iload
IloadF = fftshift(fft(Iload));
dt = t(2) - t(1);
freq = linspace(0,(N-1)/N/dt,N) - 1/dt/2.;
%semilogy(freq,abs(IloadF))
%axis([-2*Om/2/pi 2*Om/2/pi 0  max(abs(IloadF))])

jj = find(abs(freq) >= 1.5*Om/2/pi);

IloadF(jj)=0*IloadF(jj);

Iload_filtered = ifft(fftshift(IloadF));

ivals = (N/2):N;
plot(t,Iload,t,Iload_filtered)
axis([t(N/2) t(N/2*1.02) -1 1])

%The problem is here .. you need to filter out all but the fundamental
Vd = Vcc + Iload_filtered*Rload;

Comp = (max(Vd(ivals))- min(Vd(ivals)))/2/Vcc

plot(t,Vd,'k','LineWidth',3)
title('Drain Voltage')
axis([t(N/2) t(N/2*1.02) min(Vd) max(Vd)])

Pcc = Vcc*mean(Io(ivals))
Pload = mean(Iload_filtered(ivals).^2*Rload)
Pdis = mean(Vd(ivals).*Id(ivals))
Pload+Pdis
eta = Pload/Pcc

