function [Vs,Ds,Rs,Keys,theta,VM,h]=UWmulti(theta,alfa,zf,qf,F,inkey,plotflag)
%UWmulti(theta,alfa,zf,qf)
%[Vs,Ds,Rs,Keys]=UWmulti(theta,alfa,zf,qf)
%[Vs,Ds,Rs,Keys]=UWmulti(theta,alfa,zf,qf,F)
%[Vs,Ds,Rs,Keys,theta,VM,h]=UWmulti(theta,alfa,zf,qf,F,inkey,plotflag)
%[Vs,Ds,Rs,Keys,theta,VM,h]=UWmulti(theta,alfa,zf,qf,F,[],plotflag)
%[Vs,Ds,Rs,Keys,theta,VM,h]=UWmulti([],alfa,zf,qf,1,[],1)
%
%Compute the minimum energy points of the multicomponent-mountains
%Every set V,D and R is computed as the point of minimum
%energy comsumption for a legal combination of LK and HK
%Vs=[V1;V2;...],Ds=[D1; ...],R=[R1;R2....]
%Keys=[LK1,HK1;LK2,HK2;.....]
%theta: underwood roots, as from theta=UWroots(alfa,zf,qf)
%VM      : equation matrix
%h       : figure handle(s) from plotting
%inkey   : If included in argument list, compute for Keys==inkeys
%plotflag: If plotflag==1 or nargout==0, plot the regions
%
%Author I. Halvorsen 990417
%

%
%Equation matrix:  VM*[R';V]=0
%[nc-1,nc+1]=size(VM)=[number of components-1,number of components+1]
%Each row in VW corresponds to an UW-root
%


%
%Code for Vmin(D) based on output from UWmulti
%IIMIN=find(Keys(:,1)==LK | Keys(:,2)==HK)
%[DS,IS]=sort(Ds(IIMIN))
%VminD=interp1(Ds(IIMIN(IS)),Vs(IIMIN(IS),Din);
%
% clf
if (isempty(theta))
  %theta=UWroots_fsolve(alfa,zf,qf); %%% OBS OBS ny UWroots funksjon!
  %theta=UWroots(alfa,zf,qf);
  theta=UWroots_special(alfa,zf,qf); %%% OBS OBS ny UWroots funksjon (intervallhalveringsmetode)!
end

nc=size(alfa,2);
nr=size(theta,2);
if (nr~=nc-1)
   error('Wrong sizes of alfa and theta')
end
IV=nc+1;

VM=zeros(nr,nc+1);
for i=1:nr
   VM(i,1:nc)=alfa.*zf./(alfa-theta(i));	
end
VM(:,IV)= -ones(nr,1);

if nargin<5
   F=1;
end	
if nargin<6
   inkey=[];
end	
if nargin<7
   plotflag=0;
end	
   
if ~isempty(inkey)
   Keys=inkey;
else
   %Make all legal key combinations
   Keys=[];
   for K1=1:nc-1,
   	for K2=K1+1:nc,
         Keys=[Keys;[K1 K2]];
      end
   end   
   Keys=[Keys;[0 1];[nc nc+1]];
end

Ds=[];
Vs=[];
Rs=[];
for i=1:size(Keys,1)
   K1=Keys(i,1);
   K2=Keys(i,2);      
   if(K1<1 | K2>nc)      
      if(K1==0)
         Ds=[Ds;0];
         Vs=[Vs;0];
      elseif(K2>nc)
         Ds=[Ds;1];
         Vs=[Vs;1-qf];
      end   
   else	
      NRV=K2-K1-1;			%Number of variable recoveries
      IRV=(K1+1):(K2-1);	%Variable recoveries
      AR=K1:(K2-1);			%Active roots
      RR=[ones(1,K1),zeros(1,nc-K1)];
      VM0=VM*[RR';0];
      sol=-inv(VM(AR,[IRV,IV]))*VM0(AR);
      RR(IRV)=sol(1:NRV)';
      D=zf*RR';
      V=sol(NRV+1);  
            
      Ds=[Ds;D];
      Vs=[Vs;V];
      Rs=[Rs;RR];
   end
end

Ds=Ds*F;
Vs=Vs*F;
%Vs=Vs-(1-qf)*F;
Rs=Rs*F;

if plotflag | nargout==0
   ISH=ishold;
   hold on
   if ~isempty(inkey)
      h=plot(Ds,Vs,'r');
   else      
      for i=1:nr,
         II=find(Keys(:,1)==0);
         LK=0;
         for j=1:nc-i
            II=[II;find(Keys(:,1)==j & Keys(:,2)==j+i)];
            II=[II;find(Keys(:,1)==j & Keys(:,2)==j+i+1)];
         end
         II=[II;find(Keys(:,2)==nc+1)];
         h(i)=plot(Ds(II),Vs(II),'r');
      end
   end
   
%    if (~ISH)
%       hold off
%    end   
end

return;




