close all; clear all

if(0)
    fname = '../inv_osc';
    X=readtable([fname, '.txt']);
    X.Properties.VariableNames


    t = X.time;
    f = X.V_f_;
    i2 = X.V_i2_;
    i = X.I_L1_;
    u_rect1 = X.V_ul_rect1_;
    u = X.V_ul1_;
    l_ref = X.V_l_var_/100;
    r_ref = X.V_r_var_*1000;
    r = X.V_p_avg__V_i2__1000;
    p = X.V_p_avg_;
    save('inv_osc.mat', 't', 'f', 'i2', 'i', 'u_rect1', 'u', 'r', 'p', 'l_ref', 'r_ref');
    clear all
end

load('inv_osc.mat');


%wartosci freq wyliczone przez LTspice na podstawie czasu 10 przejsc przez zero
f0 = 57748.6        %f na poczatku
f45 = 59031.4       %f na srodku (45 ms)

plot(t, l_ref); grid; xlabel('t[s]'); ylabel('L[H]');
plot(t, r_ref); grid; xlabel('t[s]'); ylabel('R[Ohm]');
plot(t, f); grid; xlabel('t[s]'); ylabel('f[-]');
plot(t, u_rect1); grid; xlabel('t[s]'); ylabel('mean(u-rect) [V]');
plot(t, i2); grid; xlabel('t[s]'); ylabel('mean(i2) [A]');
%plot(t, r); grid; xlabel('t[s]'); ylabel('mean(p) [W]');

%res 
figure;
plot(t, p./i2, t, r_ref); grid; xlabel('t[s]'); ylabel('r [Ohm]');


%frequency by zero-corssing
[fx, tx] = freq_zx(u, t);

plot(tx, fx); title('freq by zero-corssing'); xlabel('t[s]'); ylabel('fx [Hz]')

C1 = 7.5e-9; C2 = C1;
C = C1*C2/(C1+C2);
L0 = mean(l_ref)
f_eq0 = 1 ./ ( 2*pi*sqrt((L0+300*C1)*C) )
L = l_ref;
f_eq = 1 ./ ( 2*pi*sqrt((L+300*C1)*C1*C2/(C1+C2)) );

kf = f(1) / fx(1)
fe = f/kf; %estymata czestotliwosci (nie zgadza sie f45)


figure;
plot(t, f_eq, t, fe, tx, fx); grid; title('freq'); ylabel('freq[Hz]'); xlabel('t[s]');
legend('eq.', 'f/u', 'zx')

%ind by eq based on fx 
Le = 1./((2*pi*fx).^2*C);
kL = Le(1)/l_ref(1)
Le_corr = Le/kL;
figure;
subplot(3,1,1); plot(t, r_ref); grid; legend('r-ref');
subplot(3,1,2); plot(t, l_ref); grid; legend('l-ref');
subplot(3,1,3); plot(tx, Le, tx, Le_corr, t, l_ref); grid; title('ind'); 
ylabel('L[H]'); xlabel('t[s]'); 
legend('est', 'est-corr', 'ref')



function [fx, tx] = freq_zx(x, t)
%     idx = find( x(2:end).*x(1:end-1)<0 );
%     tx = zeros(size(idx));
%     for i=1:numel(idx)
%         j = idx(i); % Need both the index in the idx/t_zero vector, and the f/t vector
%         tx(i) = interp1( x(j:j+1), t(j:j+1), 0.0, 'linear' );
%     end
    
    i = find( x(2:end).*x(1:end-1)<0 );
    a=(x(i+1)-x(i))./(t(i+1)-t(i)); % slope
    t0=-x(i)./a+t(i);               % x coordinates of zero crossings linear interpolated

    tx = t0(1:2:end);
    fx = 2./diff(tx)/2;
    tx = tx(2:end);
end
    
