import scipy.integrate

u = 2.
tau = 10.

# Discrétisation:
N = 101
T = 50
t = np.linspace(0, T, N)
dt = t[1] - t[0]
print('dt = {:.2f}'.format(dt))
dt/tau

# Calcul analytique :
x_th = ...

# Tracé
plot(... );
title(u'solution théorique');

x_eul = np.zeros(N)
for i in range(...):
    x_eul[...] = ...

plot(t, x_eul, '-+', label='Euler explicite')
plot(t, x_th, label='sol. analytique')
title('Comparaison des solutions')
legend(loc='lower right');

odeint = scipy.integrate.odeint
odeint?

def deriv(y,t):
    '''dérivée calculé en y et au temps t'''
    return ...

# test de la fonction :
deriv(0,0)

x_ode = odeint(deriv, 0, t)

plot(t, x_eul, '--+', label='Euler explicite')
plot(t, x_th, label='sol. analytique')
plot(t, x_ode, '--x', label='odeint')
title('Comparaison des solutions')
legend(loc='lower right');

import scipy.optimize
from scipy.optimize import fmin

fmin?

# Chargement de données  (pseudo relevé expérimental d'un diagramme de Bode)
data = np.loadtxt(u'./S2_Objets_et_NumPy/bode_data.csv', delimiter=',')
# séparation des données :
f, gain, phase = data.T
# Tracé rapide:
semilogx(f, gain);

# Tranfert du 2è ordre:
def bode_pb2(f, H0, f0, Q):
    '''Bode d'un passe bas du 2ème ordre
    renvoie gain (dB), phase (°) 
    '''
    H = H0/(1 + 2j*f/(f0*Q) - (f/f0)**2)
    Habs = np.abs(H)
    gain = 20*np.log10(Habs) # -> dB
    phase = np.angle(H)*180/np.pi
    return (gain, phase) # ceci est un 'tuple' 

# tracé rapide:
semilogx(f, bode_pb2(f,1,1e3, 10)[0], 'g-');

def bode_err((H0, f0, Q)):
    '''erreur quadratique d'ajustement du modèle'''
    gain_fit, phase_fit = bode_pb2(f, H0, f0, Q)
    err_gain = ...
    err_phase = ....
    err = np.sum(... ) + np.sum(... )
    return err
bode_err((2, 1e3, 10))

x_opt = fmin(bode_err, x0=(1, 1e3, 10))
x_opt

# Dépaquetage des paramètres :
H0, f0, Q = x_opt
print('H0 = {:.2f}'.format(H0))
print('f0 = {:.0f} kHz'.format(f0))
print(' Q = {:.2f}'.format(Q))

# diagramme de Bode complet
fig = plt.figure('Bode', figsize=(6,6))

# gain
ax1 = fig.add_subplot(2,1,1, xscale='log',
                      title='Diagramme de Bode: gain (dB)')
ax1.plot(f, gain)
ax1.plot(f, bode_pb2(f, H0, f0, Q)[0], 'r-')
# phase
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, 'b')
ax2.plot(f, bode_pb2(f, H0, f0, Q)[1], 'r-')
# ajustement des "ticks"
ax2.set_yticks(np.arange(0,-181,-45)) # (de 0 à -180°, tous les 45°)
fig.tight_layout()