function [x,y] = geom(bs,s)

global depth width height thick gap space
% Geometry file definimg the geometry to simulate

% NE = GEOM gives the number of boundary segment
% D  = GEOM(bs) gives a matrix with one column for each boundary segment
%      specified in BS.
% Row 1 contains the start parameter value.
% Row 2 contains the end parameter value.
% Row 3 contains the number of the left hand region.
% Row 4 contains the number of the right hand region.
%
% [X,Y]=GEOM(BS,S) gives coordinates of boundary points. BS specifies the
%       boundary segments and S the corresponding parameter values. BS may
%       be a scalar.

nbs = 30;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Segments definitions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==0
    x=nbs; %number of boundary segment
    return
end

d=[
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  % start parameter value
  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  % end parameter value
  0 0 0 1 1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  % left hand region
  1 1 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  % right hand region
];

bs1 = bs(:)';

if find(bs1<1 | bs1>nbs),
  error('PDE:crackg:InvalidBs', 'Non existent boundary segment number.')
end

if nargin==1,
  x=d(:,bs1);
  return
end

x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
if m==1 && n==1,
  bs=bs*ones(size(s)); % expand bs
elseif m~=size(s,1) || n~=size(s,2),
  error('PDE:crackg:SizeBs', 'bs must be scalar or of same size as s.');
end

if ~isempty(s),
    
% boundary segment 1 - right top
ii=find(bs==1);
if length(ii)
x(ii)=interp1([d(1,1),d(2,1)],[depth/2 depth/2],s(ii));
y(ii)=interp1([d(1,1),d(2,1)],[gap+height thick+height],s(ii));
end

% boundary segment 2 - right bottom
ii=find(bs==2);
if length(ii)
x(ii)=interp1([d(1,1),d(2,1)],[depth/2 depth/2],s(ii));
y(ii)=interp1([d(1,1),d(2,1)],[height 0],s(ii));
end

% boundary segment 3 - bottom
ii=find(bs==3);
if length(ii)
x(ii)=interp1([d(1,2),d(2,2)],[depth/2 -depth/2],s(ii));
y(ii)=interp1([d(1,2),d(2,2)],[0 0],s(ii));
end

% boundary segment 4 - strip right
ii=find(bs==4);
if length(ii)
x(ii)=interp1([d(1,3),d(2,3)],[width/2 width/2],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[height+thick height],s(ii));
end

% boundary segment 5 - strip bottoml
ii=find(bs==5);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[0 -width/2],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 6 - strip bottomr
ii=find(bs==6);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[width/2 0],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 7 - strip topl
ii=find(bs==7);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[-width/2 0],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 8 - strip topr
ii=find(bs==8);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[0 width/2],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 9 - left top
ii=find(bs==9);
if length(ii)
x(ii)=interp1([d(1,6),d(2,6)],[-depth/2 -depth/2],s(ii));
y(ii)=interp1([d(1,6),d(2,6)],[thick+height gap+height],s(ii));
end

% boundary segment 10 -left bottom
ii=find(bs==10);
if length(ii)
x(ii)=interp1([d(1,6),d(2,6)],[-depth/2 -depth/2],s(ii));
y(ii)=interp1([d(1,6),d(2,6)],[0 height],s(ii));
end

% boundary segment 11 - top
ii=find(bs==11);
if length(ii)
x(ii)=interp1([d(1,7),d(2,7)],[-depth/2 depth/2],s(ii));
y(ii)=interp1([d(1,7),d(2,7)],[gap+height gap+height],s(ii));
end

% boundary segment 12 - strip left
ii=find(bs==12);
if length(ii)
x(ii)=interp1([d(1,8),d(2,8)],[-width/2 -width/2],s(ii));
y(ii)=interp1([d(1,8),d(2,8)],[height height+thick],s(ii));
end

% boundary segment 13 - stripr right
ii=find(bs==13);
if length(ii)
x(ii)=interp1([d(1,3),d(2,3)],[3*width/2+space 3*width/2+space],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[height+thick height],s(ii));
end

% boundary segment 14 - stripr bottomr
ii=find(bs==14);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[3*width/2+space 2*width/2+space],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 15 - stripr bottoml
ii=find(bs==15);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[2*width/2+space width/2+space],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 16 - stripr topr
ii=find(bs==16);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[2*width/2+space 3*width/2+space],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 17 - stripr topl
ii=find(bs==17);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[width/2+space 2*width/2+space],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 18 - stripr left
ii=find(bs==18);
if length(ii)
x(ii)=interp1([d(1,8),d(2,8)],[width/2+space width/2+space],s(ii));
y(ii)=interp1([d(1,8),d(2,8)],[height height+thick],s(ii));
end

% boundary segment 19 - stripl right
ii=find(bs==19);
if length(ii)
x(ii)=interp1([d(1,3),d(2,3)],[-width/2-space -width/2-space],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[height+thick height],s(ii));
end

% boundary segment 20 - stripl bottomr
ii=find(bs==20);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[-width/2-space -2*width/2-space],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 21 - stripl bottoml
ii=find(bs==21);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[-2*width/2-space -3*width/2-space],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 22 - stripl topr
ii=find(bs==22);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[-2*width/2-space -width/2-space],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 23 - stripl topl
ii=find(bs==23);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[-3*width/2-space -2*width/2-space],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 24 - stripl left
ii=find(bs==24);
if length(ii)
x(ii)=interp1([d(1,8),d(2,8)],[-3*width/2-space -3*width/2-space],s(ii));
y(ii)=interp1([d(1,8),d(2,8)],[height height+thick],s(ii));
end

% boundary segment 25 - shieldl right
ii=find(bs==25);
if length(ii)
x(ii)=interp1([d(1,3),d(2,3)],[-3*width/2-2*space -3*width/2-2*space],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[height+thick height],s(ii));
end

% boundary segment 26 - shieldl bottom
ii=find(bs==26);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[-3*width/2-2*space -depth/2],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 27 - shieldl top
ii=find(bs==27);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[-depth/2 -3*width/2-2*space],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 28 - shieldr bottom
ii=find(bs==28);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[depth/2 3*width/2+2*space],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[height height],s(ii));
end

% boundary segment 29 - shieldr top
ii=find(bs==29);
if length(ii)
x(ii)=interp1([d(1,5),d(2,5)],[3*width/2+2*space depth/2],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[height+thick height+thick],s(ii));
end

% boundary segment 30 - shieldr left
ii=find(bs==30);
if length(ii)
x(ii)=interp1([d(1,8),d(2,8)],[3*width/2+2*space 3*width/2+2*space],s(ii));
y(ii)=interp1([d(1,8),d(2,8)],[height height+thick],s(ii));
end

end