In [1]:

    

import numpy as np
import sympy as sy
import control.matlab as cm


    

In [2]:

    

z = sy.symbols('z', real=False)
hh,a,r1,s0,s1 = sy.symbols('h,a,r1,s0,s1')
tau = 0.262
pc1 = -5-1j*5
pc2 = np.conjugate(pc1)
h = 0.04
hpt = h/tau
Km = 0.746*41.8
A2p = sy.poly((z-np.exp(h*pc1))*(z-np.exp(h*pc2)), z)
Acp = sy.poly((z-np.exp(h*pc1))*(z-np.exp(h*pc2))*(z - sy.exp(h*a)), z)
Ap = sy.poly((z-1)*(z-np.exp(-hpt)), z)
Bp = sy.poly(Km*tau*(hpt-1+np.exp(-hpt))*z + Km*tau*(1-np.exp(-hpt)-hpt*np.exp(-hpt)), z)
Rp = sy.poly(z+r1, z)
Sp = sy.poly(s0*z + s1, z)
dioph=(Ap*Rp+Bp*Sp-Acp).all_coeffs()


    

In [3]:

    

sol=sy.solve(dioph, (r1,s0,s1))
print sol[r1]
print sol[s0]
print sol[s1]

t0 = A2p.evalf(subs={z:1})/Bp.evalf(subs={z:1,})
print t0


    

    




-0.878368264722988*exp(0.04*a) + 0.115597706410771
-1.34327624390048*exp(0.04*a) + 1.52396058689998
0.972399357505391*exp(0.04*a) - 1.15308370050492
0.370876886395060

In [15]:

    

G = Km * cm.tf([1], [tau, 1, 0])
Gd = Km * cm.tf([tau*(hpt-1+np.exp(-hpt)), tau*(1-(1+hpt)*np.exp(-hpt))], [1, -(1+np.exp(-hpt)), np.exp(-hpt)], h)
Gd2 = cm.c2d(G, h)
print Gd
print Gd2


    

    




 0.001339 z + 0.001247
----------------------
z^2 - 1.807 z + 0.8071

dt = 0.03


 0.001339 z + 0.001247
----------------------
z^2 - 1.807 z + 0.8071

dt = 0.03

In [33]:

    

print A2p
print A2p.evalf(subs={z:1})
print Bp
print Bp.evalf(subs={z:1})


    

    




Poly(z**2 - 1.58555290309441*z - (-0.792776451547206 + 0.168974310731771*I)*(0.792776451547206 + 0.168974310731771*I), z, domain='EX')
0.0714939167206446
Poly(0.00133942860759726*z + 0.00124712240506047, z, domain='RR')
0.00258655101265772

In [4]:

    

0.3/(5*np.sqrt(2))


    

    
Out[4]:





0.042426406871192847

In [6]:

    

np.exp(-0.21)*np.sin(0.21)


    

    
Out[6]:





0.16897431073177133

In [7]:

    

np.exp(0.03*(-14))


    

    
Out[7]:





0.65704681981505675

In [ ]: