
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Script MATLAB ou GNU/Octave correspondant au
%    probleme du TD1 de Machines Electriques 
%
%     v2 : J. Fontchastagner / sept. 2021
%           derniere maj : sept. 2023
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% petit nettoyage
clear
clc 
close all

disp(' ')

% Donnees enonce
Un = 230;  % V
p = 2;
Sn = 6e3;  % VA
f = 50;    % Hz

Rs = 0.5;  % Ohms
Rr = 5;    % Ohms

Zch = 10*exp(j*pi/6); % Ohms
Vch = 127; % V

a0 = 0.5;    % N
a1 = 0.0032; % N/rad.s
Pfer = 300;  % W

Ie = [0:0.5:13 14 15];
E0 = [0 22.1 44.3 66.5 86.4 103.0 116.4 128.0 137.4 144.6 150.7 156.3 160.7...
      165 168 171 173 174 177 178 180 182 183 185 187 188 189 190 191];
      
Iemes = [0 1 2 3];
Iccmes = [0 6.1 11.9 18.2];

% test pour savoir si on est sous Matlab ou Octave 
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if (~isOctave)
  j=1i; % pour eviter warning dans l'editeur de Matlab...
end 

% question 1
% Vitesse nominale
Omegan = 2*pi*f/p;
Nn = Omegan*60/(2*pi);
disp(['Vitesse nominale : ',num2str(Omegan),' rad/s, soit : ',num2str(Nn),' rpm'])

% Courant nominal
% Sn = 3*Vn*In = sqrt(3)*Un*In
In = Sn/(sqrt(3)*Un);
disp(['Courant nominal : In = ',num2str(In),' A'])
disp(' ')

% Questions 2 a 4
kI = Iccmes(1:length(Iccmes))/Iemes(1:length(Iemes));
Icc = kI*Ie;

Xs = sqrt(E0.^2./Icc.^2-Rs^2);
Xs(1)=Xs(2);  % valeur en Ie = 0 A

% zone lineaire
klin = 1;
% zone lineaire : Xs ne varie pas plus d'1%
while (abs(Xs(klin+1)-Xs(klin))/Xs(klin) <= 1e-2)  
  klin = klin+1;  
end   
IeXslin = Ie(1:klin);
Xslin = sum(Xs(1:length(IeXslin)))/length(IeXslin);
disp(['Xslin = ', num2str(Xslin),' Ohms'])
kE = E0(1:length(IeXslin))/IeXslin;
disp(['Coeff. directeur E en lineaire : kE = ', num2str(kE)])
disp(['Coeff. directeur Icc : kI = ', num2str(kI)])
Xslin2 = sqrt((kE/kI)^2-Rs^2);
disp(['Xslin = ', num2str(Xslin2),' Ohms (par kI et kE)'])
disp(' ')

figure(1)
% grandeurs annexes pour trace 
%Xsprim = Xs;
Xslintrace = Xslin*ones(1,length(IeXslin)+2);
IeE0lin = Ie(1:3*klin);  % juste pour le trace
E0lin = kE*IeE0lin;
% trace
if (~isOctave) 
    % On est sous Matlab, on peut utiliser un axe a droite avec yyaxis
    yyaxis left
    plot(Ie,E0,'o-',Ie,Icc,'-','LineWidth',2,'MarkerSize',8);
    xlim([0 Ie(length(Ie))])
    ylabel('E0 (en V), Icc (en A)');
    hold on
    grid on
    plot(Iemes,Iccmes,'d','LineWidth',2)
    plot(IeE0lin,E0lin,'-.','LineWidth',2)
    xlabel('I_e (en A)')
    yyaxis right
    plot(Ie,Xs,'s-','LineWidth',2,'MarkerSize',8);
    ylim([0 (floor(Xs(1)/10)+1)*10])
    ylabel('Xs (en Ohms)')
    plot(Ie(1:length(Xslintrace)),Xslintrace,'-.','LineWidth',2)
else
    % On est sous Octave, on fait a l'ancienne
    Xsprim = 20*Xs;
    Xslintrace = 20*Xslintrace;
    p1 = plot(Ie,E0,'bo-',Ie,Icc,'r-',Ie,Xsprim,'ks-','LineWidth',2,'MarkerSize',8);
    xlim([0 Ie(length(Ie))])
    legend(p1,'E0 (en V)','Icc (en A)','20*Xs (en \Omega)','Location','northwest');
    xlabel('I_e (en A)')
    grid on
    hold on
    plot(Iemes,Iccmes,'rd','LineWidth',2)
    plot(Ie(1:length(Xslintrace)),Xslintrace,'k-.','LineWidth',2)
    plot(IeE0lin,E0lin,'b-.','LineWidth',2)
end
    
% question 5
Sch = 3*Vch^2/conj(Zch);
Pch = real(Sch);
Qch = imag(Sch);
disp(['Pch = ', num2str(Pch),' W'])
disp(['Qch = ', num2str(Qch),' VAR'])
disp(' ')

