class Torseur2D(object):
    def __init__(self, resul, mom_ref, ref):
        '''torseur de résultante `resul` et de moment `mom_ref`
        exprimé au point de coordonnées `ref`.
        '''
        # Résultante
        self.resul = resul
        # Moment exprimé en ref
        self.mom_ref = mom_ref
        # Point de référence
        self.ref = ref
        
        print("Initialisation d'un Torseur")
    # end __init__()

# Création ("Instanciation") d'un torseur
T1 = Torseur2D((2,3), 1., (0,0))
T1.mom_ref
T1.ref
T1.resul

Torseur2D?

T1.resul = [4, 5]
T1.resul

# méthode à ajouter dans la définition de la classe Torseur2D
def calc_mom(self, pt):
    '''calcule le moment du torseur en un point `pt` quelconque'''
    # Formule du changement de point: M(X) = M(O) + R^OX
    # 1) Vecteur déplacement OX = X-O:
    delta_ref = (pt[0] - self.ref[0],
                 pt[1] - self.ref[1])
    # 2) Produit vectoriel R^OX
    delta_mom = self.resul[0]*delta_ref[1] - self.resul[1]*delta_ref[0]
    # 3) Nouveau moment:
    mom = self.mom_ref + delta_mom
    return mom
# end calc_mom()

# méthode à ajouter dans la définition de la classe Torseur2D
def set_ref(self, new_ref):
    '''change le point de référence'''
    # 1) Calcul du moment au nouveau point de référence
    self.mom_ref = self.calc_mom(new_ref)
    # 2) Mise à jour du point de référence
    self.ref = new_ref
# end set_ref()

def __str__(self):
    res = 'Torseur 2D '
    res += 'de Résultante {0} et de Moment {1}'.format(str(self.resul),
                                                        str(self.mom_ref))
    res += ' exprimé au point {0}'.format(str(self.ref))
    
    return res.encode('utf-8')
# end __str__

### Opération arithmétiques
def __add__(self, other):
    '''addition de 2 torseurs: T1 + T2
    Le pt de référence utilisé est celui de T1 (opérande gauche)
    '''
    # 1) addition des Résultantes:
    resul = (self.resul[0] + other.resul[0],
             self.resul[1] + other.resul[1])
    # 2) addition des moments au point de référence de self:
    mom_ref = self.mom_ref
    other_mom = other.calc_mom(self.ref)
    mom_ref = mom_ref + other_mom
    # 3) Création de l'objet Torseur:
    T = Torseur2D(resul, mom_ref, self.ref)
    return T
# end __add__

def __neg__(self):
    '''opposé du torseur'''
    resul = (-self.resul[0], -self.resul[1])
    mom_ref = -self.mom_ref
    T = Torseur2D(resul, mom_ref, self.ref)
    return T
# end __neg__

def __sub__(self, other):
    '''soustraction de 2 torseurs: T1-T2'''
    return self + (-other)

class Torseur2D(object):
    def __init__(self, resul, mom_ref, ref):
        '''torseur de résultante `resul` et de moment `mom_ref`
        exprimé au point de coordonnées `ref`.
        '''
        # Résultante
        self.resul = resul
        # Moment exprimé en ref
        self.mom_ref = mom_ref
        # Point de référence
        self.ref = ref        
        print("[Initialisation d'un nouveau Torseur]")
        #print(self)
    # end __init__()
    
    def calc_mom(self, pt):
        '''calcule le moment du torseur en un point `pt` quelconque'''
        # Formule du changement de point: M(X) = M(O) + R^OX
        # 1) Vecteur déplacement OX = X-O:
        delta_ref = (pt[0] - self.ref[0],
                     pt[1] - self.ref[1])
        # 2) Produit vectoriel R^OX
        delta_mom = self.resul[0]*delta_ref[1] - self.resul[1]*delta_ref[0]
        # 3) Nouveau moment:
        mom = self.mom_ref + delta_mom
        return mom
    # end calc_mom()

    def set_ref(self, new_ref):
        '''change le point de référence'''
        # 1) Calcul du moment au nouveau point de référence
        self.mom_ref = self.calc_mom(new_ref)
        # 2) Mise à jour du point de référence
        self.ref = new_ref
    # end set_ref()

    def __str__(self):
        res = u'Torseur 2D '
        res += u'de Résultante {0} et de Moment {1}'.format(str(self.resul),
                                                            str(self.mom_ref))
        res += u' exprimé au point {0}'.format(str(self.ref))
        
        return res.encode('utf-8')
    # end __str__
    
    ### Opération arithmétiques
    def __add__(self, other):
        '''addition de 2 torseurs: T1 + T2
        Le pt de référence utilisé est celui de T1 (opérande gauche)
        '''
        # 1) addition des Résultantes:
        resul = (self.resul[0] + other.resul[0],
                 self.resul[1] + other.resul[1])
        # 2) addition des moments au point de référence de self:
        mom_ref = self.mom_ref
        other_mom = other.calc_mom(self.ref)
        mom_ref = mom_ref + other_mom
        # 3) Création de l'objet Torseur:
        T = Torseur2D(resul, mom_ref, self.ref)
        return T
    # end __add__
    
    def __neg__(self):
        '''opposé du torseur'''
        resul = (-self.resul[0], -self.resul[1])
        mom_ref = -self.mom_ref
        T = Torseur2D(resul, mom_ref, self.ref)
        return T
    # end __neg__
    
    def __sub__(self, other):
        '''soustraction de 2 torseurs: T1-T2'''
        return self + (-other)

