%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Script MATLAB ou GNU/Octave correspondant au
%    probleme du TD4 de Machines Electriques
%
%     v2.1 : J. Fontcchastagner / dec. 2022
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all
clc
disp(' ')

% Donnees enonce
V1 = 230;
I10 = 80;

V1cc = 51.3;
I1cc = 118;
P1cc = 3e3;

f = 50;
w = 2*pi*f;
p = 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


% Partie 1
%-----------

% question 1
L1w = V1/I10;
disp([' L1w = ', num2str(L1w),' ohms'])
disp(' ')

% question 2
S1cc = 3*V1cc*I1cc;
Q1cc = sqrt(S1cc^2-P1cc^2);
cosphicc = P1cc/S1cc;
sinphicc = Q1cc/S1cc;

I1ccC = I1cc*(cosphicc-j*sinphicc);
I2C = I1ccC-V1cc/(j*L1w);

Q2cc = Q1cc-3*V1cc^2/L1w;
Rprim2 = P1cc/(3*abs(I2C)^2);
Nprim2w = Q2cc/(3*abs(I2C)^2);

disp([' Nprim2w = ',num2str(Nprim2w),' ohms'])
disp([' Rprim2 = ',num2str(Rprim2),' ohms'])
disp(' ')

gM = Rprim2/Nprim2w;
Gamma_max = 3*p*V1^2/(2*Nprim2w*w);
Gamma_d = 2*Gamma_max/(gM+1/gM);

disp([' Gamma_max = ',num2str(Gamma_max),' N.m'])
disp([' gM = ',num2str(gM*100),' %'])
disp([' Gamma_d = ',num2str(Gamma_d),' N.m'])
disp(' ')

nfig = 1;
figure(nfig)
npoints = 200;
g = linspace(0,1,npoints);
Gamma = 2*Gamma_max./(gM./g+g/gM);
plot(g,Gamma,"linewidth",3)
xlabel('g')
ylabel('Couple (en N.m)')
grid on
hold on

% question 3
a0 = 40;
a1 = 1e-1;
a2 = 1e-2;
Omegas = w/p;
Gamma_cht = a0+a1*Omegas.*(1-g)+a2.*(Omegas.*(1-g)).^2;
plot(g,Gamma_cht,"g--","linewidth",3)
legend('Couple MAS','Couple de charge')
title(' Couples en fonction du glissement')

nfig = nfig+1;
figure(nfig)
plot(Omegas.*(1-g)*30/pi,Gamma, ...
    Omegas.*(1-g)*30/pi,Gamma_cht,"g--","linewidth",3);
grid
legend('Couple MAS','Couple de charge','location','Northwest')
xlabel('N (en tours/min)')
ylabel('Couple (en N.m)')
title(' Couples en fonction de la vitesse')

% Partie 2
%-------------

% question 1
% Point de départ obtenu par linéarisation de la formule du couple
P = [a2*Omegas^2, (a1*Omegas+2*Gamma_max/gM), (a0-2*Gamma_max/gM)];
X = roots(P);
if ( abs(X(1)) < 1 )
  X0 = X(1);
else
  X0 = X(2);
end
g0 = 1-X0;

% Résolution de l'équation non-linéaire
donne_g = @(g) a0+a1*Omegas*(1-g)+a2*(Omegas*(1-g))^2 - 2*Gamma_max/(gM/g+g/gM);
if (~isOctave)
  % On est sous Matlab
  options = optimoptions('fsolve');
  options.Display = 'off';
  g = fsolve(donne_g,g0,options);
else
  % On est sous Octave
  pkg load optim;
  g = fsolve(donne_g,g0);
end

N = (1-g)*Omegas*30/pi;
disp([' Vitesse de l''ensemble : N = ',num2str(N), ' tours/min'])
disp(' ')

Z1 = j*L1w;
Z2 = Rprim2/g+j*Nprim2w;
Zeq = Z1*Z2/(Z1+Z2);
I1C = V1/Zeq;
I1 = abs(I1C);
S1 = 3*V1*conj(I1C);
P1 = real(S1);
Gamma = 2*Gamma_max/(gM/g+g/gM);
eta = Gamma*(1-g)*Omegas/P1;
un_moins_g = 1-g;

disp([' Courant absorbé : I1 = ',num2str(I1),' A'])
disp([' Facteur de puissance : cosphi = ',num2str(P1/abs(S1))])
disp([' Rendement : eta = ',num2str(eta*100),' %'])
disp(' ')


% Partie 3
%-----------

Omega2 = 1000*pi/30;
g0 = (Omegas-Omega2)/Omegas;
Gamma_ch = a0+a1*Omega2+a2*(Omega2)^2;

% question a
kV1 = sqrt(Gamma_ch/(2*Gamma_max/(gM/g0+g0/gM)));
disp([' On doit diminuer la tension de ',num2str((1-kV1)*V1),' V'])
disp('   pour fonctionner à 1000 tours.')
disp(' ')

