Changeset 37d53225 in sasview

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

Location:

src/sas

Files:

• sascalc/corfunc/transform_thread.py (modified) (4 diffs)
• sasgui/perspectives/corfunc/corfunc.py (modified) (1 diff)
• sasgui/perspectives/corfunc/corfunc_panel.py (modified) (2 diffs)
• sasgui/perspectives/corfunc/plot_labels.py (modified) (1 diff)

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)

src/sas/sasgui/perspectives/corfunc/corfunc.py

@@ -189 +189 @@
             # Show the transformation as a curve instead of points
             new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
+        elif label == IDF_LABEL:
+            new_plot.xaxis("{x}", 'A')
+            new_plot.yaxis("{G}", '')
+
+            new_plot.xtransform = 'x'
+            new_plot.ytransform = 'y'
+            group_id = GROUP_ID_IDF
+            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
         new_plot.id = label
         new_plot.name = label

src/sas/sasgui/perspectives/corfunc/corfunc_panel.py

@@ -276 +276 @@
             return
         self._transformed_data = transforms
-        (transform1, transform3) = transforms
+        (transform1, transform3, idf) = transforms
         import numpy as np
         plot_x = transform1.x[np.where(transform1.x <= 200)]
@@ -283 +283 @@
         # No need to shorten gamma3 as it's only calculated up to x=200
         self._manager.show_data(transform3, TRANSFORM_LABEL3)
+        self._manager.show_data(idf, IDF_LABEL)
         # Only enable extract params button if a fourier trans. has been done
         if self.transform_type == 'fourier':

src/sas/sasgui/perspectives/corfunc/plot_labels.py

@@ -6 +6 @@
 TRANSFORM_LABEL1 = r"$\Gamma_1(x)$"
 TRANSFORM_LABEL3 = r"$\Gamma_3(x)$"
+
+GROUP_ID_IDF = r"$G(x)$"
+IDF_LABEL = r"$G(x)$"