%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Circuito_trifase.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc; echo off
set(0,'DefaultFigureWindowStyle','docked') %%%% Figures 'docked' or 'normal'
MainString = mfilename;     % MainString
%%%%% Unità di misura
A=1; ohm=1; V=1; H=1; rad=1; sec=1; Stampa=1;
%%%%% Parametri del sistema
R=0.1*ohm;                  % Phase resistance
L=0.05*H;                   % Phase self-inductance coefficient
M=0.04*H;                   % Phase mutual-inductance coefficient
w=50*rad/sec;               % Electric frequency
VM=100*V;                   % Amplitude of the input voltages
%%%%% System matrices
MA=[ -R/(L-M)       w;
          -w  -R/(L-M)];
MB=[  sqrt(3/2)/(L-M); 0];
MC=[  1   0; 0   1];
SYS=ss(MA,MB,MC,0);         % System SYS in the rotating frame
set(SYS,'InputName','V_M')
set(SYS,'StateName',['Iw_1'; 'Iw_2'])
set(SYS,'OutputName',['Iw_1'; 'Iw_2'])
SYS
%%%%% Ste respose of system SYS
[Y t X]=step(SYS);
X=X*VM;                 % The system response is linear with respect ...
Y=Y*VM;                 % ... to the amplitude of the input voltage VM
Iw1=X(:,1);             % Rotating two-phase current: first term
Iw2=X(:,2);             % Rotating two-phase current: second term
Lw=1.2;                 % Line width
%
%%%%% Time behaviors of currents Iw1 Iw2
figure(1); clf
plot(t,[Iw1 Iw2],'LineWidth',Lw)
xlim([0 t(end)])
title('Currents Iw1 and Iw2')
xlabel('Time [s]')
ylabel('(A)'); grid
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
%
%%%%% Trajectory (Iw1 Iw2) in the rotating frame
figure(2); clf; hold off
plot(Iw1,Iw2,'LineWidth',Lw); hold on
Iw_ss=-inv(MA)*MB*VM;       % Steady state value of currents Iw1 and Iw2
plot(Iw_ss(1),Iw_ss(2),'rx','LineWidth',Lw)
axis square; grid
title('State space trajectories (Iw1, Iw2)')
ylabel('Iw2')
xlabel('Iw1')
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
%
%%%%%  The transformation matrix Tw has the following structure:
%      Tw=[cos(w*t)  -sin(w*t);
%          sin(w*t)   cos(w*t)];
% the behavior of the currents in the two-phase static plane is the following:
Ib1=sum(([cos(w*t)  -sin(w*t)].*[Iw1';Iw2']')')';
Ib2=sum(([sin(w*t)   cos(w*t)].*[Iw1';Iw2']')')';
%
%%%%% Time behaviors of currents Ib1 Ib2
figure(3); clf
plot(t,[Ib1 Ib2],'LineWidth',Lw)
xlim([0 t(end)])
title('Currents Ib1 and Ib2')
ylabel('(A)'); grid
xlabel('Time [s]')
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
%
%%%%% Trajectory (Ib1 Ib2) in the static two-phase plane
figure(4); clf; hold off
plot(Ib1,Ib2,'LineWidth',Lw)
grid; axis square
title('State space trajectory (Ib1, Ib2)')
ylabel('Ib2')
xlabel('Ib1')
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
%
%%%%% The transformaton matrix Tb has the following form:
Tb=[
   2/sqrt(6)          0;
  -1/sqrt(6)  1/sqrt(2);
  -1/sqrt(6) -1/sqrt(2)];
I=Tb*[Ib1';Ib2'];
I1=I(1,:)'; I2=I(2,:)'; I3=I(3,:)';
%
%%%%% Time behavior of currents I1, I2 and I3
figure(5); clf
plot(t,I,'LineWidth',Lw)
xlim([0 t(end)])
title('Currents I1, I2 and I3')
ylabel('(A)'); grid
xlabel('Time [s]')
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
%
%%%%% The trajectory (I1 I2 I3) in the static three-phase plane
figure(6); clf; hold off
plot3(I1,I2,I2,'LineWidth',Lw)
axis square; grid
zlabel('I3')
ylabel('I2')
xlabel('I1')
if Stampa;  eval(['print -depsc ' MainString '_'  num2str(gcf)]); end
return