function y = geometry(x, boundary)
%GEOMETRY gives coordinates for the geometry
%   y = geometry(x, boundary) gives the y coordinate corresponding to the
%   input x coordinate at the boundary specified in the parameter boundary.
%   Valid inputs for boundary are
%       'lower_bound' - gives the lower boundary for the structure
%       'upper_bound' - gives the upper boundary for the structure
%       'lower_bound_hole1' - gives the lower boundary for the first hole
%       'upper_bound_hole1' - gives the upper boundary for the first hole
%       'lower_bound_hole2' - gives the lower boundary for the second hole
%       'upper_bound_hole2' - gives the upper boundary for the second hole
%       'plot' - plots the structure
%   Valid input values for x are in the interval 0 <= x <= 8 for the outer
%   boundary of the structure, 1.95 <= x <= 3.05 for the first hole and
%   3.8 <= x <= 4.6 for the second hole.

%Updated 2014-04-15 to correct rounding error


xp1=1;
xp2=0;
xp3=1.5;
xp4=8;
xp5=7;
yp2=2.5;
yp4=0;
r=1.2;

rh1=0.55;
rh2=0.4;
xh1=2.5;
xh2=4.2;

h1l=1.95;
h1r=3.05;
h2l=3.8;
h2r=4.6;


cx=r;
cy=yp2;

A=[xp1^2 xp1 1   0      0   0   0       0   0
   2*xp1 1   0   0      0   0   0       0   0
   0     0   0   xp3^2  xp3 1   0       0   0
   0     0   0   2*xp3  1   0   0       0   0
   xp4^2 xp4 1   0      0   0   0       0   0
   0     0   0   xp4^2  xp4 1   0       0   0
   0     0   0   0      0   0   cx^2    cx  1
   xp4   0.5 0   xp4    0.5 0  -2*xp4  -1   0
   0     0   0   0      0   0   xp4^2   xp4 1];

B=[sqrt(r^2-(xp1-cx)^2)+cy
   (cx-xp1)/sqrt(r^2-(xp1-cx)^2)
   -sqrt(r^2-(xp3-cx)^2)+cy
   (xp3-cx)/sqrt(r^2-(xp3-cx)^2)
   yp4
   yp4
   cy
   0
   yp4];
   
X=A\B;
a1=X(1);
b1=X(2);
c1=X(3);
a2=X(4);
b2=X(5);
c2=X(6);
a3=X(7);
b3=X(8);
c3=X(9);

if (strcmp(boundary, 'lower_bound') || strcmp(boundary, 'upper_bound'))
    if (x==xp2)
        y=yp2;
    elseif (x==xp4)
        y=yp4;
    elseif ((x>xp2) && (x<xp4))
        if strcmp(boundary, 'upper_bound') 
            if (x<=xp1)
                y=sqrt(r^2-(x-cx)^2)+cy;
            else
                y=a1*x.^2+b1*x+c1;
            end
        else
            if (x<=xp3)
                y=-sqrt(r^2-(x-cx)^2)+cy;
            else
                y=a2*x.^2+b2*x+c2;
            end
        end
    else
        error('Invalid x, input must be in the interval 0 <= x <= 8 for the outer boundary of the structure.');
    end
elseif (strcmp(boundary, 'lower_bound_hole1') || strcmp(boundary, 'upper_bound_hole1'))
    yh1=a3*xh1^2+b3*xh1+c3;
    if (x==h1l) || (x==h1r)
        y=yh1;
    elseif ((x>h1l) && (x<h1r))
        if strcmp(boundary, 'upper_bound_hole1') 
            y=sqrt(rh1^2-(x-xh1)^2)+yh1;
        else
            y=-sqrt(rh1^2-(x-xh1)^2)+yh1;
        end
    else
        error('Invalid x, input must be in the interval 1.95 <= x <= 3.05 for the first hole.');
    end
elseif (strcmp(boundary, 'lower_bound_hole2') || strcmp(boundary, 'upper_bound_hole2'))
    yh2=a3*xh2^2+b3*xh2+c3;
    if (x==h2l) || (x==h2r)
        y=yh2;
    elseif ((x>h2l) && (x<h2r))
        if strcmp(boundary, 'upper_bound_hole2') 
            y=sqrt(rh2^2-(x-xh2)^2)+yh2;
        else
            y=-sqrt(rh2^2-(x-xh2)^2)+yh2;
        end
    else
        error('Invalid x, input must be in the interval 3.8 <= x <= 4.6 for the second hole.');
    end
elseif strcmp(boundary, 'plot') 

    figure
    axis equal
    hold on

    x1=xp1:0.1:xp4;
    y1=a1*x1.^2+b1*x1+c1;
    plot(x1,y1,'b-')

    x2=xp3:0.1:xp4;
    y2=a2*x2.^2+b2*x2+c2;
    plot(x2,y2,'b-')

    t=(acos((xp1-cx)/r):0.01:(-acos((xp3-cx)/r)+2*pi));
    x3=r*cos(t)+cx;
    y3=r*sin(t)+cy;
    plot(x3,y3,'b-')

    tf=pi*(0:0.01:2);
    x5=rh1*cos(tf)+xh1;
    y5=rh1*sin(tf)+a3*xh1^2+b3*xh1+c3;
    plot(x5,y5,'c-')
    x6=rh2*cos(tf)+xh2;
    y6=rh2*sin(tf)+a3*xh2^2+b3*xh2+c3;
    plot(x6,y6,'c-')
    hold off
else
        error(['Invalid input for the parameter boundary. Valid inputs for boundary are ' ...
'                   ''lower_bound'' - gives the lower boundary for the structure '...
'                   ''upper_bound'' - gives the upper boundary for the structure '...
'                   ''lower_bound_hole1'' - gives the lower boundary for the first hole '...
'                   ''upper_bound_hole1'' - gives the upper boundary for the first hole '...
'                   ''lower_bound_hole2'' - gives the lower boundary for the second hole '...
'                   ''upper_bound_hole2'' - gives the upper boundary for the second hole '...
'                   ''plot'' - plots the structure']);
end

end

