Subversion Repositories pentevo

Rev

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

  1. from sympy import *
  2.  
  3. A,B,C = symbols('A B C')
  4.  
  5. beta,gamma,kappa = symbols('beta gamma kappa')
  6.  
  7. # manual!
  8. ex6 = gamma**6 - B*gamma**4 + A*C*gamma**2 - C*C
  9.  
  10. s = solve(ex6,gamma,dict=True)
  11.  
  12. # s[1] is OK
  13.  
  14. gamma_sol = simplify(s[1][gamma])
  15. print('\ngamma={}'.format(gamma_sol))
  16. beta_sol = simplify( C*gamma_sol**(-2) )
  17. print('\nbeta={}'.format(beta_sol))
  18. kappa_sol = simplify( (A-beta_sol)/gamma_sol )
  19. print('\nkappa={}'.format(kappa_sol))
  20.  
  21.  
  22. w1_sol = simplify(1/beta_sol)
  23. print('\nw1={}'.format(w1_sol))
  24. w2_sol = simplify(1/gamma_sol)
  25. print('\nw2={}'.format(w2_sol))
  26. q_sol = simplify(1/kappa_sol)
  27. print('\nq={}'.format(q_sol))
  28.  
  29. # check initial equations
  30.  
  31. print('\n1/w1 + 1/(q*w2): {}'.format(simplify(1/w1_sol + 1/(q_sol*w2_sol))))
  32. print('\n1/w2**2 + 1/(q*w1*w2): {}'.format(simplify(1/w2_sol**2 + 1/(q_sol*w1_sol*w2_sol))))
  33. print('\n1/(w1*w2**2): {}'.format(simplify(1/(w1_sol*w2_sol**2))))
  34.  
  35.  
  36. prec=50
  37.  
  38. def ABC(f1_val,f2_val,q_val,prec=50):
  39.  
  40.     w1_val=Float(f1_val)*Float(2)*pi
  41.     w2_val=Float(f2_val)*Float(2)*pi
  42.     q_val =Float(q_val)
  43.  
  44.     q,w1,w2=symbols('q,w1,w2')
  45.     A = 1/w1+1/(q*w2)
  46.     B = 1/w2**2+1/(q*w1*w2)
  47.     C = 1/(w1*w2**2)
  48.  
  49.     A_val = A.evalf(subs={w1:w1_val,w2:w2_val,q:q_val},n=prec)
  50.     B_val = B.evalf(subs={w1:w1_val,w2:w2_val,q:q_val},n=prec)
  51.     C_val = C.evalf(subs={w1:w1_val,w2:w2_val,q:q_val},n=prec)
  52.  
  53.     return (A_val,B_val,C_val)
  54.  
  55. f1_val=30000
  56. f2_val=30000
  57. q_val=1
  58. (AA,BB,CC) = ABC(f1_val,f2_val,q_val,prec)
  59. print('\nfor f1={}, f2={}, q={}:\nA={}, B={}, C={}'.format(f1_val,f2_val,q_val,AA,BB,CC))
  60.  
  61. def F1F2Q(A_val,B_val,C_val,w1_sol,w2_sol,q_sol,prec=50):
  62.     w1_val = w1_sol.evalf(subs={A:A_val,B:B_val,C:C_val},n=prec)
  63.     w2_val = w2_sol.evalf(subs={A:A_val,B:B_val,C:C_val},n=prec)
  64.     q_val  =  q_sol.evalf(subs={A:A_val,B:B_val,C:C_val},n=prec)
  65.  
  66.     f1_val = N(w1_val/Float(2)/pi,n=prec)
  67.     f2_val = N(w2_val/Float(2)/pi,n=prec)
  68.  
  69.     return (f1_val,f2_val,q_val)
  70.  
  71. (ff1,ff2,qq) = F1F2Q(AA,BB,CC,w1_sol,w2_sol,q_sol,prec)
  72. print('\nfor previous A,B,C: f1={}, f2={}, q={}'.format(ff1,ff2,qq))
  73.  
  74.  
  75.  
  76. r1 = Float(1.2e3)
  77. r2 = Float(47e3)
  78. r3 = Float(68e3)
  79. r4 = Float(180e3)
  80. c1 = Float(0.0047e-6)
  81. c2 = Float(220e-12)
  82. c3 = Float(10e-12)
  83.  
  84. AA = N((c1*r1*r2 + c3*r1*r3 + c3*r1*r4 + c3*r2*r3 + c3*r2*r4 + c3*r3*r4)/(r1 + r2),prec)
  85. BB = N(c3*(c1*r1*r2*r3 + c1*r1*r2*r4 + c1*r1*r3*r4 + c2*r1*r3*r4 + c2*r2*r3*r4)/(r1 + r2),prec)
  86. CC = N(c1*c2*c3*r1*r2*r3*r4/(r1 + r2),prec)
  87.  
  88. print('\nfor r1={}, r2={}, r3={}, r4={}, c1={}, c2={}, c3={}:\nA={}, B={}, C={}'.format(r1,r2,r3,r4,c1,c2,c3,AA,BB,CC))
  89.  
  90. (ff1,ff2,qq) = F1F2Q(AA,BB,CC,w1_sol,w2_sol,q_sol,prec)
  91. print('\nfor previous A,B,C: f1={}, f2={}, q={}'.format(ff1,ff2,qq))
  92.  
  93. (AA,BB,CC) = ABC(ff1,ff2,qq,prec)
  94. print('\nfor f1,f2,q:\nA={}, B={}, C={}'.format(AA,BB,CC))
  95.  
  96.  
  97.  
  98.  
  99.  
  100. r1 = Float(1.2e3)
  101. r2 = Float(47e3)
  102. r3 = Float(68e3)
  103. r4 = Float(185038.95471783172862)
  104. c1 = Float(5.2453655748197971039E-9)
  105. c2 = Float(1.9333777853391275606E-10)
  106. c3 = Float(10e-12)
  107.  
  108. AA = N((c1*r1*r2 + c3*r1*r3 + c3*r1*r4 + c3*r2*r3 + c3*r2*r4 + c3*r3*r4)/(r1 + r2),prec)
  109. BB = N(c3*(c1*r1*r2*r3 + c1*r1*r2*r4 + c1*r1*r3*r4 + c2*r1*r3*r4 + c2*r2*r3*r4)/(r1 + r2),prec)
  110. CC = N(c1*c2*c3*r1*r2*r3*r4/(r1 + r2),prec)
  111.  
  112. print('\nfor r1={}, r2={}, r3={}, r4={}, c1={}, c2={}, c3={}:\nA={}, B={}, C={}'.format(r1,r2,r3,r4,c1,c2,c3,AA,BB,CC))
  113.  
  114. (ff1,ff2,qq) = F1F2Q(AA,BB,CC,w1_sol,w2_sol,q_sol,prec)
  115. print('\nfor previous A,B,C: f1={}, f2={}, q={}'.format(ff1,ff2,qq))
  116.  
  117. (AA,BB,CC) = ABC(ff1,ff2,qq,prec)
  118. print('\nfor f1,f2,q:\nA={}, B={}, C={}'.format(AA,BB,CC))
  119.  
  120.