function PhasePortraitU(A)%function PhasePortrait(A,B,C,D,u,DX) - Ute
% this program is based on a program which belongs to the course 'Dynamical System Theory for Engineers' by Martin Hasler (EPFL)
% Adaptation by Dr. Ute Feldmann

% PhasePortrait(A,B,C,D,u,DX)
% Make the plot of the phase portrait of a 2D dynamical system
% A,B,C,D are the matrix describing the system
% DX around the equilibrium to be considered for the initial conditions
% DX(1,:) variation on x
% DX(2,:) variation on y
%
%
% Example 1 - System 1 with autobalancing at the input badly dimensioned
% Si=10; Xmax=5; Ks=1.5; u=1; Ki=Si*Ks/(2*Xmax*Ks-Si)
% A=[-(Ki+Ks) 0; Ki -Ks]; B=[Ki*Xmax Si-Ki*Xmax]'; C=[Ks Ks]; D=[0];
% DX=[-4 3; -4 3];
% PhasePortrait(A,B,C,D,u,DX)
% rectangle('position',[0 0 Xmax Xmax],'edgecolor','m','facecolor','none')
% text(0,0,'(0,0)','verticalalignment','top','horizontalalignment','right');
% text(Xmax,Xmax,'(X_{max},X_{max})','verticalalignment','bottom','horizontalalignment','left');
% ax=axis;
% text(ax(1),ax(3)-0.05*diff(ax(3:4)),sprintf('Input balance BAD: S_i=%4.2f [Mb/s]; X_{max}=%4.2f [Mb]; K_s=%4.2f [s^{-1}]; Ki=%4.2f [s^{-1}]',Si,Xmax,Ks,Ki),'verticalalignment','top','horizontalalignment','left');
%
%
% Example 2 - System 1 with autobalancing at the input well dimensioned
% Si=10; Xmax=8; Ks=1.5; u=1; Ki=Si*Ks/(2*Xmax*Ks-Si)
% A=[-(Ki+Ks) 0; Ki -Ks]; B=[Ki*Xmax Si-Ki*Xmax]'; C=[Ks Ks]; D=[0];
% DX=[-5 6.5; -5 6.5];
% PhasePortrait(A,B,C,D,u,DX)
% rectangle('position',[0 0 Xmax Xmax],'edgecolor','m','facecolor','none')
% text(0,0,'(0,0)','verticalalignment','top','horizontalalignment','right');
% text(Xmax,Xmax,'(X_{max},X_{max})','verticalalignment','bottom','horizontalalignment','left');
% ax=axis;
% text(ax(1),ax(3)-0.05*diff(ax(3:4)),sprintf('Input balance GOOD: S_i=%4.2f [Mb/s]; X_{max}=%4.2f [Mb]; K_s=%4.2f [s^{-1}]; Ki=%4.2f [s^{-1}]',Si,Xmax,Ks,Ki),'verticalalignment','top','horizontalalignment','left');
%
%
% Example 3 - System 2 with balancing demon
% Si=10; Ks=3; Kb=30; u=1; Xmax=4;
% A=[-(Ks+Kb) Kb; Kb -(Ks+Kb)]; B=[Si 0]'; C=[Ks Ks]; D=[0];
% DX=[-3 5; -3 5];
% PhasePortrait(A,B,C,D,u,DX)
% rectangle('position',[0 0 Xmax*1.15 Xmax],'edgecolor','b','facecolor','none')
% rectangle('position',[0 0 Xmax Xmax*1.15],'edgecolor','r','facecolor','none')
% rectangle('position',[0 0 Xmax Xmax],'edgecolor','m','facecolor','none')
% text(0,0,'(0,0)','verticalalignment','top','horizontalalignment','right');
% text(Xmax,Xmax,'(X_{max},X_{max})','verticalalignment','bottom','horizontalalignment','left');
% ax=axis;
% text(ax(1),ax(3)-0.05*diff(ax(3:4)),sprintf('Balancing demon: S_i=%4.2f [Mb/s]; K_s=%4.2f [s^{-1}]; Kb=%4.2f [s^{-1}]; X_{max}=%4.2f [Mb]',Si,Ks,Kb,Xmax),'verticalalignment','top','horizontalalignment','left');
%
%
% Example 4 - Random system
% A=randn(2,2); B=randn(2,1); C=randn(1,2); D=[0]; u=0;
% DX=[-5 5; -5 5];
% PhasePortrait(A,B,C,D,u,DX)
%
%
% Example 5 - Star node
% T=randn(2,2);
% A=-T*[1 0; 0 1]*inv(T); B=randn(2,1); C=randn(1,2); D=[0]; u=0;
% DX=[-5 5; -5 5];
% PhasePortrait(A,B,C,D,u,DX)
%
%
% Example 6 - neutral node
% T=randn(2,2);
% A=-T*[1 1; 0 1]*inv(T); B=randn(2,1); C=randn(1,2); D=[0]; u=0;
% DX=[-5 5; -5 5];
% PhasePortrait(A,B,C,D,u,DX)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameters
nStreams = 10; % Number of streams to plot
DX=[-5 5; -5 5];% state space to be plotted - Ute

