diff --git a/explore/angles.py b/explore/angles.py index cb972466..a19ab9a5 100755 --- a/explore/angles.py +++ b/explore/angles.py @@ -223,8 +223,8 @@ def print_steps(jitter, jitter_inv, view, view_inv, qc_only): ) comment("\n**** qab from qc ****") - # The indirect calculation of qab is better than directly c - # alculating qab^2 = qa^2 + qb^2 since qc can be computed + # The indirect calculation of qab is better than directly + # calculating qab^2 = qa^2 + qb^2 since qc can be computed # as qc = M31*qx + M32*qy, thus requiring only two elements # of the rotation matrix. #vprint(qab, sqrt(qa**2 + qb**2), "Direct calculation of qab") diff --git a/sasmodels/TwoYukawa/CalTYSk.py b/sasmodels/TwoYukawa/CalTYSk.py index 40f50e0c..0714c7b2 100644 --- a/sasmodels/TwoYukawa/CalTYSk.py +++ b/sasmodels/TwoYukawa/CalTYSk.py @@ -1,8 +1,8 @@ import numpy as np -from numpy import pi, mean +from numpy import mean, pi -from .Ecoefficient import TYCoeff from .CalcRealRoot import CalRealRoot +from .Ecoefficient import TYCoeff from .TInvFourier import TInvFourier # Supplied Q vector must go out to at least this value to calculate g(r). diff --git a/sasmodels/TwoYukawa/CalcRealRoot.py b/sasmodels/TwoYukawa/CalcRealRoot.py index 9db225ed..7e7287a5 100644 --- a/sasmodels/TwoYukawa/CalcRealRoot.py +++ b/sasmodels/TwoYukawa/CalcRealRoot.py @@ -3,6 +3,7 @@ from .Ecoefficient import TYCoeff from .Epoly import make_poly + def CalRealRoot(coeff: TYCoeff): Ecoefficient = coeff.Ecoefficient diff --git a/sasmodels/TwoYukawa/Ecoefficient.py b/sasmodels/TwoYukawa/Ecoefficient.py index 3182c01a..9d58f60c 100644 --- a/sasmodels/TwoYukawa/Ecoefficient.py +++ b/sasmodels/TwoYukawa/Ecoefficient.py @@ -1,7 +1,7 @@ from typing import Callable, Tuple import numpy as np -from numpy import exp, pi, cos, sin, cosh +from numpy import cos, cosh, exp, pi, sin from numpy.typing import NDArray # CalCoeff.m @@ -194,17 +194,17 @@ def ABC12(d1, d2): Ccd2_22*Cdd1_12*d1*d2**2 - Ccd1_21*Cdd2_22*d1*d2**2 - Ccd2_22*d2*k1 + Ccd1_21*d1*k2)))/ - ((d1*(( + (d1*( (-Ccd1_21)*Ccd2_12*d2 + - Ccd1_11*Ccd2_22*d2))))) + Ccd1_11*Ccd2_22*d2))) C2 = ((-((Ccd2_12*d2* (((-Cd1_1)*d1 - Cdd1_11*d1**2 - Cdd1_12*d1*d2 + k1)) - Ccd1_11*d1* (((-Cd2_2)*d2 - Cdd2_12*d1*d2 - Cdd2_22*d2**2 + k2))))) / - (((-Ccd1_21)*Ccd2_12*d1*d2 + - Ccd1_11*Ccd2_22*d1*d2))) + ((-Ccd1_21)*Ccd2_12*d1*d2 + + Ccd1_11*Ccd2_22*d1*d2)) return A, B, C1, C2 self.ABC12 = ABC12 @@ -370,10 +370,10 @@ def tau_d21(s): E1d02 = 12*c1F01*phi*sigma_d01(z1) - 12*c1F01*exp(-z1)*phi*tau_d01(z1) - E1d11 = (((12*c1F10*phi*sigma_d01(z1) + \ + E1d11 = (12*c1F10*phi*sigma_d01(z1) + \ 12*c1F01*phi*sigma_d10(z1) - \ 12*c1F10*exp(-z1)*phi*tau_d01(z1) - \ - 12*c1F01*exp(-z1)*phi*tau_d10(z1)))) + 12*c1F01*exp(-z1)*phi*tau_d10(z1)) E1d12 = (-c1F01 + 12*c1F11*phi*sigma_d01(z1) + 12*c1F01*phi*sigma_d11(z1) - 12*c1F11*exp(-z1)*phi*tau_d01(z1) - diff --git a/sasmodels/TwoYukawa/Epoly.py b/sasmodels/TwoYukawa/Epoly.py index 274bbe00..f5260f91 100644 --- a/sasmodels/TwoYukawa/Epoly.py +++ b/sasmodels/TwoYukawa/Epoly.py @@ -1,6 +1,7 @@ import numpy as np from numpy.typing import NDArray + #from numba import njit #@njit(cache=True) def make_poly(Ecoefficient:NDArray) -> NDArray: @@ -432,8 +433,8 @@ def make_poly(Ecoefficient:NDArray) -> NDArray: gE1d02*gE1d11*gE2d23*gE2d33**6) polyd2_19 = (7*gE1d02**2*gE2d31*gE2d33**6 + - 3*gE1d13**2*gE2d31*gE2d33**4*((5*gE2d13*gE2d31 + 2*gE2d11*gE2d33)) + - gE2d33*(((-9)*gE1d12*gE1d42*gE2d13*gE2d23*gE2d24**2*gE2d31**2 - + 3*gE1d13**2*gE2d31*gE2d33**4*(5*gE2d13*gE2d31 + 2*gE2d11*gE2d33) + + gE2d33*((-9)*gE1d12*gE1d42*gE2d13*gE2d23*gE2d24**2*gE2d31**2 - 3*gE1d11*gE1d42*gE2d13*gE2d24**3*gE2d31**2 + 3*gE1d42**2*gE2d13**4*gE2d31*gE2d33 - 3*gE1d32*gE1d42*gE2d13**3*gE2d23*gE2d31*gE2d33 - @@ -498,38 +499,38 @@ def make_poly(Ecoefficient:NDArray) -> NDArray: 2*gE1d11*gE1d31*gE2d13**2*gE2d33**4 + 6*gE1d12**2*gE2d13*gE2d31*gE2d33**4 + gE1d12**2*gE2d11*gE2d33**5 + gE1d11**2*gE2d13*gE2d33**5 + - gE1d22**2*gE2d13*gE2d33**3*((5*gE2d13*gE2d31 + - 2*gE2d11*gE2d33)) + - 3*gE1d33**2*gE2d13*gE2d33*((2*gE2d13**2*gE2d31**2 + - 4*gE2d11*gE2d13*gE2d31*gE2d33 + gE2d11**2*gE2d33**2)) - - gE1d22*((gE2d33**2*((gE1d32*gE2d13*((4*gE2d13*gE2d23*gE2d31 + - gE2d13*gE2d21*gE2d33 + 2*gE2d11*gE2d23*gE2d33)) + - gE1d31*gE2d13*((4*gE2d13*gE2d24*gE2d31 + - gE2d13*gE2d22*gE2d33 + 2*gE2d11*gE2d24*gE2d33)) + - gE2d33*((5*gE1d12*gE2d13*gE2d23*gE2d31 + + gE1d22**2*gE2d13*gE2d33**3*(5*gE2d13*gE2d31 + + 2*gE2d11*gE2d33) + + 3*gE1d33**2*gE2d13*gE2d33*(2*gE2d13**2*gE2d31**2 + + 4*gE2d11*gE2d13*gE2d31*gE2d33 + gE2d11**2*gE2d33**2) - + gE1d22*(gE2d33**2*(gE1d32*gE2d13*(4*gE2d13*gE2d23*gE2d31 + + gE2d13*gE2d21*gE2d33 + 2*gE2d11*gE2d23*gE2d33) + + gE1d31*gE2d13*(4*gE2d13*gE2d24*gE2d31 + + gE2d13*gE2d22*gE2d33 + 2*gE2d11*gE2d24*gE2d33) + + gE2d33*(5*gE1d12*gE2d13*gE2d23*gE2d31 + 5*gE1d11*gE2d13*gE2d24*gE2d31 + gE1d12*gE2d13*gE2d21*gE2d33 + gE1d11*gE2d13*gE2d22*gE2d33 + gE1d12*gE2d11*gE2d23*gE2d33 + - gE1d11*gE2d11*gE2d24*gE2d33)))) + - gE1d42*(((-gE2d11**2)*gE2d24**2*gE2d33**2 + + gE1d11*gE2d11*gE2d24*gE2d33)) + + gE1d42*((-gE2d11**2)*gE2d24**2*gE2d33**2 + 8*gE2d13**3*gE2d31*gE2d33**2 - - 2*gE2d11*gE2d13*gE2d33*((3*gE2d24**2*gE2d31 + - gE2d23**2*gE2d33 + 2*gE2d22*gE2d24*gE2d33)) - - gE2d13**2*((3*gE2d24**2*gE2d31**2 + - 2*gE2d24*gE2d33*((3*gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d33*((3*gE2d23**2*gE2d31 + + 2*gE2d11*gE2d13*gE2d33*(3*gE2d24**2*gE2d31 + + gE2d23**2*gE2d33 + 2*gE2d22*gE2d24*gE2d33) - + gE2d13**2*(3*gE2d24**2*gE2d31**2 + + 2*gE2d24*gE2d33*(3*gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d33*(3*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - - 6*gE2d11*gE2d33)))))))))) + - gE1d33*(((-gE1d42)*gE2d13*((3*gE2d11**2*gE2d24*gE2d33**2 + - 3*gE2d11*gE2d13*gE2d33*((3*gE2d24*gE2d31 + - gE2d22*gE2d33)) + - gE2d13**2*((3*gE2d24*gE2d31**2 + - gE2d33*((3*gE2d22*gE2d31 + - gE2d20*gE2d33)))))) + - gE2d33*((6*gE1d11*gE2d13*gE2d24**2*gE2d31**2 + + gE2d33*(gE2d22**2 - + 6*gE2d11*gE2d33))))) + + gE1d33*((-gE1d42)*gE2d13*(3*gE2d11**2*gE2d24*gE2d33**2 + + 3*gE2d11*gE2d13*gE2d33*(3*gE2d24*gE2d31 + + gE2d22*gE2d33) + + gE2d13**2*(3*gE2d24*gE2d31**2 + + gE2d33*(3*gE2d22*gE2d31 + + gE2d20*gE2d33))) + + gE2d33*(6*gE1d11*gE2d13*gE2d24**2*gE2d31**2 + 8*gE1d31*gE2d13**3*gE2d31*gE2d33 + 4*gE1d11*gE2d13*gE2d23**2*gE2d31*gE2d33 + 8*gE1d11*gE2d13*gE2d22*gE2d24*gE2d31*gE2d33 + @@ -542,148 +543,148 @@ def make_poly(Ecoefficient:NDArray) -> NDArray: 2*gE1d11*gE2d11*gE2d22*gE2d24*gE2d33**2 - 10*gE1d11*gE2d13**2*gE2d31*gE2d33**2 - 4*gE1d11*gE2d11*gE2d13*gE2d33**3 - - gE1d22*((gE2d11**2*gE2d24*gE2d33**2 + - 2*gE2d11*gE2d13*gE2d33*((4*gE2d24*gE2d31 + - gE2d22*gE2d33)) + - gE2d13**2*((6*gE2d24*gE2d31**2 + + gE1d22*(gE2d11**2*gE2d24*gE2d33**2 + + 2*gE2d11*gE2d13*gE2d33*(4*gE2d24*gE2d31 + + gE2d22*gE2d33) + + gE2d13**2*(6*gE2d24*gE2d31**2 + 4*gE2d22*gE2d31*gE2d33 + - gE2d20*gE2d33**2)))) + - 2*gE1d12*((gE2d11*gE2d33*((4*gE2d23*gE2d24*gE2d31 + + gE2d20*gE2d33**2)) + + 2*gE1d12*(gE2d11*gE2d33*(4*gE2d23*gE2d24*gE2d31 + gE2d22*gE2d23*gE2d33 + - gE2d21*gE2d24*gE2d33)) + - gE2d13*((gE2d21*gE2d33*((4*gE2d24*gE2d31 + - gE2d22*gE2d33)) + - gE2d23*((6*gE2d24*gE2d31**2 + + gE2d21*gE2d24*gE2d33) + + gE2d13*(gE2d21*gE2d33*(4*gE2d24*gE2d31 + + gE2d22*gE2d33) + + gE2d23*(6*gE2d24*gE2d31**2 + 4*gE2d22*gE2d31*gE2d33 + - gE2d20*gE2d33**2)))))))))))) + - gE1d13*((gE1d42*((3*gE2d13**2*gE2d33**2*((6*gE2d24*gE2d31**2 + - gE2d33*((4*gE2d22*gE2d31 + gE2d20*gE2d33)))) + - 3*gE2d11*gE2d33*(((-gE2d24**3)*gE2d31**2 - + gE2d20*gE2d33**2)))))) + + gE1d13*(gE1d42*(3*gE2d13**2*gE2d33**2*(6*gE2d24*gE2d31**2 + + gE2d33*(4*gE2d22*gE2d31 + gE2d20*gE2d33)) + + 3*gE2d11*gE2d33*((-gE2d24**3)*gE2d31**2 - gE2d22*gE2d23**2*gE2d33**2 - - gE2d24**2*gE2d33*((3*gE2d22*gE2d31 + gE2d20*gE2d33)) - - gE2d24*gE2d33*((3*gE2d23**2*gE2d31 + + gE2d24**2*gE2d33*(3*gE2d22*gE2d31 + gE2d20*gE2d33) - + gE2d24*gE2d33*(3*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - gE2d11*gE2d33)))))) - - gE2d13*((gE2d24**3*gE2d31**3 + - 9*gE2d24**2*gE2d31*gE2d33*((gE2d22*gE2d31 + - gE2d20*gE2d33)) + - 3*gE2d24*gE2d33*((3*gE2d23**2*gE2d31**2 + + gE2d33*(gE2d22**2 - gE2d11*gE2d33))) - + gE2d13*(gE2d24**3*gE2d31**3 + + 9*gE2d24**2*gE2d31*gE2d33*(gE2d22*gE2d31 + + gE2d20*gE2d33) + + 3*gE2d24*gE2d33*(3*gE2d23**2*gE2d31**2 + 6*gE2d21*gE2d23*gE2d31*gE2d33 + - gE2d33*((3*gE2d22**2*gE2d31 + gE2d21**2*gE2d33 + + gE2d33*(3*gE2d22**2*gE2d31 + gE2d21**2*gE2d33 + 2*gE2d20*gE2d22*gE2d33 - - 8*gE2d11*gE2d31*gE2d33)))) + - gE2d33**2*((gE2d22**3*gE2d33 + + 8*gE2d11*gE2d31*gE2d33)) + + gE2d33**2*(gE2d22**3*gE2d33 + 3*gE2d20*gE2d23**2*gE2d33 + - gE2d22*((9*gE2d23**2*gE2d31 + + gE2d22*(9*gE2d23**2*gE2d31 + 6*gE2d21*gE2d23*gE2d33 - - 6*gE2d11*gE2d33**2)))))))) + - gE2d33*((gE1d33*(((-20)*gE2d13**2*gE2d31**2*gE2d33**2 + - gE2d11*gE2d33*((6*gE2d24**2*gE2d31**2 + - 2*gE2d24*gE2d33*((4*gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d33*((4*gE2d23**2*gE2d31 + + 6*gE2d11*gE2d33**2)))) + + gE2d33*(gE1d33*((-20)*gE2d13**2*gE2d31**2*gE2d33**2 + + gE2d11*gE2d33*(6*gE2d24**2*gE2d31**2 + + 2*gE2d24*gE2d33*(4*gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d33*(4*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - - 2*gE2d11*gE2d33)))))) + - gE2d13*((4*gE2d24**2*gE2d31**3 + - 4*gE2d24*gE2d31*gE2d33*((3*gE2d22*gE2d31 + - 2*gE2d20*gE2d33)) + - gE2d33*((6*gE2d23**2*gE2d31**2 + + gE2d33*(gE2d22**2 - + 2*gE2d11*gE2d33))) + + gE2d13*(4*gE2d24**2*gE2d31**3 + + 4*gE2d24*gE2d31*gE2d33*(3*gE2d22*gE2d31 + + 2*gE2d20*gE2d33) + + gE2d33*(6*gE2d23**2*gE2d31**2 + 8*gE2d21*gE2d23*gE2d31*gE2d33 + - gE2d33*((4*gE2d22**2*gE2d31 + + gE2d33*(4*gE2d22**2*gE2d31 + gE2d21**2*gE2d33 + 2*gE2d20*gE2d22*gE2d33 - - 20*gE2d11*gE2d31*gE2d33)))))))) + - gE2d33*((2*gE1d32*((gE2d11*gE2d33*((4*gE2d23*gE2d24* -gE2d31 + gE2d22*gE2d23*gE2d33 + gE2d21*gE2d24*gE2d33)) + - gE2d13*((gE2d21*gE2d33*((4*gE2d24*gE2d31 + - gE2d22*gE2d33)) + - gE2d23*((6*gE2d24*gE2d31**2 + + 20*gE2d11*gE2d31*gE2d33)))) + + gE2d33*(2*gE1d32*(gE2d11*gE2d33*(4*gE2d23*gE2d24* +gE2d31 + gE2d22*gE2d23*gE2d33 + gE2d21*gE2d24*gE2d33) + + gE2d13*(gE2d21*gE2d33*(4*gE2d24*gE2d31 + + gE2d22*gE2d33) + + gE2d23*(6*gE2d24*gE2d31**2 + 4*gE2d22*gE2d31*gE2d33 + - gE2d20*gE2d33**2)))))) + - gE2d33*(((-gE1d22)*((gE2d11*gE2d33*((5*gE2d24* -gE2d31 + gE2d22*gE2d33)) + gE2d13*((10*gE2d24*gE2d31**2 + + gE2d20*gE2d33**2))) + + gE2d33*((-gE1d22)*(gE2d11*gE2d33*(5*gE2d24* +gE2d31 + gE2d22*gE2d33) + gE2d13*(10*gE2d24*gE2d31**2 + 5*gE2d22*gE2d31*gE2d33 + - gE2d20*gE2d33**2)))) - - gE2d33*(((-2)*gE1d11*gE2d33*((6*gE2d13* -gE2d31 + gE2d11*gE2d33)) + gE1d02*((15*gE2d24*gE2d31**2 + + gE2d20*gE2d33**2)) - + gE2d33*((-2)*gE1d11*gE2d33*(6*gE2d13* +gE2d31 + gE2d11*gE2d33) + gE1d02*(15*gE2d24*gE2d31**2 + 6*gE2d22*gE2d31*gE2d33 + - gE2d20*gE2d33**2)))))) + - gE1d31*(((-10)*gE2d13**2*gE2d31*gE2d33**2 + - gE2d11*gE2d33*((4*gE2d24**2*gE2d31 + + gE2d20*gE2d33**2))) + + gE1d31*((-10)*gE2d13**2*gE2d31*gE2d33**2 + + gE2d11*gE2d33*(4*gE2d24**2*gE2d31 + gE2d23**2*gE2d33 + - 2*gE2d22*gE2d24*gE2d33)) + - gE2d13*((6*gE2d24**2*gE2d31**2 + - 2*gE2d24*gE2d33*((4*gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d33*((4*gE2d23**2*gE2d31 + + 2*gE2d22*gE2d24*gE2d33) + + gE2d13*(6*gE2d24**2*gE2d31**2 + + 2*gE2d24*gE2d33*(4*gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d33*(4*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - - 4*gE2d11*gE2d33)))))))))))))) -+ gE1d02*((gE1d42*((gE2d24**4*gE2d31**3 + - 12*gE2d24**3*gE2d31*gE2d33*((gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d33**2*((3*gE2d23**4*gE2d31 + 4*gE2d21*gE2d23**3*gE2d33 - + gE2d33*(gE2d22**2 - + 4*gE2d11*gE2d33))))))) ++ gE1d02*(gE1d42*(gE2d24**4*gE2d31**3 + + 12*gE2d24**3*gE2d31*gE2d33*(gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d33**2*(3*gE2d23**4*gE2d31 + 4*gE2d21*gE2d23**3*gE2d33 - 8*gE2d13*gE2d21*gE2d23*gE2d33**2 + - 2*gE2d13*gE2d33**2*(((-2)*gE2d22**2 + - 5*gE2d13*gE2d31 + 2*gE2d11*gE2d33)) - - 2*gE2d23**2*gE2d33*(((-3)*gE2d22**2 + - 8*gE2d13*gE2d31 + 2*gE2d11*gE2d33)))) + - 2*gE2d24**2*gE2d33*((9*gE2d23**2*gE2d31**2 + + 2*gE2d13*gE2d33**2*((-2)*gE2d22**2 + + 5*gE2d13*gE2d31 + 2*gE2d11*gE2d33) - + 2*gE2d23**2*gE2d33*((-3)*gE2d22**2 + + 8*gE2d13*gE2d31 + 2*gE2d11*gE2d33)) + + 2*gE2d24**2*gE2d33*(9*gE2d23**2*gE2d31**2 + 18*gE2d21*gE2d23*gE2d31*gE2d33 + - gE2d33*((9*gE2d22**2*gE2d31 - 12*gE2d13*gE2d31**2 + + gE2d33*(9*gE2d22**2*gE2d31 - 12*gE2d13*gE2d31**2 + 3*gE2d21**2*gE2d33 + 6*gE2d20*gE2d22*gE2d33 - - 8*gE2d11*gE2d31*gE2d33)))) + - 4*gE2d24*gE2d33**2*((gE2d22**3*gE2d33 + - gE2d20*gE2d33*((3*gE2d23**2 - 2*gE2d13*gE2d33)) + - gE2d22*((9*gE2d23**2*gE2d31 + 6*gE2d21*gE2d23*gE2d33 - - 2*gE2d33*((4*gE2d13*gE2d31 + - gE2d11*gE2d33)))))))) - - gE2d33*((gE1d33*((4*gE2d24**3*gE2d31**3 + - 6*gE2d24**2*gE2d31*gE2d33*((3*gE2d22*gE2d31 + - 2*gE2d20*gE2d33)) + - 3*gE2d24*gE2d33*((6*gE2d23**2*gE2d31**2 + + 8*gE2d11*gE2d31*gE2d33)) + + 4*gE2d24*gE2d33**2*(gE2d22**3*gE2d33 + + gE2d20*gE2d33*(3*gE2d23**2 - 2*gE2d13*gE2d33) + + gE2d22*(9*gE2d23**2*gE2d31 + 6*gE2d21*gE2d23*gE2d33 - + 2*gE2d33*(4*gE2d13*gE2d31 + + gE2d11*gE2d33)))) - + gE2d33*(gE1d33*(4*gE2d24**3*gE2d31**3 + + 6*gE2d24**2*gE2d31*gE2d33*(3*gE2d22*gE2d31 + + 2*gE2d20*gE2d33) + + 3*gE2d24*gE2d33*(6*gE2d23**2*gE2d31**2 + 8*gE2d21*gE2d23*gE2d31*gE2d33 + - gE2d33*((4*gE2d22**2*gE2d31 - + gE2d33*(4*gE2d22**2*gE2d31 - 10*gE2d13*gE2d31**2 + gE2d21**2*gE2d33 + 2*gE2d20*gE2d22*gE2d33 - - 5*gE2d11*gE2d31*gE2d33)))) + - gE2d33**2*((gE2d22**3*gE2d33 + - 3*gE2d20*gE2d33*((gE2d23**2 - gE2d13*gE2d33)) + - 3*gE2d22*((4*gE2d23**2*gE2d31 + + 5*gE2d11*gE2d31*gE2d33)) + + gE2d33**2*(gE2d22**3*gE2d33 + + 3*gE2d20*gE2d33*(gE2d23**2 - gE2d13*gE2d33) + + 3*gE2d22*(4*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 - - gE2d33*((5*gE2d13*gE2d31 + - gE2d11*gE2d33)))))))) + - gE2d33*((3*gE1d31*((2*gE2d24**3*gE2d31**2 + - gE2d22*gE2d33**2*((gE2d23**2 - gE2d13*gE2d33)) + - gE2d24**2*gE2d33*((4*gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d24*gE2d33*((4*gE2d23**2*gE2d31 + + gE2d33*(5*gE2d13*gE2d31 + + gE2d11*gE2d33)))) + + gE2d33*(3*gE1d31*(2*gE2d24**3*gE2d31**2 + + gE2d22*gE2d33**2*(gE2d23**2 - gE2d13*gE2d33) + + gE2d24**2*gE2d33*(4*gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d24*gE2d33*(4*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - 5*gE2d13*gE2d31 - - gE2d11*gE2d33)))))) + - gE1d32*((4*gE2d23**3*gE2d31*gE2d33 + + gE2d33*(gE2d22**2 - 5*gE2d13*gE2d31 - + gE2d11*gE2d33))) + + gE1d32*(4*gE2d23**3*gE2d31*gE2d33 + 3*gE2d21*gE2d23**2*gE2d33**2 + - 3*gE2d21*gE2d33*((4*gE2d24**2*gE2d31 + + 3*gE2d21*gE2d33*(4*gE2d24**2*gE2d31 + 2*gE2d22*gE2d24*gE2d33 - - gE2d13*gE2d33**2)) + - 3*gE2d23*((6*gE2d24**2*gE2d31**2 + - gE2d33**2*((gE2d22**2 - 5*gE2d13*gE2d31 - - gE2d11*gE2d33)) + - 2*gE2d24*gE2d33*((4*gE2d22*gE2d31 + - gE2d20*gE2d33)))))) + - gE2d33*((gE2d33**2*((6*gE1d12*gE2d23*gE2d31 + + gE2d13*gE2d33**2) + + 3*gE2d23*(6*gE2d24**2*gE2d31**2 + + gE2d33**2*(gE2d22**2 - 5*gE2d13*gE2d31 - + gE2d11*gE2d33) + + 2*gE2d24*gE2d33*(4*gE2d22*gE2d31 + + gE2d20*gE2d33))) + + gE2d33*(gE2d33**2*(6*gE1d12*gE2d23*gE2d31 + 6*gE1d11*gE2d24*gE2d31 + gE1d12*gE2d21*gE2d33 + - gE1d11*gE2d22*gE2d33)) - - gE1d22*((10*gE2d24**2*gE2d31**2 + - 2*gE2d24*gE2d33*((5*gE2d22*gE2d31 + - gE2d20*gE2d33)) + - gE2d33*((5*gE2d23**2*gE2d31 + + gE1d11*gE2d22*gE2d33) - + gE1d22*(10*gE2d24**2*gE2d31**2 + + 2*gE2d24*gE2d33*(5*gE2d22*gE2d31 + + gE2d20*gE2d33) + + gE2d33*(5*gE2d23**2*gE2d31 + 2*gE2d21*gE2d23*gE2d33 + - gE2d33*((gE2d22**2 - + gE2d33*(gE2d22**2 - 12*gE2d13*gE2d31 - - 2*gE2d11*gE2d33))))))))))))))) + 2*gE2d11*gE2d33)))))))) polyd2_18 = ((-3)*gE1d13*gE1d42*gE2d13*gE2d23*gE2d24**2*gE2d31**3 - gE1d12*gE1d42*gE2d13*gE2d24**3*gE2d31**3 + diff --git a/sasmodels/TwoYukawa/TFourier.py b/sasmodels/TwoYukawa/TFourier.py index 7eda2947..05a935e6 100644 --- a/sasmodels/TwoYukawa/TFourier.py +++ b/sasmodels/TwoYukawa/TFourier.py @@ -1,6 +1,7 @@ -from numpy import exp, pi, arange, linspace, abs, ceil, log2, interp +from numpy import abs, arange, ceil, exp, interp, linspace, log2, pi from numpy.fft import fft + def TFourier(x, y, deltaX): """ Compute the Fourier transform of a function y(x) with sampling interval deltaX. diff --git a/sasmodels/TwoYukawa/TInvFourier.py b/sasmodels/TwoYukawa/TInvFourier.py index 21829935..20978e1c 100644 --- a/sasmodels/TwoYukawa/TInvFourier.py +++ b/sasmodels/TwoYukawa/TInvFourier.py @@ -1,6 +1,7 @@ -from numpy import exp, ceil, log2, pi, arange, interp +from numpy import arange, ceil, exp, interp, log2, pi from numpy.fft import fft + def TInvFourier(x, y, deltaX): """ Inverse Fourier transform implementation diff --git a/sasmodels/jitter.py b/sasmodels/jitter.py index caadff3d..f7dece02 100755 --- a/sasmodels/jitter.py +++ b/sasmodels/jitter.py @@ -46,6 +46,7 @@ """ import argparse +import warnings import numpy as np from numpy import arccos, arctan2, cos, degrees, exp, log, pi, radians, sin, sqrt @@ -569,8 +570,7 @@ def draw_labels(axes, view, jitter, text): # TODO: zdir for labels is broken, and labels aren't appearing. for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): - zdir = np.asarray(zdir).flatten() - axes.text(p[0], p[1], p[2], label, zdir=zdir) + axes.text(p[0], p[1], p[2], label, zdir=tuple(zdir)) # Definition of rotation matrices comes from wikipedia: # https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations @@ -615,15 +615,17 @@ def apply_jitter(jitter, points): Apply the jitter transform to a set of points. Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. + + Jitter is (dtheta, dphi, dpsi, convention) where convention is "zyz" or "xyz". + If convention is missing it defaults to "xyz". """ if jitter is None: return points - # Hack to deal with the fact that azimuthal_equidistance uses euler angles - if len(jitter) == 4: - dtheta, dphi, dpsi, _ = jitter + convention = jitter[3] if len(jitter) == 4 else "xyz" + dtheta, dphi, dpsi = jitter[:3] + if convention == "zyz": points = Rz(dphi)@Ry(dtheta)@Rz(dpsi)@points else: - dtheta, dphi, dpsi = jitter points = Rx(dphi)@Ry(dtheta)@Rz(dpsi)@points return points @@ -808,12 +810,15 @@ def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): ## increase pd_n for testing jitter integration rather than simple viz #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] - pars = dict( - theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n, - phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n, - psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n, - ) - pars.update(calculator.pars) + # Some models do not accept theta, phi or psi + pars = calculator.pars.copy() + if 'theta' in pars: + pars.update(theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n) + if 'phi' in pars: + pars.update(phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n) + if 'psi' in pars: + pars.update(psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n) + # compute the pattern qx, qy = calculator.qxqy @@ -849,7 +854,7 @@ def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): else: axes.pcolormesh(qx, qy, Iqxy) -def build_model(model_name, n=150, qmax=0.5, **pars): +def build_model(model_name, n=150, num_angles=3, qmax=0.5, **pars): """ Build a calculator for the given shape. @@ -881,10 +886,15 @@ def build_model(model_name, n=150, qmax=0.5, **pars): # stuff the values for non-orientation parameters into the calculator calculator.pars = pars.copy() calculator.pars.setdefault('background', 1e-3) + if num_angles > 0: + calculator.pars.update(theta=0, phi=0) + if num_angles == 3: + calculator.pars.update(psi=0) + # fix the data limits so that we can see if the pattern fades # under rotation or angular dispersion - Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) + Iqxy = calculator(**calculator.pars) Iqxy = log(Iqxy) vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') calculator.limits = vmin, vmax+1 @@ -908,52 +918,54 @@ def select_calculator(model_name, n=150, size=(10, 40, 100)): d_factor = 0.06 # for paracrystal models if model_name == 'sphere': calculator = build_model( - 'sphere', n=n, radius=c) + 'sphere', n=n, num_angles=0, radius=c) a = b = c elif model_name == 'sc_paracrystal': a = b = c dnn = c radius = 0.5*c calculator = build_model( - 'sc_paracrystal', n=n, + 'sc_paracrystal', n=n, num_angles=0, dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius, background=0) elif model_name == 'fcc_paracrystal': a = b = c # nearest neigbour distance dnn should be 2 radius, but I think the # model uses lattice spacing rather than dnn in its calculations + warnings.warn("nearest neighbour distance may be wrong here") dnn = 0.5*c radius = sqrt(2)/4 * c calculator = build_model( - 'fcc_paracrystal', n=n, + 'fcc_paracrystal', n=n, num_angles=0, dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius, background=0) elif model_name == 'bcc_paracrystal': a = b = c # nearest neigbour distance dnn should be 2 radius, but I think the # model uses lattice spacing rather than dnn in its calculations + warnings.warn("nearest neighbour distance may be wrong here") dnn = 0.5*c radius = sqrt(3)/2 * c calculator = build_model( - 'bcc_paracrystal', n=n, + 'bcc_paracrystal', n=n, num_angles=0, dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius, background=0) elif model_name == 'cylinder': calculator = build_model( - 'cylinder', n=n, qmax=0.3, radius=b, length=c) + 'cylinder', n=n, num_angles=2, qmax=0.3, radius=b, length=c) a = b elif model_name == 'ellipsoid': calculator = build_model( - 'ellipsoid', n=n, qmax=1.0, + 'ellipsoid', n=n, num_angles=2, qmax=1.0, radius_polar=c, radius_equatorial=b) a = b elif model_name == 'triaxial_ellipsoid': calculator = build_model( - 'triaxial_ellipsoid', n=n, qmax=0.5, + 'triaxial_ellipsoid', n=n, num_angles=3, qmax=0.5, radius_equat_minor=a, radius_equat_major=b, radius_polar=c) elif model_name == 'parallelepiped': calculator = build_model( - 'parallelepiped', n=n, length_a=a, length_b=b, length_c=c) + 'parallelepiped', n=n, num_angles=3, length_a=a, length_b=b, length_c=c) else: raise ValueError("unknown model %s"%model_name) @@ -1092,6 +1104,9 @@ def _update(val, axis=None): limit = [0, 0.5, 5, 5][dims] jitter = [0 if v < limit else v for v in jitter] axes.cla() + axes.set_xlabel('x axis') + axes.set_ylabel('y axis') + axes.set_zlabel('z axis') ## Visualize as person on globe #draw_sphere(axes, radius=0.5) diff --git a/sasmodels/models/two_yukawa.py b/sasmodels/models/two_yukawa.py index f1f22c4d..00d110d9 100644 --- a/sasmodels/models/two_yukawa.py +++ b/sasmodels/models/two_yukawa.py @@ -94,7 +94,7 @@ from numpy import inf # TODO: pep8 says packages and modules should not use camel case -from sasmodels.TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF +from sasmodels.TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk # If you want a customized version of two_yukawa as a plugin (for example, # because you want to use the high precision polynomial root solver from mpmath) @@ -104,9 +104,9 @@ # https://github.com/SasView/sasmodels/tree/master/sasmodels/models/two_yukawa.py # https://github.com/SasView/sasmodels/tree/master/sasmodels/TwoYukawa if 0: + import importlib.util import sys from pathlib import Path - import importlib.util # Remove existing TwoYukawa from sys.modules to force a reload remove = [modname for modname in sys.modules if modname.startswith('TwoYukawa.') or modname == 'TwoYukawa'] @@ -121,7 +121,7 @@ sys.modules['TwoYukawa'] = module # Override sasmodels library symbols with the local symbols. - from TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF + from TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk name = "two_yukawa" title = "User model for two Yukawa structure factor (S(q))"