import math, numpy as np
from scipy.optimize import root
from delay_controller import critical_delay

def char(z,a,G,tau):
    lam=complex(z[0],z[1])
    q=lam+a+G*np.exp(-lam*tau)
    return [q.real,q.imag]

def roots(a,G,tau):
    # Multi-start Newton solves across the relevant low-frequency roots.
    out=[]
    for re in np.linspace(-3,2,9):
      for im in np.linspace(-12,12,49):
        sol=root(char,[re,im],args=(a,G,tau))
        if sol.success and np.linalg.norm(char(sol.x,a,G,tau))<1e-7:
          z=complex(*sol.x)
          if all(abs(z-w)>1e-5 for w in out): out.append(z)
    return out

def abscissa(a,G,tau):
    return max(z.real for z in roots(a,G,tau))

for G in [1.25,1.5,2,3]:
 t=critical_delay(1,G)
 print(G,t,abscissa(1,G,t),sorted(roots(1,G,t),key=lambda z:z.real)[-3:])