% question 6 et suivantes
E5A = interp1(Ie,E0,5);
Xs5A = interp1(Ie,Xs,5);
Sch5A = 3*Zch*(E5A/abs(Rs+j*Xs5A+Zch))^2;
Pch5A = real(Sch5A);
Qch5A = imag(Sch5A);
disp(['Pch(Ie=5A) = ', num2str(Pch5A),' W'])
disp(['Qch(Ie=5A) = ', num2str(Qch5A),' VAR'])
disp(' ')

% Calcul de Ie par modele lineaire
IeModeleLin = Vch/kE*abs(Rs+j*Xslin+Zch)/abs(Zch);
disp(['Ie (par modele lineaire) = ', num2str(IeModeleLin),' A'])
disp(' ')

% Calcul de Ie par BE

% methode 1
Vcalc = E0*abs(Zch)./abs(Rs+j*Xs+Zch);
Ieval = interp1(Vcalc,Ie,127);
disp(['Ie (par Behn Eshenbourg, methode 1) = ', num2str(Ieval),' A'])
Vchval = interp1(Ie,Vcalc,Ieval);
Vchlin = interp1(Ie,Vcalc,IeModeleLin);

figure(2)
plot(Ie,Vcalc,'bo-','LineWidth',2,'MarkerSize',8);
xlim([0 Ie(length(Ie))])
hold on
grid on
plot([Ieval Ieval 0],[0 Vchval Vchval],'--r','LineWidth',2)
plot([IeModeleLin IeModeleLin 0],[0 Vchlin Vchlin],'--g','LineWidth',2)
xlabel('I_e (en A)')
ylabel('V_s (en V)')
title(' Tension de sortie en fonction du courant d''excitation')

% methode 2
Ecalc = abs(Rs+j*Xs+Zch)/abs(Zch)*Vch;
Ieval2 = interp1(Ecalc-E0,Ie,0);
disp(['Ie (par Behn Eshenbourg, methode 2) = ', num2str(Ieval2),' A'])

figure(3)
plot(Ie,E0,'bs-',Ie,Ecalc,'r-','LineWidth',2,'MarkerSize',8)
xlim([0 Ie(length(Ie))])
grid on
hold on
Eval = interp1(Ie,Ecalc,Ieval2);
plot([Ieval2 Ieval2 0],[0 Eval Eval],'--r','LineWidth',2)
xlabel('I_e (en A)')
title(' FEM en fonction du courant d''excitation')

% methode 3 (iterative "a la main")
Ietest = IeModeleLin;
Xstest = interp1(Ie,Xs,Ietest);
Etest = interp1(Ie,E0,Ietest);
Ecalctest = abs(Rs+j*Xstest+Zch)/abs(Zch)*Vch;
iter = 1;
figure(4)
plot(Ie,E0,'bs-','LineWidth',2)
xlim([0 Ie(length(Ie))])
grid on
hold on
plot([Ietest Ietest],[0, Ecalctest],'r-','LineWidth',2)
while (abs(Ecalctest-Etest)/Etest > 1e-3)
  iter = iter+1;
  Ietestprec = Ietest;
  Ietest = interp1(E0,Ie,Ecalctest);
  plot([Ietestprec Ietest],[Ecalctest Ecalctest],'r-','LineWidth',2)
  Xstest = interp1(Ie,Xs,Ietest);
  Etest = interp1(Ie,E0,Ietest);
  Ecalctestprec = Ecalctest;
  Ecalctest = abs(Rs+j*Xstest+Zch)/abs(Zch)*Vch;
  plot([Ietest Ietest],[Ecalctestprec Ecalctest], 'r-','LineWidth',2) 
end
xlabel('I_e (en A)')
title('Procede iteratif')
ylabel('E (en V)')
Ieval3 = Ietest;
disp(['Ie (par Behn Eshenbourg, methode 3) = ', num2str(Ieval3),' A / (',...
     num2str(iter),' iterations)'])
disp(' ')

% question 9
Gammaf = a0+a1*Omegan;
Pmec = Gammaf*Omegan;
Pjs = 3*Rs*Vch^2/abs(Zch)^2;
Pjr = Rr*Ieval^2;
Pu = Pch;
eta = Pu/(Pu+Pjs+Pjr+Pmec+Pfer);
disp(['Pertes fer : Pfer = ', num2str(Pfer),' W'])
disp(['Pertes meca : Pmec = ', num2str(Pmec),' W'])
disp(['Pertes Joule stator : PJs = ', num2str(Pjs),' W'])
disp(['Pertes Joule rotor : PJr = ', num2str(Pjr),' W'])
disp(' ')
disp(['Puissance utile : Pu = ', num2str(Pch),' W'])
disp(' ')
disp(['Rendement : eta = ', num2str(eta*100),' %'])
disp(' ')