% Compute the equilibrium
Xeq=[0;0];%-A\B*u -Ute

% The eigenvectors
[EV,Ev]=eig(A);
% sort them in leading and non leading
[sv,idx]=sort(diag(abs(real(Ev))));
EV=EV(:,idx);
display('Eigenwerte:'),a=diag(Ev), 
Ev=diag(a(idx));

% Initial conditions
XR=Xeq*[1 1]+DX;
xxi=linspace(XR(1,1),XR(1,2),nStreams)';
yyi=linspace(XR(2,1),XR(2,2),nStreams)';
yyi([1 end])=[];
Xi=[xxi XR(2,1)*ones(size(xxi)); XR(1,2)*ones(size(yyi)) yyi; xxi XR(2,2)*ones(size(xxi)); XR(1,1)*ones(size(yyi)) yyi; 0 0];


% The stream
[X,Y]=meshgrid(linspace(XR(1,1),XR(1,2),50),linspace(XR(2,1),XR(2,2),50));
XYDOT=A*[X(:) Y(:)]';%+repmat(B*u,[1 size(X(:),1)]); - Ute
XDOT=reshape(XYDOT(1,:),size(X));
YDOT=reshape(XYDOT(2,:),size(Y));
XYstreamfw=stream2(X,Y,XDOT,YDOT,Xi(:,1),Xi(:,2),[0.05 1500]); 
XYstreambw=stream2(X,Y,-XDOT,-YDOT,Xi(:,1),Xi(:,2),[0.05 1500]); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PLOTTING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
% Position of the frames
  hlow=1.5;
  hhig=1.5;
  vlow=1.5;
  vhig=1.5;

% 2D Box size
  wd=(17.5-hlow-hhig);
  ht=wd;
  
  % Printout the final size
  [wd+hlow+hhig ht+vlow+vhig];

% Fonts sizes
  fnt=14;
  fnl=18;
  fnlb=18;

% Open the Figure
figure(...
       'units','centimeters',...
       'position',[hlow vlow wd+hlow+hhig ht+vlow+vhig], ...
       'PaperUnits','centimeters', ...
       'PaperOrientation','portrait', ...
       'PaperPosition',[hlow vlow wd+hlow+hhig ht+vlow+vhig], ...
       'PaperType','A4');
%       'Color','none' ...
      
  
% prepare the axes
ha=axes(...
    'units','centimeters',...
    'position',[hlow vlow wd ht],...
    'fontsize',fnt,...
    'fontname','helvetica',...
    'fontangle','normal',...
    'dataaspectratio',[1 1 1], ...
    'box','on', ...
    'visible','on', ...
    'Color','none' ...
    );
