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

global depth width height gap
% 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 = 6;

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

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

d=[
  0 0 0 0 0 0 % start parameter value
  1 1 1 1 1 1 % end parameter value
  1 0 0 1 1 1 % indice of left hand region
  0 1 1 0 1 1 % indice of 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
ii=find(bs==1);
if ~isempty(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/2 height/2],s(ii));
end

% boundary segment 2 - left
ii=find(bs==2);
if ~isempty(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)],[-height/2 height/2],s(ii));
end

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

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

% boundary segment 5 - strip top
ii=find(bs==5);
if ~isempty(ii)
x(ii)=interp1([d(1,5),d(2,5)],[-width/2 width/2],s(ii));
y(ii)=interp1([d(1,5),d(2,5)],[gap/2 gap/2],s(ii));
end

% boundary segment 6 - strip bottom
ii=find(bs==6);
if ~isempty(ii)
x(ii)=interp1([d(1,6),d(2,6)],[-width/2 width/2],s(ii));
y(ii)=interp1([d(1,6),d(2,6)],[-gap/2 -gap/2],s(ii));
end

end