function [np,xs,ys]=ptop(R,p,q,xf,xd,xb,npx)
% ptop - plate to plate calculation
% function [np,xs,ys]=ptop(R,p,q,xf,xd,xb,npx)
% Returns:
%  np number of plates
%  xs, ys vectors to show steps up the column
% Parameters:
%  R Reflux ratio: L/D
%  p pressure kPa
%  q quality of feed: 1=sat vapor
%  xf, xd, xb mol fraction of 1st compound in feed, 
%     distillate, bottom
%  npx max number of plates
if nargin <7
   npx=50;
end
np=0;
[LoV,DxDoV,LpoVp,BxboVp,xo,yo]=oplines(R,q,xf,xd,xb);
xc=xb;
y=xb; xs=[]; ys=[];
while xc<xo(2)
   xs=[xs,xc];
   ys=[ys,y];
   [t,y]=bubpt2(p,xc);
   xs=[xs,xc];
   ys=[ys,y];
   xc=(y+BxboVp)/LpoVp;
   np=np+1;
   if np>npx
      disp('Too many plates in the bottom')
      return
   end
end
xc=(y-DxDoV)/LoV;
while xc<xd
   xs=[xs,xc];
   ys=[ys,y];
   [t,y]=bubpt2(p,xc);
   xs=[xs,xc];
   ys=[ys,y];
   xc=(y-DxDoV)/LoV;
   np=np+1;
   if np>npx
      disp('Too many plates in the top')
      return
   end
end
x2np=xs(2*np);
np=np+(x2np-xd)/(x2np-xc);
   
