#!/bin/env python
# -*- coding : utf-8 -*-

from brian2 import *

defaultclock.dt = 0.01*msecond

# 1- Je teste si j'ai le bon nombre de paramètres
if(len(sys.argv) != 4):
    print("Usage {} dt1(ms) dt2(ms) g/tau(mvolt)".format(sys.argv[0]))
    sys.exit(-1)

# 2- Je récupère les paramètres dt1, dt2, g
dt1 = float(sys.argv[1])*msecond
dt2 = float(sys.argv[2])*msecond
g = float(sys.argv[3])*mvolt

# 3- Je définis les temps d'émission des 3 spikes
t0 = 20 * msecond
t1 = t0 + dt1
t2 = t0 + dt2
tmax = t0 + 3 * t2 + 10*ms

# 4- Je définis mes constantes Vt,Vr,Vs,tau
...

# 5- Je me définis (cf annexe pour savoir comment définir une
#    fonction) une fonction qui retourne
#    le courant à injecter étant donnée une période de décharge
#    et une fonction  qui retourne le potentiel initial étant
#    un déphasage pour le train de spikes
def get_current(T):
    ....
def get_initial_potential(delta, I):
    ...

# 6- Je définis mes équations et DEUX groupes
#    Le groupe G est formé de 3 neurones
#    Le groupe Gd est formé d'un neurone (le détecteur)
....

# 7- Je connecte les neurones de G au neurone de G

# 8- Je calcule le courant à injecter pour que les neurones
#    déchargent avec une période de 40 ms.
I0 = get_current(40 * msecond)

# 9- J'initialise les potentiels des neurones 
G.V[0] = get_initial_potential(t0,I0)
G.V[1] = get_initial_potential(t1,I0)
G.V[2] = get_initial_potential(t2,I0)
Gd.V = Vs

# 10- J'initialise les courants des neurones
G.I = I0
Gd.I = 0

# 11- Je monitore les spikes de G et Gd
#     ainsi que le potentiel de Gd
spikemon = SpikeMonitor(G)
spikemond = SpikeMonitor(Gd)
statemon = StateMonitor(Gd,'V',record=0)

# 12- On lance la simulation suffisemment longtemps
run(tmax)

# 13- Et on affiche les résultats
figure()
subplot(211)
for i,ts in zip(spikemon.i, spikemon.t):
    plot([ts/msecond, ts/msecond],[i,i+1],'k')
for i,ts in zip(spikemond.i, spikemond.t):
    plot([ts/msecond,ts/msecond],[len(G),len(G)+1],'r')
ylim([-1, 4])
xlim([0, tmax/msecond])

subplot(212)
plot(statemon.t, statemon.V.T)
xlim([0, tmax/second])

show()