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

fprintf(1,'appel de parab,  nargin=%d\n',nargin);
if nargin>0, bs, end;
if nargin>1, s, end;
nbs=2;

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

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

bs1=bs(:)';
               bs1,
if find(bs1<1 | bs1>nbs),
  error('Non existent boundary segment number')
end

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

s,
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);
x(ii)=2*s(ii)-1;           y(ii)=( 2*s(ii)-1).^2;
disp('ii x y=');disp([ii,x(ii),y(ii)]);

% boundary segment 2
ii=find(bs==2);
x(ii)=1-2*s(ii);        y(ii)=ones(size(s(ii)));  % y=1

end
