function [fr, varargout] = feedcalc(ss, diameter, nflutes, varargin)
% FEEDCALC Calculates the cutting speed and spindle RPM for milling operations.
%
%       FR = FEEDCALC(SS, DIA, NFLUTES {,UNITS}) returns the range of feedrates for
%       machining a full-diameter slot, 0.25" (6.35 mm) deep in alloyed aluminum,
%       brass, bronze, and abrasive plastics using a HSS end mill. 
%
%       SS  - Spindle speed (RPM).
%       DIA - Mill diameter (in. or mm).
%       NFLUTES - Number of flutes on endmill.
%       UNITS   - (Optional.) Specifies working units 'in' or 'mm'. (Default is imperial.)
%       FR - Feed rate. Units are either in/min of mm/min depending on UNITS parameter.
%
%       [FR, SS] = FEEDCALC([], DIA, NFLUTES {,UNITS}) returns the ranges
%       of feedrates and spindle speed for a specified mill diameter and
%       number of flutes.
%
%       Results are suggested starting points only. They may be increased
%       or decreased depending on machine conditions, depth of cut,
%       required finish, coolant, etc.
%
%       Example:
%           y = feedcalc(1200, 10, 2, 'mm') calculates the feedrate (in mm/min) for a
%           10-mm, 2-fluted end mill rotating at 1200 RPM.

% -----------------------------------------------------------------------------------
%       Ross Wagner
%       McGill University, Dept. Biomedical Engineering
%       Ver 2.0
% -----------------------------------------------------------------------------------

SS_MAX = 4200;      % Max RPM of milling machine
SFM    = [100 200]; % Speed range for alloyed aluminum (surface-feet/minute)

% Setup chip load data for HSS mill in alloyed aluminum
% Table taken from the Melin Tool Co. catalogue (www.endmill.com)
%
% Column order: Mill dia. (in.), minimum chip load (in.), maximum chip load (in.)

HSSdata = [1/32, .0002 .0004; ...     % Travers Tool tech support recommends .0004 
           1/16, .0002 .0005; ...
           3/32, .0002 .0005; ...
           1/8 , .0002 .001 ; ...
           3/16, .0002 .001 ; ...
           1/4 , .0005 .002 ; ...
           5/16, .0005 .002 ; ...
           3/8,  .0005 .003 ; ...
           7/16, .0005 .003 ; ...
           1/2 , .0005 .003 ; ...
           9/16, .0005 .004 ; ...
           5/8 , .0005 .004 ; ...
          11/16, .0005 .004 ; ...
           3/4 , .001  .006 ; ...
          13/16, .001  .006 ; ...
           7/8 , .001  .006 ; ...
          15/16, .001  .006 ; ...
          1    , .002  .008 ; ...
           9/8 , .002  .012 ; ...
           5/4 , .002  .012 ; ...
          11/8 , .002  .012 ; ...
           3/2 , .003  .012 ; ...
          13/8 , .003  .012 ; ...
           7/4 , .003  .012 ; ...
          15/8 , .003  .012 ; ...
          2    , .003  .012 ; ...
          17/8 , .003  .012 ; ...
           9/4 , .003  .012 ; ...
          19/8 , .003  .012 ; ...
           5/2 , .003  .012 ; ...
           21/8, .003  .012 ; ...
           11/4, .003  .012 ; ...
           23/8, .003  .012 ; ...
           3   , .003  .012];

if ~isempty(ss)
  % Check spindle speed
  if ss > SS_MAX
    error(['Spindle speed exceeds milling machine capabilities of ' num2str(SS_MAX) ' RPM.'])
  end
end

% Process units argument and convert tool diameter to same units as table data.

dia_inch = diameter;      % Tool diameter in table units
convfact = 1;             % Conversion factor to get sizes into user preferred units
units = 'in.';
if length(varargin) > 0
  switch lower(varargin{1})
    case 'in'
      % Nothing to do. Use default values.
    case 'mm'
      dia_inch = diameter / 25.4;
      convfact = 25.4;
      units = 'mm';
    otherwise
      error('Not a valid units argument.')
  end
end

% Check tool diameter against table data
if dia_inch < HSSdata(1,1,1)
  str = ['Tool diameter must be larger than ' num2str(HSSdata(1,1,1)*convfact) ' ' units];
  error(str)
elseif dia_inch > HSSdata(end,1,1)
  dia_inch = HSSdata(end, 1, 1);
  str = ['Warning: Tool diameter beyond look-up table range. Using data for ' num2str(dia_inch*convfact) ' ' units];
  disp(str) 
end

%% Find index in table corresponding to largest tool diameter <= to diameter of
%% tool specified by user.
%idx = find(dia_inch >= HSSdata(:,1), 1, 'last');

% Find index in table corresponding to closest tool diameter specified by user.
[err, idx] = min(abs(HSSdata(:,1) - dia_inch));

% Find chip load range for this tool diameter
ChipLoad = HSSdata(idx, 2:3);

% Compute allowable range of spindle speeds for specified tool diameter
ss_allowed = 3.82 * SFM / dia_inch;   % Theoretically allowed tool RPM range
idx = find(ss_allowed > SS_MAX);      % Limit allowed RPM to within machine's abilities
for i = 1:length(idx)
  ss_allowed(idx(i)) = SS_MAX;
end

if isempty(ss)
  % User did not specify desired spindle RPM. Return the allowable range.
  varargout = {ss_allowed};
  
  % Calculate min and max feed rate based on extreme values of chip load and
  % spindle speed ranges.
  fr = ss_allowed .* ChipLoad * nflutes * convfact;   % in./min or mm/min depending on user preference  
else
  % Check whether user-specified RPM is greater than theoretically allowed
  if ss > ss_allowed(end)
    disp(['Warning: You are exceeding the recommeded spindle RPM range of: [' num2str(ss_allowed(1)) ', ' num2str(ss_allowed(2)) ']' ])
  end
  
  % Advise user if they could be milling at a faster speed, within the milling machines' abilities.
  if any(ss < ss_allowed)
    disp('Suggestion: Material SFPM allows for faster milling')
    disp(['by increasing spindle speed up to ' num2str(ss_allowed(2)) ' RPM.' ])
  end  
  fr = ss * ChipLoad * nflutes * convfact;   % in./min or mm/min depending on user preference
end