In [1]:

    

import numpy as np
import sympy as sy
from sympy.utilities.codegen import codegen
import control.matlab as cm
import re


    

In [43]:

    

z = sy.symbols('z', real=False)
r1,s0,s1 = sy.symbols('r1,s0,s1', real=True)
h = 0.04
tau = 0.262
pc1 = -5-1j*5
pc2 = np.conjugate(pc1)
po1 = -4
hpt = h/tau
Km = 0.746*41.8
A2p = sy.poly(sy.simplify(sy.expand((z-np.exp(h*pc1))*(z-np.exp(h*pc2)))), z)
Acp = A2p*sy.poly((z-np.exp(h*po1)), 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()
print dioph


    

    




[1.0*r1 + 0.0905485642505141*s0 + 0.598553639437276, -1.85841144421398*r1 + 0.0860565404440641*s0 + 0.0905485642505141*s1 - 1.17944710048824, 0.858411444213975*r1 + 0.0860565404440641*s1 + 0.571209063848815]

In [44]:

    

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.632898354797954
0.379296078794350
-0.324459627611956
0.370876886395060

In [45]:

    

# Reorganize solution expression for matlab code generation
sol_expr = ('RST_DC_lab', [Bp.all_coeffs()[0], Bp.all_coeffs()[1],
                           Ap.all_coeffs()[1], Ap.all_coeffs()[2],
                           sol[r1], sol[s0], sol[s1], A2p.subs(z, 1)/Bp.subs(z,1), h,np.exp(h*po1) ])


    

In [46]:

    

# Export to matlab code
[(m_name, m_code)] = codegen(sol_expr, 'octave')


    

In [47]:

    

m_code = m_code.replace("out1", "b0").replace("out2", "b1").replace("out3", "a1").replace("out4", "a2")
m_code = m_code.replace("out5", "r1").replace("out6", "s0").replace("out7", "s1").replace("out8", "t0")
m_code = m_code.replace("out9", "h").replace("out10", "obsPole")
m_code = m_code.replace("function ", "% function ")
m_code = m_code.replace("end", "")
print m_code
with open("/home/kjartan/Dropbox/undervisning/tec/MR2007/labs/dc_rst_design.m", "w") as text_file:
    text_file.write(m_code)


    

    




% function [b0, b1, a1, a2, r1, s0, s1, t0, h, b00] = RST_DC_lab()
  %RST_DC_LAB  Autogenerated by sympy
  %   Code generated with sympy 0.7.6.1
  %
  %   See http://www.sympy.org/ for more information.
  %
  %   This file is part of 'project'

  b0 = 0.0905485642505141;
  b1 = 0.0860565404440641;
  a1 = -1.85841144421398;
  a2 = 0.858411444213975;
  r1 = -0.632898354797954;
  s0 = 0.37929607879435;
  s1 = -0.324459627611956;
  t0 = 0.37087688639506;
  h = 0.04;
  b00 = 0.852143788966211;

In [48]:

    

m_code.replace??


    

In [10]:

    

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.09055 z + 0.08606
----------------------
z^2 - 1.858 z + 0.8584

dt = 0.04


  0.09055 z + 0.08606
----------------------
z^2 - 1.858 z + 0.8584

dt = 0.04

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 [42]:

    

0.746*41.8


    

    
Out[42]:





31.182799999999997

In [ ]: