from math import *
from matplotlib import pyplot as plt
# L'équadiff est x''+ w0/Q x'+ w0² x = 0
# On peut l'écrire aussi x''+ 2 lam x'+ w0² x = 0
w0=1

# La précision demandée : si 0.01, on cherche la date à partir de laquelle 
# |x| reste toujours inférieur à 1% du max de x.
fractMax=0.0005

# à entrer en dur, le reste est calculé à partir de ça
# listQ=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.75,1.,1.5,2.,2.5]
listQ=[0.1+i/200 for i in range(600)]
# nombre de facteurs de qualités à traiter
N=len(listQ)

# création de listes annexes : lambdas et discriminants
listLam=[w0/(2*q) for q in listQ]
listDeltaSur4=[lam**2-w0**2 for lam in listLam]

# Fabrique la fonction de t solution de l'équadiff et vérifiant les CI passées en paramètre
def creeFonctionSolution(i, x0,xPt0):
    deltaSur4=listDeltaSur4[i]
    lam=listLam[i]
    if deltaSur4 ==0: # c'est alors (At+B)exp(-lam t)
        B=x0
        A=xPt0+B*lam
        return lambda t:(A*t+B)*exp(-lam*t)
    elif deltaSur4 >0: # apériodique, somme d'exp décroissantes
        r1=-lam+sqrt(deltaSur4)
        r2=-lam-sqrt(deltaSur4)
        A=(r2*x0-xPt0)/(r2-r1)
        B=(xPt0-r1*x0)/(r2-r1)
        return lambda t:A*exp(r1*t)+B*exp(r2*t)
    # sinon, pseudopériodique exp(-lam t)[A cos w t + B sin w t]
    w=sqrt(-deltaSur4)
    A=x0
    B=(xPt0+lam*x0)/w
    return lambda t:exp(-lam*t)*(A*cos(w*t)+B*sin(w*t))
    
# On les fixe ici : une interface utilisateur pourrait les demander
x0=1
xPt0=0

listFonct=[creeFonctionSolution(i,x0,xPt0) for i in range(N)]

def dateFin(i):
    fonction=listFonct[i]
    lam=listLam[i]
    # on est sûr qu'au delà, le signal est négligeable (théorie assez compliquée)
    majorant=max(1/lam,2*lam/w0**2)*10*abs(log(fractMax))
    # 1ère passe pour rechercher le max
    maximum=0
    t=0
    while(t<majorant):
        xt=abs(fonction(t))
        if xt>maximum:
            maximum=xt
        t+=majorant/10000
    # 2nde passe : la date d'arrêt est obtenue en trouvant la dernière valeur de
    # t qui dépasse en valeur absolue le seuil recherché
    # on commence donc par la fin (t valant majorant)
    while(t>0):
        xt=abs(fonction(t))
        if xt>fractMax*maximum:
            return t
        t-=majorant/10000
        
listDates=[dateFin(i) for i in range(N)]

def qDureeMin():
    indice=0
    dateMin=listDates[0]
    for i in range(1,N):
        if listDates[i]<dateMin :
            dateMin=listDates[i]
            indice=i
    return listQ[indice]
    
Qmin=qDureeMin()

plt.plot(listQ,listDates,'b-',label='Durée du transitoire ordre 2 (w0.t en fonction de Q)\n pour  un seuil de '+str(fractMax*100) +'%\n\n Qmin = '+str(Qmin))
plt.legend()
plt.show()