function [x,y]=program3_boundary(bs,s)
%PROGRAM3_GEOMETRY	Gives geometry data for the program3_geometry PDE model.
%
%   NE=PROGRAM3_GEOMETRY gives the number of boundary segments
%
%   D=PROGRAM3_GEOMETRY(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 regions.
%   Row 4 contains the number of the right-hand regions.
%
%   [X,Y]=PROGRAM3_GEOMETRY(BS,S) gives coordinates of boundary points. BS specifies the
%   boundary segments and S the corresponding parameter values. BS may be
%   a scalar.

nbs=4;

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

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

bs1=bs(:)';

if find(bs1<1 | bs1>nbs),
  error('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('bs must be scalar or of same size as s');
end

if ~isempty(s),

% boundary segment 1
ii=find(bs==1);
if length(ii)
t=[1 -0
0 1];
nth=100;
th=linspace(-3.1415926535897931,-1.5707963267948966,nth);
rr=t*[cos(th);sin(th)];
theta=pdearcl(th,rr,s(ii),d(1,1),d(2,1));
rr=t*[cos(theta);sin(theta)];
x(ii)=rr(1,:)+(0);
y(ii)=rr(2,:)+(0);
end

% boundary segment 2
ii=find(bs==2);
if length(ii)
t=[1 -0
0 1];
nth=100;
th=linspace(-1.5707963267948966,0,nth);
rr=t*[cos(th);sin(th)];
theta=pdearcl(th,rr,s(ii),d(1,2),d(2,2));
rr=t*[cos(theta);sin(theta)];
x(ii)=rr(1,:)+(0);
y(ii)=rr(2,:)+(0);
end

% boundary segment 3
ii=find(bs==3);
if length(ii)
t=[1 -0
0 1];
nth=100;
th=linspace(0,1.5707963267948966,nth);
rr=t*[cos(th);sin(th)];
theta=pdearcl(th,rr,s(ii),d(1,3),d(2,3));
rr=t*[cos(theta);sin(theta)];
x(ii)=rr(1,:)+(0);
y(ii)=rr(2,:)+(0);
end

% boundary segment 4
ii=find(bs==4);
if length(ii)
t=[1 -0
0 1];
nth=100;
th=linspace(-4.7123889803846897,-3.1415926535897931,nth);
rr=t*[cos(th);sin(th)];
theta=pdearcl(th,rr,s(ii),d(1,4),d(2,4));
rr=t*[cos(theta);sin(theta)];
x(ii)=rr(1,:)+(0);
y(ii)=rr(2,:)+(0);
end

end
