% process and plot results from the lamp I&V measurements

% transfmr_L1_v
% transfmr_L1_i
% 
% mains_L1_v
% mains_L1_i
% 
% mains_L2_v
% mains_L2_i
% 
% mains_L3_v
% mains_L3_i


names = { ...
 'mains_L2'
 'mains_L3'  };
nl = size(names,1);

for c=1:nl,
    
    vdata = csvread([names{c},'_v.csv']);
    idata = csvread([names{c},'_i.csv']);
    t  = vdata(:,1);
    tt = idata(:,1);
    if (t~=tt),
        error('there seem to be different times for i&v measurements');
    end
    n = length(t);
    v = vdata(:,2) * 10;   % x10 probe
    i = idata(:,2) /  3.3;  % 3.3Ohm shunt
    
    figure
    plot(t, v, 'b', t, i*1000, 'r')
    title( [ 'v(t) & i(t) for   ', strrep(names{c},'_',' ') ] )
    xlabel('time')
    ylabel('v(t)/V  and  i(t)/mA')
    drawnow
    
    % assume integer number of periods was recorded each time
    % calculate powers and indices
    p = sum(v.*i)/n;
    vrms = sqrt(sum(v.^2)/n);
    irms = sqrt(sum(i.^2)/n);
    vmean = sum(abs(v))/n;
    imean = sum(abs(i))/n;
    vff = vrms/vmean;
    iff = irms/imean;
    % still assuming integer number of cycles!
    % do FFT
    ncyc = 2;  % how many cycles...
    nh = 10; % highest hamronic to plot
    vspec = fft(v)/n;
    ispec = fft(i)/n;
    vfd = [ vspec(1); 2*vspec(1+ncyc*(1:1:nh)) ];
    ifd = [ ispec(1); 2*ispec(1+ncyc*(1:1:nh)) ];
    % now we have fundamental quantities we can get reactive power
    s = 0.5 * vfd(2)*conj(ifd(2));
    q = imag(s);
    pff = real(s)/abs(s);  % fundamental freq PF
    pfh = p/abs(s);  % PF using (harmonic)actual power and fundamental current
    
    fprintf( [ '\nTest:  ', names{c}, ...
        '\n\tP    : ', num2str(p), 'W', '\tQ : ', num2str(q), 'VAr', ...
            '\t pf(fund): ', num2str(pff), '\t pf(harm): ', num2str(pfh),  ...
        '\n\tVrms : ', num2str(vrms), 'V', '\tVmean : ', num2str(vmean), 'V', ...
            '\tVpk/sqrt(2): ', num2str(max(abs(v))/sqrt(2)), '\t form: ', num2str(vff), ...
        '\n\tIrms : ', num2str(irms), 'A', '\tImean : ', num2str(imean), 'A', ...
            '\tIpk/sqrt(2): ', num2str(max(abs(i))/sqrt(2)), '\t form: ', num2str(iff), '\n', ...
        ])
    
    figure
    title( [ 'Current spectrum of   ', strrep(names{c},'_',' ') ] )
    subplot(2,1,1)
    bar( 0:nh, abs(ifd) );
    xlabel('harmonic order')
    ylabel('abs( i(h) ), Amps')
    subplot(2,1,2)
    bar( 0:nh, angle(ifd)*180/pi );
    xlabel('harmonic order')
    ylabel('arg( i(h) ), degrees')
    drawnow
    
end

disp(' ')