T1 = Torseur2D((1,1),  0, (0,0))
T2 = Torseur2D((1,1),  1, (0,0))
T3 = Torseur2D((1,1), -1, (0,0))
print(T3)

# Changer le point de référence :
T3.set_ref((-1,0))
print(T3)

print(-T1)

print(T1+T2) # torseurs déjà exprimés au même pt

print(T1+T3) # torseurs exprimés en des pts différents

# Importer le module numpy avec un nom raccourci :
import numpy as np

# Créer un "array" à partir d'une liste
l = np.array([10, 20, 30, 40]) # l n'est plus une liste mais un tableau NumPy

print(type(l))
l

l[0]

l+l

l*10

l**2

np.exp(l)
np.sqrt(l)

# On a déjà vu la fonction "constructeur" qui transforme tout type de listes (même imbriquées) en tableau
l = [1, 6, 10]
np.array(l)

# L'équivalent de `range` :
np.arange(10)

print(np.zeros(5))
print(np.ones(5))

# avec des listes imbriquées
np.array([[1,2,3], [4,5,6]])

# Matrice identité 3x3 :
np.eye(3)

a = np.zeros( (4,3) )
# attention aux parenthèses pour que (3,3) soit bien *un seul argument*
a

print(l.ndim)
a.ndim

a.shape

a = np.array([[1,2,3], [4,5,6], [7,8,9],[10,11,12]])
a

# sélectionner 1 élément
a[0,0]

# sélectionner la 1ère colonne 
a[:, 0]

# Méthode 1 : indexer avec la syntaxe start:stop:step
a[:, 0:2]

# Méthode 2 : indexer avec une liste (plus général)
a[:, [0, 1] ]

# sélectionner la 2ème (-> index 1) ligne
a[1, :]

# Étape 1 : créé le tableau booléen qui va sélectionner les éléments
a5 = a > 5
a5

# Étape 2 : utiliser le tableau booléen pour indexer :
a[a5]

# Écriture compacte (1 et 2 en une fois)
a[a>5]

print(np.array([1, 2, 3]).dtype) # tableau d'entier
print(np.array([1., 2., 3.]).dtype) # tableau de flottants

a = np.array([1, 2, 3])
# le tableau créé prend par défaut le type des donnnées d'entrée : ici `int`
a/2

a = np.array([1, 2, 3], dtype=float)
a/2

# Fabrication d'une matrice diagonale à partir d'une liste des éléments diagonaux
A = np.diag([1,2])
B = np.ones((2,2))*3
print(A)
print(B)
A * B

np.dot(A,B)

np.inner?

# Fabrication d'une matrice diagonale à partir d'une liste des éléments diagonaux
a = np.diag([1,2,3])
a

np.linalg.det(a) # det(a) vaut 1*2*3

# Consultons l'aide:
np.loadtxt?

# Lecture du fichier:
conso = np.loadtxt('./S2_Objets_et_NumPy/RTE_conso_2010-08_2011-07.csv',
                   delimiter=',', # les champs sont séparés par des virgules
                   skiprows=4, # ne pas lire les 4 premières lignes (du blabla)
                   usecols=[1]) # ne lire que la colonne *2* (comptage à partir de 0 !)

N = conso.shape[0]
print(u'{} pts chargés, soit {:.1f} jours'.format(N, N/24./4))

# Conversion MW -> GW
conso = conso/1000.

print(conso.min())
print(conso.max())
print(conso.mean())
print(conso.std())

# Méthode 1 : boucle for (très déconseillée pour les grands tableaux)

N80 = 0
for c in conso:
    if c >= 80:
        # on incrémente le compteur
        N80 += 1
N80

conso_80 = conso >= 80 # (GW)
conso_80

# Méthode 2 : indexaction logique pour ne sélectionner que les éléments vérifiant conso >= 80
len(conso[conso_80])

# Méthode 3 : jouer avec le transtypage implicite  (bool -> int)
# On calcule la "somme"
conso_80.sum()

np.savetxt?

# Importer l'API de Matplotlib:
import matplotlib.pyplot as plt
# Remarques :
# * cet import a déjà été réalisé en mode "Pylab"
# * il est par contre nécessaire dans un script

# Pour un travail interactif, on s'autorisera aussi
from matplotlib.pyplot import *

# Création des "données" : un sinus
x = np.linspace(0, 4*np.pi)
y = np.sin(x)

