Changeset 7e0b281 in sasmodels for explore/bccpy.py

Timestamp:

Apr 17, 2017 6:34:52 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

64ca163

Parents:

cb038a2

Message:

adjust paracrystal equations to better match terminology in Matsuoka paper

File:

• explore/bccpy.py (modified) (6 diffs)

explore/bccpy.py

@@ -47 +47 @@
 SLD_SOLVENT = 6.3
 
+# Note: using Matsuoka 1990; this is different from what
+# is in the sasmodels/models code  (see bcc vs bcc_old).
+# The difference is that the sign of phi and theta seem to be
+# negative in the old vs. the new, yielding a pattern that is
+# swapped left to right and top to bottom.
+def sc(qa, qb, qc):
+    return qa, qb, qc
+
+def bcc(qa, qb, qc):
+    a1 = (+qa + qb + qc)/2
+    a2 = (-qa - qb + qc)/2
+    a3 = (-qa + qb - qc)/2
+    return a1, a2, a3
+
+def bcc_old(qa, qb, qc):
+    a1 = (+qa + qb - qc)/2.0
+    a2 = (+qa - qb + qc)/2.0
+    a3 = (-qa + qb + qc)/2.0
+    return a1, a2, a3
+
+def fcc(qa, qb, qc):
+    a1 = ( 0. + qb + qc)/2
+    a2 = (-qa + 0. + qc)/2
+    a3 = (-qa + qb + 0.)/2
+    return a1, a2, a3
+
+def fcc_old(qa, qb, qc):
+    a1 = ( qa + qb + 0.)/2
+    a2 = ( qa + 0. + qc)/2
+    a3 = ( 0. + qb + qc)/2
+    return a1, a2, a3
+
+KERNEL = bcc
+
 def kernel(q, dnn, d_factor, theta, phi):
     """
@@ -56 +90 @@
     qc = q*cos(theta)
 
-    if 0: # sc
-        a1, a2, a3 = qa, qb, qc
-        dcos = dnn
-    if 1: # bcc
-        a1 = +qa - qc + qb
-        a2 = +qa + qc - qb
-        a3 = -qa + qc + qb
-        dcos = dnn/2
-    if 0: # fcc
-        a1 = qb + qa
-        a2 = qa + qc
-        a3 = qb + qc
-        dcos = dnn/2
-
-    arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2)
-    exp_arg = exp(-arg)
-    den = [((exp_arg - 2*cos(dcos*a))*exp_arg + 1.0) for a in (a1, a2, a3)]
-    Sq = -expm1(-2*arg)**3/np.prod(den, axis=0)
-    return Sq
+    a1, a2, a3 = KERNEL(qa, qb, qc)
+
+    # Note: paper says that different directions can have different distortion
+    # factors.  Easy enough to add to the code.  This would definitely break
+    # 8-fold symmetry.
+    arg = -0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2)
+    exp_arg = exp(arg)
+    den = [((exp_arg - 2*cos(dnn*a))*exp_arg + 1.0) for a in (a1, a2, a3)]
+    Zq = -expm1(2*arg)**3/np.prod(den, axis=0)
+    return Zq
 
 
@@ -85 +110 @@
     def integrand(theta, phi):
         evals[0] += 1
-        Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
-        return Sq*sin(theta)
+        Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi)
+        return Zq*sin(theta)
     ans = dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi)
     print("dblquad evals =", evals[0])
@@ -134 +159 @@
     Aw = w[None, :] * w[:, None]
     sin_theta = np.fmax(abs(sin(Atheta)), 1e-6)
-    Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
     print("gauss %d evals ="%n, n**2)
-    return np.sum(Sq*Aw*sin_theta)*8/(4*pi)
+    return np.sum(Zq*Aw*sin_theta)*8/(4*pi)
 
 
@@ -148 +173 @@
     phi = np.linspace(0, pi/2, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
-    Sq *= abs(sin(Atheta))
+    Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq *= abs(sin(Atheta))
     dx, dy = theta[1]-theta[0], phi[1]-phi[0]
-    print("rect", n, np.sum(Sq)*dx*dy*8/(4*pi))
-    print("trapz", n, np.trapz(np.trapz(Sq, dx=dx), dx=dy)*8/(4*pi))
-    print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi))
-    print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi))
+    print("rect", n, np.sum(Zq)*dx*dy*8/(4*pi))
+    print("trapz", n, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*8/(4*pi))
+    print("simpson", n, simps(simps(Zq, dx=dx), dx=dy)*8/(4*pi))
+    print("romb", n, romb(romb(Zq, dx=dx), dx=dy)*8/(4*pi))
     print("gridded %d evals ="%n, n**2)
 
@@ -168 +193 @@
     #phi = np.linspace(0, pi/2, n)
     Atheta, Aphi = np.meshgrid(theta, phi)
-    Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
-    Sq *= abs(sin(Atheta))
-    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6)))
+    Zq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi)
+    Zq *= abs(sin(Atheta))
+    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
     pylab.axis('tight')
-    pylab.title("BCC S(q) for q=%g, dnn=%g d_factor=%g" % (q, dnn, d_factor))
+    pylab.title("%s Z(q) for q=%g, dnn=%g d_factor=%g"
+                % (KERNEL.__name__, q, dnn, d_factor))
     pylab.xlabel("theta (degrees)")
     pylab.ylabel("phi (degrees)")