axis(reshape(XR',[1 4]));   


% Plot the eigenvectors (at maximum are two)
  if isreal(Ev(1)), % otherwise the eigenvectors are complex and conjugate, reminder we sorted them thus the first one is leading
     % radius for plotting
     rr=sqrt(sum(diff(DX,1,2).^2));
     % Eigenvector extrema
     EVR=EV*rr*[-1 1 0 0; 0 0 -1 1];
     % Arrows points
     dd=sign(diag(Ev));
     AP=1e-1*EV*rr*[1-dd(1)*0.15:dd(1)*0.1:1+dd(1)*0.15 0 0 0 0; 0 0 0 0 1-dd(2)*0.15:dd(2)*0.1:1+dd(2)*0.15];
     
     % color choice, red unstable green stable
     cs=(1+sign(diag(Ev)))/2;
     clr=[cs 1-cs zeros(size(cs))]*0.75;
     for n=1:2,
         line(Xeq(1)+EVR(1,(n-1)*2+(1:2)),Xeq(2)+EVR(2,(n-1)*2+(1:2)),'linewidth',3,'color',clr(n,:));
         if n==2,
             arrowline(Xeq(1)+AP(1,(n-1)*4+(3:4)),Xeq(2)+AP(2,(n-1)*4+(3:4)),'color',clr(n,:),'arrowsize',500);
             arrowline(Xeq(1)-AP(1,(n-1)*4+(3:4)),Xeq(2)-AP(2,(n-1)*4+(3:4)),'color',clr(n,:),'arrowsize',500);
         end;
         arrowline(Xeq(1)+AP(1,(n-1)*4+(1:2)),Xeq(2)+AP(2,(n-1)*4+(1:2)),'color',clr(n,:),'arrowsize',500);
         arrowline(Xeq(1)-AP(1,(n-1)*4+(1:2)),Xeq(2)-AP(2,(n-1)*4+(1:2)),'color',clr(n,:),'arrowsize',500);
     end;
  end
  

% Plot the stream
  for n=1:prod(size(XYstreamfw)),
      line(XYstreamfw{n}(:,1),XYstreamfw{n}(:,2),'color','k');
      idf=round(size(XYstreamfw{n},1)/4.*[1 1.05]);
      if all(idf) & (abs(diff(idf)>1)),
         arrowline(XYstreamfw{n}(idf,1),XYstreamfw{n}(idf,2),'color','k','arrowsize',250);
      end;
  end

  for n=1:prod(size(XYstreambw)),
      line(XYstreambw{n}(:,1),XYstreambw{n}(:,2),'color','k');
      idf=round(size(XYstreambw{n},1)/4.*[1 0.95]);
      if all(idf) & (abs(diff(idf))>1),
         arrowline(XYstreambw{n}(idf,1),XYstreambw{n}(idf,2),'color','k','arrowsize',250);
      end;
  end
  

% Plot the equilibrium
  line(Xeq(1),Xeq(2),'linestyle','none','marker','.','markersize',28,'markerfacecolor','k','markeredgecolor','k');
  

% Axis labels
  ax=axis;
  lt=[ax(2)-diff(ax(1:2))*0.02 ax(3)+diff(ax(3:4))*0.02];
  text(lt(1),lt(2),'x_1', ...
      'fontname','times',...
      'fontangle','italic',...
      'fontsize',fnl, ...
      'fontweight', 'normal', ...
      'interpreter','tex',...
      'HorizontalAlignment','right',...
      'VerticalAlignment','bottom' ...
      );
  
  lt=[ax(1)+diff(ax(1:2))*0.02 ax(4)-diff(ax(3:4))*0.01];
  text(lt(1),lt(2),'x_2', ...
      'fontname','times',...
      'fontangle','italic',...
      'fontsize',fnl, ...
      'fontweight', 'normal', ...
      'interpreter','tex',...
      'HorizontalAlignment','left',...
      'VerticalAlignment','top' ...
      );


% END
  return;