# Tracé
plt.plot(x,y) # ligne (bleue par défaut)
plt.plot(x,y, 'r+') # marqueurs "+" rouges

# Annotations:
plt.title('Un sinus')
plt.xlabel('temps t')
plt.grid(True);

# Vecteur temps en jours:
t = np.arange(N)/24./4 # attention au '24.'
t[:10]

# %pylab qt # afficher la figure dans une fenêtre interactive: nécessite IPython 1.0 (prévu pour juillet 2014)
plt.plot(t, conso)

plt.xlabel('temps (j)')
plt.ylabel(u'Puissance consommée (GW)');

# Zoom "programmé"
plt.plot(t, conso);
plt.title('Zoom sur 2 semaines consommation')
plt.xlabel('temps (j)')

plt.xlim(99, 113);
# il existe aussi plt.ylim

plt.hist(conso, bins=50, normed=False)
# (mettre normed=True pour tracer une densité de probabilité)
plt.xlabel(u'Puissance consommée (GW)');

# Chargement de données  (pseudo relevé expérimental d'un diagramme de Bode)
data = np.loadtxt('./S2_Objets_et_NumPy/bode_data.csv', delimiter=',')
data.shape # nombre de lignes, nombre de colonnes

# Séparation des données
f = data[:,0]
gain = data[:,1]
phase = data[:,2]

# ou en abrégé :
# f, gain, phase = data.T

# Affichage (tracé en échelle linéaire):
plt.plot(f, gain);

# Tracé en échelle semilog pour les abscisses (la fréquence)
semilogx(f, gain)
# Annotations:
title('Diagramme de Bode')
xlabel(u'fréquence (Hz)')
ylabel('gain (dB)');

fig = plt.figure('Bode', figsize=(6,6))
ax1 = fig.add_subplot(2,1,1, xscale='log',
                      title='Diagramme de Bode: gain (dB)')
ax1.plot(f, gain)
ax2 = fig.add_subplot(2,1,2, sharex=ax1,  # sharex -> synchronise les échelles. Pratique pour zoomer !
                      title=u'Phase (°)',
                      xlabel=u'fréquence (Hz)')
ax2.plot(f, phase, 'g')
# ajustement des "ticks"
ax2.set_yticks([-180, -135, -90, -45, 0])
# ou de façon plus automatique :
# ax2.set_yticks(range(0,-181,-45)) # (de 0 à -180°, tous les 45°)

fig.tight_layout() # "étire" les graphiques dans la fenêtre

# Observation du gain à basse fréquence:
print(f[:5])
print(gain[:5])
gain[:5].mean()

# dB -> linéaire
H0_dB = 6 # dB
print(10**(H0_dB/20.)) # Attention au point de "20."
H0 = 2

gain.max()

# Fréquence canonique: c'est environ la fréquence de résonance
i_max = gain.argmax()
f[i_max]

f0 = 1e3 

# Facteur de qualité : devinette
Q = 10

# On construit une fonction de transfert

# Dans cet exemple, on a déjà un vecteur de fréquence. À défaut, on pourrait utiliser :
# f = np.logspace(2, 4, 100)# de 100 à 10 kHz

# Tranfert du 2è ordre:
H_fit = H0/(1 + 2j*f/(f0*Q) - (f/f0)**2)
H_fit[:5] # (tableau de nombres complexes)

# Calcul du gain et de la phase
Habs = np.abs(H_fit) # `abs` donne le module du nombre complexe

gain_fit = 20*np.log10(Habs) # -> dB
gain_fit[:5]

phase_fit = np.angle(H_fit) * 180 / np.pi
phase_fit[:5]

# Par défaut, les graphiques se superposent:
semilogx(f, gain, '-+', label=u'données')
semilogx(f, gain_fit, label=u'modèle')
legend(loc='lower left');

# Fabrication de la chaîne de caractère de titre:
title_str = u'Modélisation du 2è ordre: H0={:.1f}, f0={:.1f} kHz, Q={:.1f}'.format(H0, f0/1000, Q)
# La syntaxe du "string formatting" : http://docs.python.org/2/library/string.html#formatstrings
title(title_str);

# Black :
plot(phase, gain)
plot(phase_fit, gain_fit)
plot(-180,0, 'r*', markersize=15) # point critique

# ajustement des "ticks"
xticks([-180, -135, -90, -45, 0])

xlabel(u'phase (°)')
ylabel('Gain (dB)');

# Diagramme de Nyquist

# Étape 1 : on reconstruit la réponse fréquentielle complexe
gain_lin = 10**(gain/20.) # dB -> lin
H = gain_lin*np.exp(1j*phase/180*np.pi)

# Étape 2 : tracé (utilisant les attributs .real et .image des tableaux de nombres complexes)
plot(H.real, H.imag)
plot(H_fit.real, H_fit.imag)
plot(-1,0, 'r*', markersize=10) # point critique

ax = gca()
ax.set_aspect('equal') # Même échelle en x et y
ylim(-11,1);