%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CSE Your Homework Assignment Project 8    %
%  Elastoplastic Torsion: Twist and Stress               %
%  problem3.m Dianne P. O'Leary   04/04                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Validate the dist_to_ellipse function
%  by plotting the distances on a grid of points in
%  the ellipse.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables:
% myalpha, mybeta = parameters defining the ellipse
%    (x/myalpha)^2 + (y/mybeta)^2 = 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myalpha mybeta

myalpha = 2;
mybeta = 1;

grid = -2:.1:2;

m = length(grid);

for i=1:m,
for j=1:m,
   if (grid(i)/myalpha)^2 + (grid(j)/mybeta)^2 >= 1
     d(i,j) = 0;
   else
     d(i,j) = dist_to_ellipse([grid(i);grid(j)],10*eps);
   end
end
end

contour(grid,grid,d)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test a point outside the ellipse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(sprintf('The distance from (5,1) is computed to be %f',...
      dist_to_ellipse([5,1])))