Changeset ac995be in sasmodels for sasmodels/sesans.py

Timestamp:

Jun 12, 2018 10:29:02 AM (6 years ago)

Author:

Adam Washington <adam.washington@…>

Branches:

master

Children:

b032200

Parents:

d7af1c6

Message:

Get log sesans roughly working

File:

• sasmodels/sesans.py (modified) (5 diffs)

sasmodels/sesans.py

@@ -14 +14 @@
 import numpy as np  # type: ignore
 from numpy import pi  # type: ignore
-from scipy.special import j0, j1
+from scipy.special import j1
 
 
@@ -39 +39 @@
 
     # transform arrays
-    _H = None  # type: np.ndarray
-    _H0 = None # type: np.ndarray
+    _H = None   # type: np.ndarray
+    _H0 = None  # type: np.ndarray
 
     def __init__(self, z, SElength, lam, zaccept, Rmax):
         # type: (np.ndarray, float, float) -> None
-        #import logging; logging.info("creating SESANS transform")
         self.q = z
         self._set_hankel(SElength, lam, zaccept, Rmax)
@@ -52 +51 @@
         G0 = np.dot(self._H0, Iq)
         G = np.dot(self._H.T, Iq)
+        G0 = G[0]  # FIXME: This is a kludge
         P = G - G0
         return P
@@ -57 +57 @@
     def _set_hankel(self, SElength, lam, zaccept, Rmax):
         # type: (np.ndarray, float, float) -> None
-        # Force float32 arrays, otherwise run into memory problems on some machines
+        # Force float32 arrays, otherwise run into memory problems on
+        # some machines
         SElength = np.asarray(SElength, dtype='float32')
 
@@ -64 +65 @@
         q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
         # q = np.arange(q_min, q_max, q_min, dtype='float32')
-        q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2),
-                             dtype=np.float32))
-        q = np.hstack([0], q)
+        # q = np.exp(np.arange(np.log(q_min), np.log(q_max), np.log(2),
+        #                      dtype=np.float32))
+        q = np.exp(np.linspace(np.log(q_min), np.log(q_max), 10*SElength.size,
+                               dtype=np.float32))
+        q = np.hstack([[0], q])
 
-        H0 = np.pi * (q[1:]**2 - q[-1]**2)
+        H0 = np.pi * (q[1:]**2 - q[:-1]**2) * (q[1:] - q[:-1])
 
         # repq = np.tile(q, (SElength.size, 1)).T
         H = np.outer(q, SElength)
         j1(H, out=H)
-        H *= q
+        H *= q.reshape((-1, 1))
         H = H[1:] - H[:-1]
         H *= 2 * np.pi / SElength
 
         lam = np.asarray(lam, dtype=np.float32)
-        reptheta = np.outer(q, lam)
+        reptheta = np.outer(q[1:], lam)
         reptheta /= np.float32(2*np.pi)
         np.arcsin(reptheta, out=reptheta)
         # reptheta = np.arcsin(repq*replam/2*np.pi)
         mask = reptheta > zaccept
-        H[mask] = 0
+        # H[mask] = 0
 
         # H = np.zeros((q.size, SElength.size), dtype=np.float32)
         # H0 = q * 0
-        assert(H.shape == (q.size, SElength.size))
+        assert(H.shape == (q.size-1, SElength.size))
 
-        self.q_calc = q
+        self.q_calc = q[1:]
         self._H, self._H0 = H, H0