% Data processing for Test 3 CSV Scope data from Tektronix MSO 4034
% Outputs Rp and Cp values as a .mat file
% Katherine A. Kim
% University of Illinois at Urbana-Champaign
% 1/31/2013

clear all
close all
clc

ndat = 14;  %Number of data files

for j = 1:ndat % Number of data files
    T   = zeros(1,10000);
    CH1 = zeros(1,10000);
    CH2 = zeros(1,10000);

    filename = ['test3_500hz_' num2str(j) '.csv']
    datfile = fopen(filename);

    for k = 1:15    %ignore the first 19 lines
        header = fgetl(datfile);
    end
    header
    devline = fgetl(datfile);

    ind = 1;
    while devline~=-1         
%         disp(devline)
        X = textscan(devline, '%f %f %f','delimiter',',');
        if ~isempty(X{3})
            T(ind) = X{1};
            CH1(ind) = X{2};
            CH2(ind) = X{3};
        end
        devline = fgetl(datfile);
        ind = ind+1;
    end
    V = CH1;
    I = CH2;
    
    %Filter Waveform
    [B,A] = butter(3,0.01);

    Vfil = filter(B,A,V);
    Ifil = filter(B,A,I);
    
    n = length(Vfil);

    start = 300;
    Tfil = T(start:n);
    Vfil = Vfil(start:n);
    Ifil = Ifil(start:n);
    
    figure(3+j)
    subplot(211)
    plot(T,V,Tfil,Vfil,'r')
    subplot(212)
    plot(T,I,Tfil,Ifil,'r')

    Vamp(j) = max(Vfil)-min(Vfil);
    Iamp(j) = max(Ifil)-min(Ifil);
    Zmag(j) = Vamp(j)/Iamp(j);
    
    a = fft(Vfil);
    a = a(2:end);
%     a(abs(a) < 1e-10) = 0;
    b = fft(Ifil);
    b = b(2:end);
%     b(abs(b) < 1e-10) = 0;

    ph_sh = abs(angle(a)-angle(b));
    PS(j) = median(ph_sh);
end

load 'LRdat' % Rs, Ls
Rp = Zmag-Rs;

Idat = load('I500Hz.dat');
Vdat = load('V500Hz.dat');

figure(1)
subplot(211)
plot(Vdat(2:end),-Idat(2:end),'^-k','LineWidth',1.5)
xlabel('Voltage (V)')
ylabel('Current (A)')
grid on
title('Current vs. Voltage')

subplot(212)
semilogy(Vdat(2:end),Zmag(2:end),'o-','LineWidth',1.5)
xlabel('Voltage (V)')
ylabel('50 Hz Impedance ( \Omega)')
grid on
title('Low Freq Impedance vs. Voltage')

%% 100 kHz

for j = 1:ndat % Number of data files
    T   = zeros(1,10000);
    CH1 = zeros(1,10000);
    CH2 = zeros(1,10000);

    filename = ['test3_100khz_' num2str(j) '.csv']
    datfile = fopen(filename);

    for k = 1:15    %ignore the first 19 lines
        header = fgetl(datfile);
    end
    header
    devline = fgetl(datfile);

    ind = 1;
    while devline~=-1         
%         disp(devline)
        X = textscan(devline, '%f %f %f','delimiter',',');
        if ~isempty(X{3})
            T(ind) = X{1};
            CH1(ind) = X{2};
            CH2(ind) = X{3};
        end
        devline = fgetl(datfile);
        ind = ind+1;
    end
    V = CH1;
    I = CH2;
    
    %Filter Waveform
    [B,A] = butter(3,0.01);

    Vfil = filter(B,A,V);
    Ifil = filter(B,A,I);
    
    n = length(Vfil);

    start = 300;
    Tfil = T(start:n);
    Vfil = Vfil(start:n);
    Ifil = Ifil(start:n);
    
    figure(3+ndat+j)
    subplot(211)
    plot(T,V,Tfil,Vfil,'r')
    subplot(212)
    plot(T,I,Tfil,Ifil,'r')

    Vamp(j) = max(Vfil)-min(Vfil);
    Iamp(j) = max(Ifil)-min(Ifil);
    Zmag(j) = Vamp(j)/Iamp(j);
    
    a = fft(Vfil);
    a = a(2:end);
%     a(abs(a) < 1e-10) = 0;
    b = fft(Ifil);
    b = b(2:end);
%     b(abs(b) < 1e-10) = 0;

    ph_sh = abs(angle(a)-angle(b));
    PS(j) = median(ph_sh);
end

Idat = load('I100kHz.dat');
Vdat = load('V100kHz.dat');

% figure(3)
% plot(Vdat,Zmag,'--o')
% figure(4)
% plot(VdatC,Zmag,'--o')

%% Calculate Capacitance

f = 100e3;
w = 2*pi*f;
Z_M = Zmag;
C =sqrt(Rp.*(2*Rs+Rp)+Rs.^2-Z_M.^2)./(Rp.*w.*sqrt(Z_M.^2-Rs.^2));

VdatC = Vdat;
IdatC = Idat;
i = 1;
while i<=size(C,2) %Filters out bad measurements
    if (isreal(C(i))==0)
        if i==1
            C = C(2:end);
            VdatC = VdatC(2:end);
            IdatC = IdatC(2:end);
        elseif i == size(Z_M,2)
            C = C(1:end-1);
            VdatC = VdatC(1:end-1);
            IdatC = IdatC(1:end-1);
        else
            C = [C(1:i-1) C(i+1:end)];
            VdatC = [VdatC(1:i-1) VdatC(i+1:end)];
            IdatC = [IdatC(1:i-1) IdatC(i+1:end)];
        end
    else
        i = i+1;
    end
end 

figure(2)
subplot(211)
semilogy(Vdat(2:end),Zmag(2:end),'o-b','LineWidth',1.5)
grid on
xlabel('Voltage (V)')
ylabel('50 kHz Impedance ( \Omega)')
title('High Freq Impedance vs. Voltage')

subplot(212)
semilogy(VdatC(2:end),C(2:end),'s-','LineWidth',1.5,'Color',[0 0.5 0])
grid on
xlabel('Voltage (V)')
ylabel('Capacitance (F)')
title('Capacitance vs. Voltage')

%Find Breakdown Voltage
[Zmax,ibd] = max(Zmag);
Vbd = Vdat(ibd)

save('Cdat.mat','C','VdatC','IdatC','Vbd');