Subversion Repositories pentevo

Rev

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

  1. #!/usr/bin/env python
  2.  
  3. import sys
  4. import math
  5. import cmath
  6. import random
  7.  
  8. #e24 = [1.0, 1.1, 1.2, 1.3, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.3, 4.7, 5.1, 5.6, 6.2, 6.8, 7.5, 8.2, 9.1];
  9. #e6  = [1.0,                1.5,                2.2,                3.3,                4.7,                6.8               ];
  10.  
  11.  
  12.  
  13. # these formulae calculated symbolically by the function above and are set here in comments
  14. # as reference for further numerical calculations.
  15. # see mfb3_*play.py for ways to calculate
  16. """
  17. r1,r2,r3,c3 -- select arbitrary/as needed,
  18. calculate c1,c2,r4
  19.  
  20. (see mfb3.png)
  21.  
  22. w1   - 1st order,
  23. w2,q - 2nd order
  24.  
  25. 1/w1 + 1/(q*w2)       = A = (c1*r1*r2 + c3*r1*r3 + c3*r1*r4 + c3*r2*r3 + c3*r2*r4 + c3*r3*r4)/(r1 + r2)
  26. 1/w2**2 + 1/(q*w1*w2) = B = c3*(c1*r1*r2*r3 + c1*r1*r2*r4 + c1*r1*r3*r4 + c2*r1*r3*r4 + c2*r2*r3*r4)/(r1 + r2)
  27. 1/(w1*w2**2)          = C = c1*c2*c3*r1*r2*r3*r4/(r1 + r2)
  28.  
  29. w1 = (2**(1/3)*A*C - 2**(1/3)*B**2/3 + B*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)/3 - (18*A*B*C - 4*B**3 - 54*C**2 + 6*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(2/3)/6)/(C*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))
  30. w2 = 2**(1/6)*sqrt(3)/sqrt((3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))
  31. q  = 2**(5/6)*sqrt(3)*sqrt((3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(6*(A*(3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)) - 3*2**(1/3)*C*(9*A*B*C - 2*B**3 - 27*C**2 + 3*sqrt(3)*sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)))
  32.  
  33. r4 = (-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3))
  34. c1 = -c3*(r1 + r2 + r3)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))/(r1*r2) + (A - c3*r3)*(r1 + r2)/(r1*r2)
  35. c2 = C*(r1 + r2)/(c3*r3*(-c3*(r1 + r2 + r3)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))**2 + (A - c3*r3)*(r1 + r2)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))))
  36.  
  37. 1. select r1,r2,r3,c3,w1,w2,q (-r3/(r1+r2) makes k)
  38. 2. calc A,B,C from w1,w2,q
  39. 3. calc r4,c1,c2 from A,B,C,r1,r2,r3,c3
  40. 4. fit r4,c1,c2 into Exx
  41. 5. check resulting w1,w2,q,k
  42. """
  43.  
  44.  
  45. def select_pvalues_list(values='e24'):
  46. # returns [1.0, ..., 9.1] list (for e24)
  47.  
  48.     if values=='e24':
  49.         r = [1.0, 1.1, 1.2, 1.3, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.3, 4.7, 5.1, 5.6, 6.2, 6.8, 7.5, 8.2, 9.1]
  50.     elif values=='e12':
  51.         r = [1.0,      1.2,      1.5,      1.8,      2.2,      2.7,      3.3,      3.9,      4.7,      5.6,      6.8,      8.2     ]
  52.     elif values=='e6':
  53.         r = [1.0,                1.5,                2.2,                3.3,                4.7,                6.8               ]
  54.     else:
  55.         sys.stderr.write('wrong argument given!\n')
  56.         raise NameError('stop here')
  57.         sys.exit(1)
  58.  
  59.     r = [r[-1]*0.1] + r + [r[0]*10.0]
  60.  
  61.     return r
  62.  
  63.  
  64.  
  65.  
  66. def select_two_nearest(val,val_exp,pvalues_name='e24'):
  67.  
  68.     def my_lt(v1,v2):
  69.         return (v1<v2) and (not math.isclose(v1,v2,rel_tol=1e-3))
  70.  
  71.     def my_gt(v1,v2):
  72.         return (v1>v2) and (not math.isclose(v1,v2,rel_tol=1e-3))
  73.  
  74.  
  75.     pvalues = select_pvalues_list(pvalues_name)
  76.  
  77.     # select range
  78.     if my_lt(val, min(pvalues[1:-1])):
  79.         while my_lt(val, min(pvalues[1:-1])):
  80.             if my_lt(val, max(pvalues[1:-1])*0.1):
  81.                 val     = val     * 10.0;
  82.                 val_exp = val_exp *  0.1;
  83.             else:
  84.                 break;
  85.  
  86.     elif my_gt(val, max(pvalues[1:-1])):
  87.         while my_gt(val, max(pvalues[1:-1])):
  88.             if my_gt(val, min(pvalues[1:-1])*10.0):
  89.                 val     = val     *  0.1;
  90.                 val_exp = val_exp * 10.0;
  91.             else:
  92.                 break;
  93.  
  94.  
  95.     # select nearest values
  96.     for i in range(len(pvalues)-1):
  97.         vdn = pvalues[i];
  98.         vup = pvalues[i+1];
  99.         if vdn<=val and val<=vup:
  100.             break;
  101.  
  102.     vdn_exp = val_exp
  103.     vup_exp = val_exp
  104.  
  105.     if vdn<pvalues[1]:
  106.         vdn *= 10.0
  107.         vdn_exp *= 0.1
  108.  
  109.     if vup>pvalues[-2]:
  110.         vup *= 0.1
  111.         vup_exp *= 10.0
  112.  
  113.     return [ (vdn,vdn_exp), (vup,vup_exp) ];
  114.  
  115. def select_one_nearest(val,pvalues_name='e24'):
  116.  
  117.     lval = select_two_nearest(val,1.0,pvalues_name)
  118.  
  119.     (vdn,vdn_exp) = lval[0]
  120.     (vup,vup_exp) = lval[1]
  121.  
  122.     vdn *= vdn_exp
  123.     vup *= vup_exp
  124.  
  125.     if math.fabs( (vdn-val)/val ) <= math.fabs( (vup-val)/val ):
  126.         return vdn
  127.     else:
  128.         return vup
  129.  
  130. def get_neighbour_value(val,dir='up',pvalues_name='e24'):
  131.    
  132.     pvalues = select_pvalues_list(pvalues_name)
  133.     exp = 1.0
  134.  
  135.     found_v = None
  136.     found_idx = None
  137.     for idx,v in enumerate(pvalues[1:-1]):
  138.         if math.isclose(val,v,rel_tol=1e-3):
  139.             found_v = v
  140.             found_idx = idx
  141.  
  142.             #print(' found {} {}'.format(found_v,found_idx))
  143.             break
  144.    
  145.     if found_v is None:
  146.         raise NameError('found_v is None!')
  147.         sys.exit(1)
  148.  
  149.     if dir=='up':
  150.         if math.isclose(found_v,pvalues[-2],rel_tol=1e-3):
  151.             next_v = pvalues[1]
  152.             exp = 10.0
  153.         else:
  154.             tmp = pvalues[1:-1]
  155.             #print('  tmp {}'.format(tmp))
  156.             next_v = tmp[found_idx+1]
  157.        
  158.     elif dir=='down':
  159.         if math.isclose(found_v,pvalues[1],rel_tol=1e-3):
  160.             next_v = pvalues[-2]
  161.             exp = 0.1
  162.         else:
  163.             tmp = pvalues[1:-1]
  164.             next_v = tmp[found_idx-1]
  165.  
  166.     else:
  167.         raise NameError('dir is neither "up" nor "down"!')
  168.         sys.exit(1)
  169.  
  170.     return (next_v,exp)
  171.  
  172.  
  173. def val_exp_le( val1, val1_exp, val2, val2_exp ):
  174.  
  175.     return (val1_exp < val2_exp) or ( math.isclose(val1_exp, val2_exp, rel_tol=1e-3) and ( (val1 < val2) or math.isclose(val1, val2, rel_tol=1e-3) ) )
  176.  
  177.  
  178. def select_values_in_range(vrange, vrange_exp, pvalues='e24'):
  179. # range is [min,max]
  180. # returns list of values already multiplied by exp
  181.  
  182.     vmin=vrange[0]
  183.     vmax=vrange[1]
  184.  
  185.     if vmax<=vmin:
  186.         raise NameError('wrong vrange: [0] must be bigger than [1]')
  187.         sys.exit(1)
  188.  
  189.     # expand interval by nearest preferred values
  190.     vmin_vals = select_two_nearest(vmin,vrange_exp, pvalues_name=pvalues)
  191.     vmax_vals = select_two_nearest(vmax,vrange_exp, pvalues_name=pvalues)
  192.     #print(vmin,vmax,vrange_exp)
  193.     #print(vmin_vals,vmax_vals)
  194.  
  195.  
  196.     # start building values list
  197.  
  198.     vmin_val_0, vmin_exp_0 = vmin_vals[0]
  199.     vmin_val_1, vmin_exp_1 = vmin_vals[1]
  200.     vmax_val_0, vmax_exp_0 = vmax_vals[0]
  201.     vmax_val_1, vmax_exp_1 = vmax_vals[1]
  202.  
  203.     if math.isclose(vmin_val_1*vmin_exp_1, vmin*vrange_exp, rel_tol=1e-3):
  204.         curr_v, curr_exp = vmin_vals[1]
  205.     else:
  206.         curr_v, curr_exp = vmin_vals[0]
  207.  
  208.     if math.isclose(vmax_val_0*vmax_exp_0, vmax*vrange_exp, rel_tol=1e-3):
  209.         vmax_v, vmax_exp = vmax_vals[0]
  210.     else:
  211.         vmax_v, vmax_exp = vmax_vals[1]
  212.  
  213.     values_range = []
  214.  
  215.     #print('vmaxs: {} {}'.format(vmax_v, vmax_exp))
  216.     while val_exp_le( curr_v, curr_exp, vmax_v, vmax_exp ):
  217.      
  218.         #print('currs: {} {}'.format(curr_v, curr_exp))
  219.  
  220.         values_range = values_range + [ (curr_v, curr_exp) ]
  221.  
  222.         next_v, exp_update = get_neighbour_value( curr_v, dir='up', pvalues_name=pvalues )
  223.        
  224.         #print('nexts: {} {}\n'.format(next_v, exp_update))
  225.  
  226.         curr_v = next_v
  227.         curr_exp *= exp_update
  228.  
  229.    
  230.     return values_range
  231.  
  232.  
  233.  
  234.  
  235. def calc_c1c2r4_from_w1w2q_r1r2r3c3(w1,w2,q,r1,r2,r3,c3):
  236.  
  237.     A = 1.0/w1 + 1.0/(q*w2)      
  238.     B = 1.0/(w2*w2) + 1.0/(q*w1*w2)
  239.     C = 1.0/(w1*w2*w2)          
  240.    
  241. #    print('A={}'.format(A))
  242. #    print('B={}'.format(B))
  243. #    print('C={}'.format(C))
  244.  
  245.     try:
  246.         r4 = (-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3))
  247.         c1 = -c3*(r1 + r2 + r3)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))/(r1*r2) + (A - c3*r3)*(r1 + r2)/(r1*r2)
  248.         c2 = C*(r1 + r2)/(c3*r3*(-c3*(r1 + r2 + r3)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))**2 + (A - c3*r3)*(r1 + r2)*((-(-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))/(3*(math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)) - (math.sqrt(-4*((-3*B*r1*r2*(r1 + r2)*(r1 + r2 + r3) + 6*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) - 3*r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**2/(c3**2*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**2))**3 + ((-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + 2*(c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**2)/2 + (-27*B*r1*r2*(A - c3*r3)*(r1 + r2)**2/c3 + 27*C*r1*r2*(r1 + r2)**2/c3 + 27*r1*r2*r3*(A - c3*r3)**2*(r1 + r2)**2)/(2*c3**2*r1*(r2 + r3)*(r1 + r2 + r3)**2) - (9*c3*r1*r2*r3*(r1 + r2 + r3) - 18*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))*(B*r1*r2*(r1 + r2)*(r1 + r2 + r3) - 2*c3*r1*r2*r3*(A - c3*r3)*(r1 + r2)*(r1 + r2 + r3) + r1*(A - c3*r3)**2*(r1 + r2)**2*(r2 + r3))/(2*c3**3*r1**2*(r2 + r3)**2*(r1 + r2 + r3)**3) + (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))**3/(c3**3*r1**3*(r2 + r3)**3*(r1 + r2 + r3)**3))**(1/3)/3 - (c3*r1*r2*r3*(r1 + r2 + r3) - 2*r1*(A - c3*r3)*(r1 + r2)*(r2 + r3))/(3*c3*r1*(r2 + r3)*(r1 + r2 + r3)))))
  249.     except ValueError:
  250.         return (-1,-1,-1)
  251.  
  252.     return (c1,c2,r4)
  253.  
  254. def remove_imaginary_or_error(x):
  255.  
  256.     if abs(x.real)>1e10*abs(x.imag) and x.real>0:
  257.         return x.real
  258.  
  259.     return -1
  260.  
  261.  
  262. def calc_mfb3_parameters_from_values(r1,r2,r3,r4,c1,c2,c3):
  263.  
  264.     A = (c1*r1*r2 + c3*r1*r3 + c3*r1*r4 + c3*r2*r3 + c3*r2*r4 + c3*r3*r4)/(r1 + r2)
  265.     B = c3*(c1*r1*r2*r3 + c1*r1*r2*r4 + c1*r1*r3*r4 + c2*r1*r3*r4 + c2*r2*r3*r4)/(r1 + r2)
  266.     C = c1*c2*c3*r1*r2*r3*r4/(r1 + r2)
  267.  
  268.     k = -r3/(r1+r2)
  269.    
  270.     w1 = (2**(1/3)*A*C - 2**(1/3)*B**2/3 + B*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)/3 - (18*A*B*C - 4*B**3 - 54*C**2 + 6*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(2/3)/6)/(C*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))
  271.     w2 = 2**(1/6)*cmath.sqrt(3)/cmath.sqrt((3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))
  272.     q  = 2**(5/6)*cmath.sqrt(3)*cmath.sqrt((3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))/(6*(A*(3*2**(2/3)*A*C - 2**(2/3)*B**2 + (2**(1/3)*B - (9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3))*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)) - 3*2**(1/3)*C*(9*A*B*C - 2*B**3 - 27*C**2 + 3*cmath.sqrt(3)*cmath.sqrt(C**2*(4*A**3*C - A**2*B**2 - 18*A*B*C + 4*B**3 + 27*C**2)))**(1/3)))
  273.  
  274.     w1 = remove_imaginary_or_error(w1)
  275.     w2 = remove_imaginary_or_error(w2)
  276.     q  = remove_imaginary_or_error(q)
  277.  
  278.     return (w1,w2,q,k)
  279.  
  280.  
  281. def randomize_mfb3(r1,r2,r3,r4,c1,c2,c3, c_precision=0.05, r_precision=0.01, iters=1000):
  282.  
  283.     (w1_0, w2_0, q_0, k_0) = calc_mfb3_parameters_from_values(r1,r2,r3,r4,c1,c2,c3)
  284.  
  285.     d_w1 = 0.0
  286.     d_w2 = 0.0
  287.     d_q  = 0.0
  288.     d_k  = 0.0
  289.  
  290.     for i in range(iters):
  291.         dr1 = random.uniform(1.0-r_precision,1.0+r_precision)
  292.         dr2 = random.uniform(1.0-r_precision,1.0+r_precision)
  293.         dr3 = random.uniform(1.0-r_precision,1.0+r_precision)
  294.         dr4 = random.uniform(1.0-r_precision,1.0+r_precision)
  295.  
  296.         dc1 = random.uniform(1.0-c_precision,1.0+c_precision)
  297.         dc2 = random.uniform(1.0-c_precision,1.0+c_precision)
  298.         dc3 = random.uniform(1.0-c_precision,1.0+c_precision)
  299.  
  300.         rr1 = r1 * dr1
  301.         rr2 = r2 * dr2
  302.         rr3 = r3 * dr3
  303.         rr4 = r4 * dr4
  304.  
  305.         rc1 = c1 * dc1
  306.         rc2 = c2 * dc2
  307.         rc3 = c3 * dc3
  308.  
  309.         (w1,w2,q,k) = calc_mfb3_parameters_from_values(rr1,rr2,rr3,rr4,rc1,rc2,rc3)
  310.  
  311.         d_w1 = max(d_w1, abs(w1-w1_0))
  312.         d_w2 = max(d_w2, abs(w2-w2_0))
  313.         d_q  = max(d_q,  abs( q- q_0))
  314.         d_k  = max(d_k,  abs( k- k_0))
  315.  
  316.     return (d_w1/w1_0, d_w2/w2_0, d_q/q_0, d_k/k_0)
  317.  
  318.  
  319.  
  320. def calc_filter( w1_set=2*math.pi*30000, w2_set=2*math.pi*30000, q_set=1.0, k_set=(-0.4), r1_range=[1.0,20.0], r1_exp=1e3, r2_range=[1.0,20.0], r2_exp=1e3, r_pvalues='e24', c3_range=[56.0,4700.0], c3_exp=1e-12, c_pvalues='e12'):
  321.  
  322.     c3_values = select_values_in_range( c3_range, c3_exp, pvalues=c_pvalues )
  323.  
  324.     r1_values = select_values_in_range( [x for x in r1_range], r1_exp, pvalues=r_pvalues )
  325.     r2_values = select_values_in_range( [x for x in r2_range], r2_exp, pvalues=r_pvalues )
  326.  
  327.     print( 'c3 values for sweep: {}'.format([x*y for (x,y) in c3_values]) )
  328.     print( 'r1 values for sweep: {}'.format([x*y for (x,y) in r1_values]) )
  329.     print( 'r2 values for sweep: {}'.format([x*y for (x,y) in r2_values]) )
  330.  
  331.     print( '\ncalculating filter for: f1={}, f2={}, q={}, k={}'.format(w1_set/(2*math.pi),w2_set/(2*math.pi),q_set,k_set) )
  332.  
  333.  
  334.     f_values=[]
  335.  
  336.     for r1 in [x*y for (x,y) in r1_values]:
  337.         for r2 in [x*y for (x,y) in r2_values]:
  338.             for c3 in [x*y for (x,y) in c3_values]:
  339.  
  340.                 r3_ideal = (-1) * k_set * (r1+r2)
  341.                 if r3_ideal<=0.0:
  342.                     continue
  343.  
  344.                 r3 = select_one_nearest(r3_ideal, r_pvalues)
  345.  
  346.                 print('\ntry c3={:1.1e}, r1={:.1f}, r2={:.1f}:'.format(c3,r1,r2))
  347.  
  348.                 (c1_ideal,c2_ideal,r4_ideal) = calc_c1c2r4_from_w1w2q_r1r2r3c3(w1_set,w2_set,q_set,r1,r2,r3,c3)
  349.                 if r4_ideal<=0.0 or c1_ideal<=0.0 or c2_ideal<=0.0:
  350.                     print(' FAILURE!')
  351.                     continue
  352.  
  353.  
  354.  
  355.                 print(' ideal values: r3={:.1f}, r4={:.1f}, c1={:1.4e}, c2={:1.4e}'.format(r3, r4_ideal, c1_ideal, c2_ideal))
  356.  
  357.                 r4 = select_one_nearest(r4_ideal, r_pvalues)
  358.                 c1 = select_one_nearest(c1_ideal, c_pvalues)
  359.                 c2 = select_one_nearest(c2_ideal, c_pvalues)
  360.                 print(' nearest values: r3={:.1f}, r4={:.1f}, c1={:1.4e}, c2={:1.4e}'.format(r3, r4, c1, c2))
  361.  
  362.                 (w1_real,w2_real,q_real,k_real) = calc_mfb3_parameters_from_values(r1,r2,r3,r4,c1,c2,c3)
  363.                 if w1_real<0 or w2_real<0 or q_real<0:
  364.                     print(' FAILURE!')
  365.                     continue
  366.  
  367.                 print(' real parameters: f1={:.1f}, f2={:.1f}, q={:.4f}, k={:.4f}'.format(w1_real/(2*math.pi),w2_real/(2*math.pi),q_real,k_real) )
  368.  
  369.                 f_values += [ {'r1':r1,'r2':r2,'r3':r3,'r4':r4, 'c1':c1,'c2':c2,'c3':c3, 'w1':w1_real,'w2':w2_real,'q':q_real,'k':k_real} ]
  370.  
  371.  
  372.    
  373.     #print(f_values)
  374.  
  375.     # sort values by cutoff freq match
  376.     print('\n\n\n sorted by neares cutoff freqs:')
  377.     for e in  sorted(f_values, key = lambda tup: (tup['w1']-w1_set)**2/(w1_set*w1_set) + (tup['w2']-w2_set)**2/(w2_set*w2_set) ):
  378.         print('  r1={:.1f}, r2={:.1f}, r3={:.1f}, r4={:.1f}, c1={:1.1e}, c2={:1.1e}, c3={:1.1e}, f1={:.1f}, f2={:.1f}, q={:.4f}, k={:.4f}'.format(e['r1'],e['r2'],e['r3'],e['r4'],e['c1'],e['c2'],e['c3'],e['w1']/(2*math.pi),e['w2']/(2*math.pi),e['q'],e['k']) )
  379.         (p_w1, p_w2, p_q, p_k) = randomize_mfb3(e['r1'],e['r2'],e['r3'],e['r4'], e['c1'],e['c2'],e['c3'])
  380.         e['p_w1']=p_w1
  381.         e['p_w2']=p_w2
  382.         e['p_q'] =p_q
  383.         e['p_k'] =p_k
  384.  
  385.         print('   percent_f1={:.2f}%, percent_f2={:.2f}%, percent_q={:.2f}%, percent_k={:.2f}%'.format(100*e['p_w1'],100*e['p_w2'],100*e['p_q'],100*e['p_k']))
  386.  
  387.     # sort values by q match
  388.     print('\n\n\n sorted by q:')
  389.     for e in  sorted(f_values, key = lambda tup: abs((tup['q']-q_set)/q_set) ):
  390.         print('  r1={:.1f}, r2={:.1f}, r3={:.1f}, r4={:.1f}, c1={:1.1e}, c2={:1.1e}, c3={:1.1e}, f1={:.1f}, f2={:.1f}, q={:.4f}, k={:.4f}'.format(e['r1'],e['r2'],e['r3'],e['r4'],e['c1'],e['c2'],e['c3'],e['w1']/(2*math.pi),e['w2']/(2*math.pi),e['q'],e['k']) )
  391.         print('   percent_f1={:.2f}%, percent_f2={:.2f}%, percent_q={:.2f}%, percent_k={:.2f}%'.format(100*e['p_w1'],100*e['p_w2'],100*e['p_q'],100*e['p_k']))
  392.  
  393.  
  394.  
  395. def main():
  396.  
  397.     calc_filter()
  398.  
  399.     #print(" #### AY ####")
  400.     #calc_filter( wc_set=2*math.pi*30000, q_set=1.0, k_set=(-1.75), r1_range=[20.0,30.0], r1_exp=1e3, r1_div=3.0, r_pvalues='e24', c1_range=[33.0,4700.0], c1_exp=1e-12, c_pvalues='e12')
  401.  
  402. #    print(" #### beep ####")
  403. #    calc_filter( wc_set=2*math.pi*30000, q_set=1.0, k_set=(-0.388), r1_range=[4.7,10.0], r1_exp=1e3, r1_div=1.0, r1_add=750, r_pvalues='e24', c1_range=[33.0,4700.0], c1_exp=1e-12, c_pvalues='e12')
  404.  
  405.  
  406. if __name__=="__main__":
  407.     main()
  408.  
  409.  
  410.  
  411.