Changeset 9b90bf8 in sasview for src/sas/sascalc

Timestamp:

Jul 19, 2017 10:04:45 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:

a309667

Parents:

d78b5cb

Message:

Use DCT to calculate IDF

File:

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

src/sas/sascalc/corfunc/transform_thread.py

@@ -37 +37 @@
             # ----- 1D Correlation Function -----
             gamma1 = dct((iqs-background)*qs**2)
-            gamma1 = gamma1 / gamma1.max()
+            Q = gamma1.max()
+            gamma1 /= Q
 
             if self.check_if_cancelled(): return
@@ -52 +53 @@
 
             # ----- 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
+            # dmax = 200.0 # Max real space value to calculate IDF up to
+            # dstep = 0.5 # Evaluate the IDF in steps of dstep along the real axis
+            # 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
+            # x_scale = qmax / 0.6
 
-            xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale)
-            idf = np.zeros(len(xgamma))
+            # xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale)
+            # idf = np.zeros(len(xgamma))
+            idf = dct(-qs**4 * (iqs-background))
+            idf[0] = trapz(-qs**4 * (iqs-background), qs)
+            idf /= Q
 
             # 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
+            # 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:
@@ -99 +103 @@
         transform1 = Data1D(xs, gamma1)
         transform3 = Data1D(xs[xs <= 200], gamma3)
-        idf = Data1D(xgamma, idf)
+        idf = Data1D(xs, idf)
 
         transforms = (transform1, transform3, idf)