% 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

filename = input('Data file: ','s');
fid = fopen(filename,'r');       % open the file
if fid < 0
   error(['File ' filename ' not found']);
end
npts = fscanf(fid,'%d',1);       % get the number of points
dpts = zeros(npts,4);            % array to hold the data

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

fclose(fid);

stype = input('Type of interpolation (1=Polygon, 2=Inv Dis Squ): ');
rn = input('Range minimum: ');
rx = input('Range maximum: ');
n = input('No. grid elements: ');

x = linspace(rn,rx,n);
y = linspace(rn,rx,n);
z = linspace(rn,rx,n);
v = zeros(n,n,n);

switch stype
   case {1}
      fprintf('The type is Polygon\n');
      title_str = 'Lab4 - Polygon interpolation';
   case {2}
      fprintf('The type is Inverse Distance Squared\n');
      title_str = 'Lab4 - Inverse Distance Squared Interpolation';
end

for i=1:n
   for j=1:n
      for k=1:n
         dx2 = (x(j)-dpts(:,1)).^2;
         dy2 = (y(i)-dpts(:,2)).^2;
         dz2 = (z(k)-dpts(:,3)).^2;
         if stype == 1                % Polygon interpolation
            [val nearest] = min(dx2 + dy2 + dz2);
            v(i,j,k) = dpts(nearest,4);
	 elseif stype == 2            % Inverse distance squared interpolation
            dp = 1./(eps + dx2 + dy2 + dz2);
            fp = dpts(:,4)./(eps + dx2 + dy2 + dz2);
            v(i,j,k) = sum(fp)/sum(dp);
         end
      end
   end
end

clf;
isosurface(x,y,z,v,23);
grid
xlabel('x');
ylabel('y');
zlabel('z');
title(title_str);
axis equal;
axis([rn rx rn rx rn rx]);
grid on;