function [xtp,xbp,theta_top,theta_bot]=UWpinch(alfa,zf,qf,R,V)
%function [xtp,xbp]=UWpinch(alfa,zf,qf,R,V)
%
%Based on Underwoods (JIP 1946)  equation 6
%
% Author(s): I. Halvorsen 999999
%
%

%
% Note: This function is a prototype under development.
%       The general problem is to decide the valid pinch
%       calculations.
%
%       Underwoods formula gets into trouble since some roots approach 
%       alfa-values
%
%       Note the simplification (by manipulating UW original)
%       use :  
%       
%
%       Observation: The pinch zone composition in one section is constant
%                    in one distribution zone !!!!!!
%
%                    E.g operating in zone AB/ABC, the pinch in the top
%                    is contant, but the feed pinch in the bottom is
%                    varying
%

nc=size(alfa,2);% Number of components
nt=size(alfa,1);% Number of data sets

%(Adjust recoveries a little to avoid zero and one)
eps=1.0e-6;
for i=1:length(R)
  if(R(i)==0)
    R(i)=R(i)+eps;
  elseif(R(i)==1)
    R(i)=R(i)-eps;
  end
end



%% For the top section:
%disp('Top')
w=zf.*R;
theta_top=UWallroots(alfa,w,V);
D=sum(w);

%Underwoods original expression:
abcd=w/V;
L=V-D;
m=L/D/(L/D+1);

for(k=1:nt),
  x=zeros(1,nc);
  for j=1:nc,	
    %Underwoods original expression:
    %for i=1:nc,	
    %  x(i)=abcd(k,i)*theta_top(k,j)/(m*(alfa(k,i)-theta_top(k,j)));
    %end
    
    %Simplified expression (vectorized):
    xtp(j,:)=w/L.*(theta_top(k,j)./(alfa(k,:)-theta_top(k,j)));
    
    %disp([theta_top(k,j),xtp(j,:),sum(xtp(j,:))]);
  end	

  II=find(~isfinite(xtp(:)));
  xtp(II)=0;
  II=find(abs(xtp(:))>10);
  xtp(II)=0;
end





%% For the bottom section:
%disp('Bottom')

%
%
  

w=zf.*(R-1);
V=V-(1-qf);

theta_bot=UWallroots(alfa,w,V);

abcd=w/V;
B=-sum(w);
L=V+B;
m=V/B/(V/B+1);


for(k=1:nt),
  x=zeros(1,nc);
  for j=1:nc,	
    %for i=1:nc,	
    %  x(i)=abcd(k,i)*theta_bot(k,j)/(m*(alfa(k,i)-theta_bot(k,j)));
    %end
    xbp(j,:)=w/L.*(theta_bot(k,j)./(alfa(k,:)-theta_bot(k,j)));
%    disp([theta_bot(k,j),xbp(j,:),sum(xbp(j,:))]);
  end	

  II=find(~isfinite(xbp(:)));
  xbp(II)=0;
  II=find(abs(xbp(:))>10);
  xbp(II)=0;

end



%%







