1import math, numpy as np
 2from scipy.optimize import root
 3from delay_controller import critical_delay
 4
 5def char(z,a,G,tau):
 6    lam=complex(z[0],z[1])
 7    q=lam+a+G*np.exp(-lam*tau)
 8    return [q.real,q.imag]
 9
10def roots(a,G,tau):
11    # Multi-start Newton solves across the relevant low-frequency roots.
12    out=[]
13    for re in np.linspace(-3,2,9):
14      for im in np.linspace(-12,12,49):
15        sol=root(char,[re,im],args=(a,G,tau))
16        if sol.success and np.linalg.norm(char(sol.x,a,G,tau))<1e-7:
17          z=complex(*sol.x)
18          if all(abs(z-w)>1e-5 for w in out): out.append(z)
19    return out
20
21def abscissa(a,G,tau):
22    return max(z.real for z in roots(a,G,tau))
23
24for G in [1.25,1.5,2,3]:
25 t=critical_delay(1,G)
26 print(G,t,abscissa(1,G,t),sorted(roots(1,G,t),key=lambda z:z.real)[-3:])