%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   circuito_trifae.M                                                          	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc; echo off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definizione delle unità di misura
A=1; ohm=1; V=1; H=1; rad=1; sec=1;
% Definizione dei parametri del sistema
%
R=0.1*ohm;					% Resistenza di ogni singolo circuito
L=0.05*H;				   % Autoinduttanza di ogni singolo circuito
M=0.04*H;				   % Mutuainduttanza di ogni singolo circuito
w=50*rad/sec;				% Frequenza di rete
VM=100*V						% Ampiezza della tensione in ingresso
% 
% Si definisce "SYS" come il sistema lineare (LS) tempo continuo (TC) 
% caratterizzato dalle seguenti matrici MA, MB e MC:
%
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);
set(SYS,'InputName','V_M')
set(SYS,'StateName',['Iw_1';  'Iw_2'])
set(SYS,'OutputName',['Iw_1';  'Iw_2'])
SYS
%
% La risposta del sistema SYS ad un gradino di tensione in ingresso 
% si ottiene utilizzando il comando "step"
[Y t X]=step(SYS);
X=X*VM;					% La risposta del sistema è lineare rispetto 
Y=Y*VM;					% all'ampiezza VM della tensione di ingresso
Iw1=X(:,1);				% Corrente bifase routante: prima componente
Iw2=X(:,2);				% Corrente bifase routante: seconda componente
%
figure(1); clf
plot(t,[Iw1 Iw2])
title('Correnti Iw1 e Iw2')
ylabel('(A)'); grid
%
figure(2); clf; hold off
plot(Iw1,Iw2)
axis square; grid
title('Trajettoria nello spazio degli stati (Iw1, Iw2)')
ylabel('Iw2')
xlabel('Iw1')
%
% I poli del sistema SYS sono i seguenti:
pole(SYS)
%
% L'ampiezza delle correnti Iw1 e Iw2 a regime si ottiene imponendo 
% (d(Iw1)/dt)=0 e (d(Iw2)/dt)=0 :
%
figure(2); hold on
Iw_ss=-inv(MA)*MB*VM
plot(Iw_ss(1),Iw_ss(2),'rx')
%
% La matrice di trasformazione Tw ha la seguente struttura:
%      Tw=[cos(w*t)  -sin(w*t); 
%          sin(w*t)   cos(w*t)];
% L'andamento delle correnti sul piano bifase statico è il seguente:
Ib1=sum(([cos(w*t)  -sin(w*t)].*[Iw1';Iw2']')')';
Ib2=sum(([sin(w*t)   cos(w*t)].*[Iw1';Iw2']')')';
%
figure(3); clf
plot(t,[Ib1 Ib2])
title('Correnti Ib1 e Ib2')
ylabel('(A)'); grid
%
figure(4); clf; hold off
plot(Ib1,Ib2)
grid; axis square
title('Trajettoria nello spazio degli stati (Ib1, Ib2)')
ylabel('Ib2')
xlabel('Ib1')
%
% La matrice di trasformazione Tb ha la forma seguente:
% 
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,:)';
%
figure(5); clf
plot(t,I)
title('Corrent1 I1, I2 e I3')
ylabel('(A)'); grid
%
figure(6); clf; hold off
plot3(I1,I2,I2)
axis square; grid
%title('Trajettoria nello spazio degli stati (Ib1, Ib2)')
zlabel('I3')
ylabel('I2')
xlabel('I1')
return
