...

function [No_pts_clip,x_clip,y_clip]= ...
         genarea(No_pts_orig,x_orig, ...
         y_orig,y_clip_plane)
%
% [No_pts_clip,x_clip,y_clip]= ...
%        genarea(No_pts_orig,x_orig, ...
%        y_orig,y_clip_plane)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function generates a new polygon
% by clipping the reference polygon at
% the prescribed horizontal plane.
%
% No_pts_orig - number of coordinates defining
%               the reference polygon
% x_orig      - vector of x coordinates in
%               the reference polygon
% y_orig      - vector of y coordinates in
%               the reference polygon
% y_clip_plane- the y coordinate of the 
%               horizontal clipping plane 
%
% No_pts_clip - number of coordinates defining
%               the clipped polygon
% x_clip      - vector of x coordinates in
%               the clipped polygon
% y_clip      - vector of y coordinates in
%               the clipped polygon
%
% User m functions called:  clip
%----------------------------------------------

Epsilon=1e-8;

%...Note: Define the clipping plane using the
%...      normal form:
%...      
%...      alpha   beta*x   gamma*y
%...      ----- + ------ + ------- = 0
%...        d        d        d    
%...      
%...      where:
%...      
%...      d=SIGN(beta)*SQRT(beta^2+gamma^2)
%...
Nform(1)=y_clip_plane; % perpendicular distance
Nform(2)=0;            % run of normal vector
Nform(3)=-1;           % rise of normal vector

x_clip=x_orig; y_clip=y_orig;
NoClipPts=No_pts_orig;

%...Clip the polygon
[NoNewPts,IsetClipPoly,NoInterpPts, ...
  x_clip,y_clip]=clip(x_clip,y_clip,NoClipPts, ...
  Nform,Epsilon);

%...Clip complete, reorder arrays
Loop=NoClipPts+NoInterpPts;
xtmp=x_clip; ytmp=y_clip;
clear x_clip y_clip;
for j=1:NoNewPts
  itmp=IsetClipPoly(j);
  x_clip(j)=xtmp(itmp); y_clip(j)=ytmp(itmp);
end
No_pts_clip=NoNewPts;