function [x,y]=UnitSquare(bs,s)
%RINGG        Gives geometry data for the unit square PDE model
%
%   NE=UnitSquare gives the number of boundary segment
%
%   D=UnitSquare(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]=UnitSquare(BS,S) gives coordinates of boundary points. BS specifies the
%   boundary segments and S the corresponding parameter values. BS may be
%   a scalar.

% Copyright 2011 The MathWorks, Inc.
% $Revision: 1.8.4.1 $   $Date: 2011/12/05 04:27:58 $

nbs=4;

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

%the direction of the boundary is counterclockwise 
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
];
%0: the outside region

bs1=bs(:)';

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

if ~isempty(s),

    % boundary segment 1
    ii=find(bs==1);
    if length(ii)
        x(ii)= s(ii);
        y(ii)= 0;
    end

    % boundary segment 2
    ii=find(bs==2);
    if length(ii)
        x(ii)= 1;
        y(ii)= s(ii);
    end
    
    % boundary segment 3
    ii=find(bs==3);    
    if length(ii)
        x(ii)= 1-s(ii);
        y(ii)= 1;
    end
    
    % boundary segment 4
    ii=find(bs==4);    
    if length(ii)
        x(ii)= 0;
        y(ii)= 1-s(ii);
    end
end