function [CT] = calculate_CT(Seconds, Throttle, RPM, b, h, gramsMeas, varargin)
% This function calculates the parameter for a zero intercept linear fit 
% between RPM^2 and thrust of the form:
%   Thrust = (RPM^2)*CT
% Output is given in the form of returned CT value and plots.
%
% Argument definitions:
%   b = length of test rig from pin joint to scale contact (same units as h)
%   h = length of test rig from pin joint to center motor axis (same units as b)
%   gramsMeas = user measured values direct from scale, MUST BE IN GRAMS!
%   These values must correspond to ALL unique throttle values, even those
%   outside of "LowCutoff" and "HighCutoff."
% (Note: values of gramsMeas corresponding to throttle setting outside 
% "LowCutoff" and "HighCutoff" are not used in calculation, therefore if
% measurements are only available for the range of "LowCutoff" to
% "HighCutoff", the vector can be padded with zeros in the appropriate
% locations.
%
% See also calculate_CR_B, calculate_CQ.

% Assign variable length input arguments
if ((nargin-5)==0) % if no input arguments are given for cutoff range:
LowCutoff = min(Throttle); % Use full range of input throttle data
HighCutoff = max(Throttle);
else
LowCutoff = varargin{1}; % Ignore data below this throttle value
HighCutoff = varargin{2}; % Ignore data above this throttle value
end

% Find unique Throttle Setting values (Not counting 0)
ThrottleU = unique(Throttle(Throttle~=0));

% Check length of inputs match
lengthGM = length(gramsMeas);
lengthTU = length(ThrottleU);
if lengthGM ~= lengthTU
    msgbox({['Length of measured grams values vector (' num2str(lengthGM) ')'];
        [' and number of unique Throttle values (' num2str(lengthTU) ')']; ...
       ' do not match! Check input data!'},'Argument Error','warn');
    return
end

% Find unique/filtered Throttle Setting values between low/high cutoff limits
ThrottleUF = ThrottleU(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff);

% Find coresponding gramsMeasurments
gramsMeasU = gramsMeas(ThrottleU>=LowCutoff & ThrottleU<=HighCutoff); % Unique scale values

% Calculate arms ratio
armRatio = b/h;
g = 9.81; % Gravity (m/s^2)

% Calculate thrust generated by motors
newtonsThrust = g*armRatio*gramsMeasU/1000; % Vector of thrusts in Newtons


% Find an average RPM value for each unique/filtered Throttle Setting
RPMAverages100 = zeros(length(ThrottleUF),1); % preallocate
for i = 1:length(ThrottleUF) % For all of the unique and filtered throttle settings:
    RPMUF = RPM(Throttle==ThrottleUF(i));  % get RPM values cooresponding to unique/filtered Throttle value
    if length(RPMUF)>=120 % check to make sure there are enough data points for filtering
        RPMfilt100 = RPMUF(100:end); % RPMUF minus the first 100 data points (to remove transient RPM fluctuation as steady state value is reached)
    else RPMfilt100 = RPMUF((fix(length(RPMUF)/2)):end); % If there's a limited number of points, just take the second half of the values as steady state  
    end
    RPMAverages100(i,1) = mean(RPMfilt100); % Take the average RPM value from these filtered data points
end

RPMsq = RPMAverages100.^2;
CT = RPMsq\newtonsThrust; % This is a technique for forcing a zero intercept using matrix multiplication techniques
CTp = [CT; 0];

% Plot the linear relationship and display calculated CT value
figure
plot(RPMsq,newtonsThrust,'*',RPMsq,polyval(CTp,RPMsq))
xlabel('RPM^2 (rev^2/min^2)','FontSize',12)
ylabel('Thrust (N)','FontSize',12)
title({['Linear Fit of Thrust vs. RPM^2']; ... 
    [' '];['CT: ' num2str(CT) ' (N/RPM^2)']},'FontSize',14)
xlims = get(gca,'xlim');
ylims = get(gca,'ylim');
grid on

% Plot Residuals for assesing general quality of fit
predictedThrust = polyval(CTp, RPMsq); % predict thrust for each throttle setting of interest
residuals = newtonsThrust - predictedThrust; % find absolute error between prediction and actual
figure
plot(predictedThrust,residuals,'r*',[min(predictedThrust) max(predictedThrust)],[0 0])
title('CT Residuals Plot','FontSize',14)
xlabel('Predicted Thrust (N)','FontSize',12)
ylabel('Prediction Error (N)','FontSize',12)
grid on
axis tight
end
