% Polygon and Inverse Distance Squared Interpolation
% Example of .DAT file
% 3
% 12.2    1.34     16.34     29.1
% 10.1    2.10     15.34     22.34
% 9.21    4.34     12.45     19.23
% Must have four columns - they are, in order, x,y,z position and magnitude
%
% Description:
%   For your information: this is a functional solution to the interpolation 
%   problem whereby a user-defined interpolation function is passed as an
%   argument and used to compute values at points v(i,j,k).  This has the
%   advantage of eliminating the "if" to test for type of interpolation.
%   doing so makes this code useful for any unforseen interpolation methods -
%   in such cases a user merely must write an interpolation function and
%   pass it in.
%
% Inputs:
%   rn = minimum of range
%   rx = maximum of range
%   n = number of slices
%   filename = name of file
%   interp_func = the interpolation function
% Output:
%   A 3D rendering of an isosurface based on the "interp_func" interpolation
%   method.
%
function interp_f(filename, rn, rx, n, interp_func);
   fid = fopen(filename,'r');       % open the data file
   if fid < 0
      error(['File ' filename ' not found']);
   end
   npts = fscanf(fid,'%d',1);       % get the number of data points
   dpts = zeros(npts,4);            % create array to hold the data

   for i = 1:npts                   % read in coordinates and values
      dpts(i,:) = fscanf(fid,'%f',4);
   end

   fclose(fid);

   % Create the volume matrix which holds the processed data and prepares it
   % for viewing
   x = linspace(rn,rx,n);
   y = linspace(rn,rx,n);
   z = linspace(rn,rx,n);
   v = interp_func(x,y,z,dpts);

   % View the interpolated data
   hpatch = patch(isosurface(x,y,z,v,23));
   isonormals(x,y,z,v,hpatch);
   set(hpatch,'FaceColor','green','EdgeColor','none');
   camlight right;
   lighting phong
   view(45,45);

   % Pretty up the figure
   grid on;
   xlabel('x');
   ylabel('y');
   zlabel('z');
   title('Interpolation');
   axis([rn rx rn rx rn rx]);
   axis equal;
end