import cmath

Rs = 0.4
Ls = 0
Rf = 350
Lm = 0.15
Rr = 0.5
Lr = 1.5e-2
p = 2
Vs = 230
f = 50

w = 2*cmath.pi*f
g = 0.03


Zs = Rs+1j*Ls*w
Zm = Rf*1j*Lm*w/(Rf+1j*Lm*w)
Zr = 1j*Lr*w+Rr/g

Zeq = Zs+Zm*Zr/(Zm+Zr)

print(" ")
print(" ")
print("Zs =",Zs,"Ohm")
print("Zm =",Zm,"Ohm")
print("Zr =",Zr,"Ohm")
print("Zeq =",Zeq,"Ohm")
print(" ")

Is = Vs/Zeq
print("Courant stator : Is =",Is,"A  / efficace :",abs(Is),"A")
print(" ")

deltaV = Zs*Is
print("Chute de tension stator : deltaV =",deltaV,"V  / efficace :",abs(deltaV),"V")
print(" ")

E = Vs-deltaV
print("Fem : E =",E,"V  / efficace :",abs(E),"V")
print(" ")

Im = E/Zm
print("Courant magnetisant : Im =",Im,"A / efficace :",abs(Im),"A")
print(" ")

Ir = Is-Im
print("Courant rotor : Ir =",Ir,"A  / efficace :",abs(Ir),"A")
print(" ")

S = 3*Vs*Is.conjugate()
P = S.real
Q = S.imag
cosphi1 = P/abs(S)
cosphi2 = cmath.cos(cmath.phase(Is)).real
print("Puissance apparente : S =",S,"VA")
print("Puissance active : P =",P,"W")
print("Puissance reactive : Q =",Q,"VAR")
print("Facteur de puissance : cosphi =",cosphi1,"(P/module(S)), ou :",cosphi2,"cos(phase de Is)")
print(" ")

Pjs = 3*Rs*abs(Is)**2
Pf = 3*abs(E)**2/Rf
Pjr = 3*Rr*abs(Ir)**2
Pertes = Pjs+Pf+Pjr
Pu = P-Pertes
print("Pertes Joule statoriques : Pjs =",Pjs,"W")
print("Pertes fer : Pf =",Pf,"W")
print("Pertes Joule rotoriques : Pjr =",Pjr,"W")
print("Puissance utile : Pu =",Pu,"W")
print(" ")

W = (1-g)*w/p
N = W*30/cmath.pi
C = Pu/W
print("Vitesse : Omega =",W,"rad/s")
print("Vitesse : N =",N,"rpm")
print("Couple : C =",C,"N.m")
print(" ")

eta1 = C*W/P
eta2 = Pu/P
eta3 = (P-Pertes)/P
eta4 = Pu/(Pu+Pertes)
print("Rendement : eta =",eta1,"%")
print("    ou    : eta =",eta2,"%")
print("    ou    : eta =",eta3,"%")
print("    ou    : eta =",eta4,"%")
print(" ")


#####################

# Partie Bonus

import matplotlib.pyplot as plt

list_C = [0]
list_g = [0]
list_N = [60*f/p]
npoints = 1000
for k in range(1,npoints):
    g = k*1.0/npoints
    Zr = 1j*Lr*w+Rr/g
    Zeq = Zs+Zm*Zr/(Zm+Zr)
    Is = Vs/Zeq
    deltaV = Zs*Is
    E = Vs-deltaV
    Im = E/Zm
    Ir = Is-Im
    S = 3*Vs*Is.conjugate()
    P = S.real
    Pjs = 3*Rs*abs(Is)**2
    Pf = 3*abs(E)**2/Rf
    Pjr = 3*Rr*abs(Ir)**2
    Pertes = Pjs+Pf+Pjr
    Pu = P-Pertes
    W = (1-g)*w/p
    N = W*30/cmath.pi
    C = Pu/W
    list_C += [C]
    list_g += [g*1e2]
    list_N += [N]

plt.figure(1)
plt.plot(list_g,list_C)
plt.xlabel('Glissement (en %)')
plt.ylabel('Couple (en N.m)')
plt.title('Courbe $C(g)$')
plt.grid(True)
plt.show()

plt.figure(2)
plt.plot(list_N,list_C)
plt.xlabel('Vitese (en rpm)')
plt.ylabel('Couple (en N.m)')
plt.title('Courbe $C(N)$')
plt.grid(True)
plt.show()