Changeset 5b5ea20 in sasmodels for explore/jitter.py

Timestamp:

Oct 30, 2017 10:55:42 AM (7 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:

904cd9c

Parents:

87a6591

Message:

[possibly] correct weighting for the azimuthal equidistant projection

File:

• explore/jitter.py (modified) (4 diffs)

explore/jitter.py

@@ -197 +197 @@
             return Rz(longitude)*Ry(latitude)
         def weight(theta_i, phi_j, wi, wj):
+            # Weighting for each point comes from the integral:
+            #     \int\int I(q, lat, log) sin(lat) dlat dlog
+            # We are doing a conformal mapping from disk to sphere, so we need
+            # a change of variables g(theta, phi) -> (lat, long):
+            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta)
+            # giving:
+            #     dtheta dphi = det(J) dlat dlong
+            # where J is the jacobian from the partials of g. Using
+            #     R = sqrt(theta^2 + phi^2),
+            # then
+            #     J = [[x/R, Y/R], -y/R^2, x/R^2]]
+            # and
+            #     det(J) = 1/R
+            # with the final integral being:
+            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi
+            #
+            # This does approximately the right thing, decreasing the weight
+            # of each point as you go farther out on the disk, but it hasn't
+            # yet been checked against the 1D integral results. Prior
+            # to declaring this "good enough" and checking that integrals
+            # work in practice, we will examine alternative mappings.
+            #
+            # The issue is that the mapping does not support the case of free
+            # rotation about a single axis correctly, with a small deviation
+            # in the orthogonal axis independent of the first axis.  Like the
+            # usual polar coordiates integration, the integrated sections
+            # form wedges, though at least in this case the wedge cuts through
+            # the entire sphere, and treats theta and phi identically.
             latitude = sqrt(theta_i**2 + phi_j**2)
-            w = wi*wj if latitude < 180 else 0
-            return w
+            w = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return w*wi*wj if latitude < 180 else 0
     elif PROJECTION == 'azimuthal_equal_area':
         # https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection
@@ -209 +237 @@
             return Rz(longitude)*Ry(latitude)
         def weight(theta_i, phi_j, wi, wj):
-            R = sqrt(theta_i**2 + phi_j**2)
-            return wi*wj if R <= 180 else 0
-            #return wi*wj if latitude <= 180 else 0
+            latitude = sqrt(theta_i**2 + phi_j**2)
+            w = sin(radians(latitude))/latitude if latitude != 0 else 1
+            return w*wi*wj if latitude < 180 else 0
     elif SCALED_PHI == 10:  # random thrashing
         def rotate(theta_i, phi_j):
@@ -234 +262 @@
                   for wi, theta_i in zip(weights, dtheta*t)
                   for wj, phi_j in zip(weights, dphi*t)])
+    #print(max(w), min(w), min(w[w>0]))
     points = points[:, w>0]
     w = w[w>0]
+    w /= max(w)
 
     if 0: # Kent distribution
@@ -504 +534 @@
 
     ## use gaussian distribution unless testing integration
-    #dist = 'rectangle'
-    dist = 'gaussian'
+    dist = 'rectangle'
+    #dist = 'gaussian'
 
     ## initial view