Changeset 37d53225 in sasview for src/sas/sascalc/corfunc

Timestamp:

Jul 18, 2017 10:55:35 AM (7 years ago)

Author:

lewis

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

7dda833

Parents:

626f833

Message:

Implement computation of IDF

File:

• src/sas/sascalc/corfunc/transform_thread.py (modified) (4 diffs)

src/sas/sascalc/corfunc/transform_thread.py

@@ -15 +15 @@
         self.extrap_fn = extrap_fn
 
+    def check_if_cancelled(self):
+        if self.isquit():
+            self.update("Fourier transform cancelled.")
+            self.complete(transforms=None)
+            return True
+        return False
+
     def compute(self):
         qs = self.extrapolation.x
@@ -27 +34 @@
         self.ready(delay=0.0)
 
-        if self.isquit():
-            return
+        if self.check_if_cancelled(): return
         try:
+            # ----- 1D Correlation Function -----
             gamma1 = dct((iqs-background)*qs**2)
             gamma1 = gamma1 / gamma1.max()
 
+            if self.check_if_cancelled(): return
+
+            # ----- 3D Correlation Function -----
             # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx
             # trapz uses the trapezium rule to calculate the integral
@@ -39 +49 @@
             gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1
             gamma3 = np.array(gamma3)
+
+            if self.check_if_cancelled(): return
+
+            # ----- Interface Distribution function -----
+            dmax = 200.0 # Max real space value to calculate IDF up to
+            dstep = 0.5
+            qmax = 1.0 # Max q value to integrate up to when calculating IDF
+
+            # Units of x axis depend on qmax (for some reason?). This scales
+            # the xgamma array appropriately, since qmax was set to 0.6 in
+            # the original fortran code.
+            x_scale = qmax / 0.6
+
+            xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale)
+            idf = np.zeros(len(xgamma))
+
+            # nth moment = integral(q^n * I(q), q=0, q=inf)
+            moment = np.zeros(5)
+            for n in range(5):
+                integrand = qs**n * (iqs-background)
+                moment[n] = trapz(integrand[qs < qmax], qs[qs < qmax])
+                if self.check_if_cancelled(): return
+
+            # idf(x) = -integral(q^4 * I(q)*cos(qx), q=0, q=inf) / 2nd moment
+            # => idf(0) = -integral(q^4 * I(q), 0, inf) / (2nd moment)
+            #  = -(4th moment)/(2nd moment)
+            idf[0] = -moment[4] / moment[2]
+            for i in range(1, len(xgamma)):
+                d = xgamma[i]
+
+                integrand = -qs**4 * (iqs-background) * np.cos(d*qs)
+                idf[i] = trapz(integrand[qs < qmax], qs[qs < qmax])
+                idf[i] /= moment[2]
+                if self.check_if_cancelled(): return
+
+            xgamma *= x_scale
+
         except Exception as e:
             import logging
@@ -53 +100 @@
         transform1 = Data1D(xs, gamma1)
         transform3 = Data1D(xs[xs <= 200], gamma3)
+        idf = Data1D(xgamma, idf)
 
-        transforms = (transform1, transform3)
+        transforms = (transform1, transform3, idf)
 
         self.complete(transforms=transforms)