nfig = nfig+1;
figure(nfig)
npoints = 200;
g = linspace(0,1,npoints);
Gamma = 2*Gamma_max./(gM./g+g/gM);
plot(g,Gamma,"b-","linewidth",3)
grid on
hold on
Gamma_max2 = 3*p*(kV1*V1)^2/(2*Nprim2w*w);
Gamma2 = 2*Gamma_max2./(gM./g+g/gM);
plot(g,Gamma2,"r-","linewidth",3)
Gamma_cht = a0+a1*Omegas.*(1-g)+a2.*(Omegas.*(1-g)).^2;
plot(g,Gamma_cht,"g--","linewidth",3)

% question b
gMprim = (Gamma_max/Gamma_ch+sqrt((Gamma_max/Gamma_ch)^2-1))*g0;
Rsup = gMprim*Nprim2w-Rprim2;
disp([' On doit jouter une résistance de Rsup = ',num2str(Rsup),' ohms'])
disp(['   au rotor, pour fonctionner à 1000 tours/min'])
disp([' Et dans ce cas :'])


nfig = nfig+1;
figure(nfig)
g=0:1/(npoints-1):1;
Gamma = 2*Gamma_max./(gM./g+g/gM);
plot(g,Gamma,"b-","linewidth",3)
hold on
grid on
Gamma2 = 2*Gamma_max./(gMprim./g+g/gMprim);
plot(g,Gamma2,"r-","linewidth",3)
Gamma_cht = a0+a1*Omegas.*(1-g)+a2.*(Omegas.*(1-g)).^2;
plot(g,Gamma_cht,"g--","linewidth",3)

Z1 = j*L1w;
Z2 = (Rprim2+Rsup)/g0+j*Nprim2w;
Zeq = Z1*Z2/(Z1+Z2);
I1C = V1/Zeq;
I1 = abs(I1C);
S1 = 3*V1*conj(I1C);
P1 = real(S1);
Q1 = imag(S1);
cosphi = P1/abs(S1);
Gamma = 2*Gamma_max/(gMprim/g0+g0/gMprim);
eta = Gamma*(1-g0)*Omegas/P1;
disp([' Courant absorbé : I1 = ',num2str(I1),' A'])
disp([' Puissance active : P1 = ',num2str(P1),' W'])
disp([' Puissance active : Q1 = ',num2str(Q1),' VAR'])
disp([' Facteur de puissance : cosphi = ',num2str(cosphi)])
disp([' Rendement : eta = ',num2str(eta*100),' %'])
disp(' ')

% question c :
% en V/f constant
Omega0 = Omega2;
donne_w = @(w) Gamma_ch - 2*Gamma_max/((1-Omega0*p/w)/(gM*2*pi*f/w)+(gM*2*pi*f/w)/(1-Omega0*p/w));
w0 = p*Omega0;
if (~isOctave)
  % On est sous Matlab
  w2 = fsolve(donne_w,w0,options);
  else
  % On est sous Octave
  w2 = fsolve(donne_w,w0);
  end
f2 = w2/2/pi;
V12 = V1/f*f2;

disp([' Loi V/f constant :'])
disp([' Tension : V1 = ',num2str(V12),' V'])
disp([' Fréquence : f = ',num2str(w2/(2*pi)),' Hz'])

nfig = nfig+1;
figure(nfig)
Omega = linspace(0,Omegas,npoints);
g=1-Omega/Omegas;
Gamma = 2*Gamma_max./(gM./g+g/gM);
plot(Omega*30/pi,Gamma,"b-","linewidth",3)
hold on
grid on
Gamma_cht = a0+a1*Omega+a2.*(Omega).^2;
plot(Omega*30/pi,Gamma_cht,"g--","linewidth",3)

Omegas2 = w2/p;
Omega2 = linspace(0,Omegas2,npoints);
g=1-Omega2/Omegas2;
gMprim = Rprim2/(Nprim2w/w*w2);
Gamma2 = 2*Gamma_max./(gMprim./g+g/gMprim);
plot(Omega2*30/pi,Gamma2,"r-","linewidth",3)

g=1-Omega0/Omegas2;
Z1 = j*L1w/w*w2;
Z2 = (Rprim2)/g+j*Nprim2w/w*w2;
Zeq = Z1*Z2/(Z1+Z2);
I1C = V12/Zeq;
I1 = abs(I1C);
S1 = 3*V12*conj(I1C);
P1 = real(S1);
Q1 = imag(S1);
cosphi = P1/abs(S1);
Gamma = 2*Gamma_max/(gMprim/g+g/gMprim);
eta = Gamma*(1-g)*Omegas2/P1;
disp([' Courant absorbé : I1 = ',num2str(I1),' A'])
disp([' Puissance active : P1 = ',num2str(P1),' W'])
disp([' Puissance active : Q1 = ',num2str(Q1),' VAR'])
disp([' Facteur de puissance : cosphi = ',num2str(cosphi)])
disp([' Rendement : eta = ',num2str(eta*100),' %'])
disp(' ')


