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

clear all
close all
clc

fmeas = [500, 1000, 2000, 5000, 10000, 20000, 50000, 100000, 200000, 300000,500000, 1000000, 2000000, 3000000]
ndat = length(fmeas);  %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 = ['test2_' num2str(fmeas(j)) 'Hz.csv']

    datfile = fopen(filename);

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

    ind = 1;
    %for i=1:10000
    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;
    

    [B,A] = butter(2,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(1+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)
end

Rs = min(Zmag)
Rp = max(Zmag)-Rs
Z_M = Zmag(8)
w = 2*pi*fmeas(8)
C =sqrt(Rp.*(2*Rs+Rp)+Rs.^2-Z_M.^2)./(Rp.*w.*sqrt(Z_M.^2-Rs.^2))

Z_M = Zmag(end)
w = 2*pi*fmeas(end)
Ls = sqrt(Z_M^2-(Rs+Rp/(1+(Rp*C*w)^2)))/w + Rp^2*C/(1+(Rp*C*w)^2)

figure(1)
loglog(fmeas,Zmag,'-b^','LineWidth',1.5)
xlabel('Frequency (Hz)')
ylabel('Impedance Magnitude (\Omega)')

save('LRdat','Rs','Ls');
