%Prelab 1 Matlab
%This MATLAB file takes a set of frequency, peak to peak input/output, and
%phase measurements and matches them against a known transfer function to
%find the likely cutoff frequency, 1/RC. This frequency then can be used to
%find the unknown capacitance given the resistance is known.
%
%Written by Stephen Maloney
%

clear all; clc; close all;

%Input data from LABVIEW - given in the prelab
input = [120, 4.1202, 3.7137, -23.7380; 160, 4.0984, 3.6122, -27.6350; 200, ...
    4.0939, 3.4766, -32.5860; 240, 4.0995, 3.3042, -37.1420; 280, 4.1080, ...
    3.1847, -38.8100; 320, 4.0994, 3.0345, -41.8310; 360, 4.0994, 2.8817, ...
    -45.8470; 400, 4.0833, 2.7716, -47.4650; 440, 4.0875, 2.5484, -52.3820];

%The frequency range we are interested in for the transfer function
w = 2*pi*input(:, 1)';

%The range to try for the cutoff frequency
wc = 2*pi*linspace(300, 400, 100);

%Build up the magnitude plot - output/input
mag = (input(:, 3)./input(:, 2))';

%Build up the phase plot
phase = input(:, 4)';

mError = zeros(1, length(wc));
pError = zeros(1, length(wc));

for idx = 1:length(wc),
    %Calculate the frequency response given wc(idx) as the cutoff freq
    trial = 1./(1+j*w./(wc(idx)));
    
    %Find the mean square error for this generated response and the given
    %response. Notice that dividing by n is not necessary as every trial
    %would be divided by the same n
    mError(idx) = sum((abs(trial) - mag).^2);
    pError(idx) = sum((angle(trial)*180/pi - phase).^2);
    
    %For debug - this will show the phase fit in progress
    semilogx(w, angle(trial)*180/pi, w, phase);
    pause(.01);
end

%Display the output
[minVal, minIdx] = min(mError);
disp('Magnitude Match Cutoff frequency (Hz):');
FcM = wc(minIdx)/(2*pi)

[minVal, minIdx] = min(pError);
disp('Phase Match Cutoff frequency (Hz):');
FcP = wc(minIdx)/(2*pi)

%If the resistor value is known, the cutoff frequency can be used to find
%the unknown capacitance; Tau = RC = 1/Fc, C = 1/R*fc
R = 1e3;

disp('Magnitude Fit Unkown Capacitor Value (F):');
CM = 1/(R*FcM)

disp('Phase Fit Unkown Capacitor Value (F):');
CM = 1/(R*FcP)