Subversion Repositories pentevo

Rev

Rev 1246 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed

  1. from sympy import *
  2.  
  3. r1,r2,r3,r4,c1,c2,c3 = symbols('r1 r2 r3 r4 c1 c2 c3', real=True, positive=True)
  4.  
  5. ui,uo,u1,u2,u3 = symbols('ui uo u1 u2 u3')
  6.  
  7. # currents through appropriate elements.
  8. ir1,ir2,ir3,ir4 = symbols('ir1 ir2 ir3 ir4')
  9. ic1,ic2,ic3 = symbols('ic1 ic2 ic3')
  10.  
  11. # angular frequency (omega)
  12. w = symbols('w', real=True, positive=True)
  13.  
  14.  
  15. # make set of equations to calculate H
  16. eqs=[]
  17. #
  18. # u1 node, currents and voltages
  19. eqs += [ Eq( ir1, (ui-u1)/r1 ) ]
  20. eqs += [ Eq( ir2, (u2-u1)/r2 ) ]
  21. eqs += [ Eq( u1, (ir1+ir2)*(1/(I*w*c1)) ) ]
  22. #
  23. # u2 node
  24. eqs += [ Eq( ir3, (uo-u2)/r3 ) ]
  25. eqs += [ Eq( ir4, (u3-u2)/r4 ) ]
  26. eqs += [ Eq( u2, (ir3+ir4-ir2)*(1/(I*w*c2)) ) ]
  27. #
  28. # u3 node
  29. eqs += [ Eq( ir4, (uo-u3)/(1/(I*w*c3)) ) ]
  30. eqs += [ Eq( u3, 0 ) ]
  31.  
  32. # solve the set
  33. u_solve = solve( eqs, [u1,u2,u3,ir1,ir2,ir3,ir4,uo], dict=True, domain=S.Complexes )
  34.  
  35. #print(u_solve)
  36.  
  37. if len(u_solve)!=1:
  38.     sys.stderr.write("Many or no solutions: {} !\n".format(u_solve))
  39.     exit(1)
  40.  
  41. u_expr = u_solve[0]
  42.  
  43. h_expr = u_expr[uo]/ui
  44.  
  45. print(h_expr)
  46.  
  47. (n,d)=fraction(h_expr)
  48. d2=(d/n)*(-r3/(r1+r2))
  49. co=simplify(d2).collect(w)
  50.  
  51. print(co)
  52.  
  53.  
  54. w1,w2,q=symbols('w1,w2,q')
  55.  
  56. #w1 - 1st order, w2,q - 2nd order
  57. hd=(1+I*w/w1)*(1+I*w/(w2*q)-w**2/w2**2)
  58.  
  59. hdo = collect(expand(hd),w)
  60.  
  61. print('w**3: {}'.format(hdo.coeff(w,3)))
  62. print('w**2: {}'.format(hdo.coeff(w,2)))
  63. print('w**1: {}'.format(hdo.coeff(w,1)))
  64. print('w**0: {}'.format(hdo.coeff(w,0)))
  65.