Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • src/sas/sascalc/corfunc/transform_thread.py

    r626f833 r37d53225  
    1515        self.extrap_fn = extrap_fn 
    1616 
     17    def check_if_cancelled(self): 
     18        if self.isquit(): 
     19            self.update("Fourier transform cancelled.") 
     20            self.complete(transforms=None) 
     21            return True 
     22        return False 
     23 
    1724    def compute(self): 
    1825        qs = self.extrapolation.x 
     
    2734        self.ready(delay=0.0) 
    2835 
    29         if self.isquit(): 
    30             return 
     36        if self.check_if_cancelled(): return 
    3137        try: 
     38            # ----- 1D Correlation Function ----- 
    3239            gamma1 = dct((iqs-background)*qs**2) 
    3340            gamma1 = gamma1 / gamma1.max() 
    3441 
     42            if self.check_if_cancelled(): return 
     43 
     44            # ----- 3D Correlation Function ----- 
    3545            # gamma3(R) = 1/R int_{0}^{R} gamma1(x) dx 
    3646            # trapz uses the trapezium rule to calculate the integral 
     
    3949            gamma3.insert(0, 1.0) # Gamma_3(0) is defined as 1 
    4050            gamma3 = np.array(gamma3) 
     51 
     52            if self.check_if_cancelled(): return 
     53 
     54            # ----- Interface Distribution function ----- 
     55            dmax = 200.0 # Max real space value to calculate IDF up to 
     56            dstep = 0.5 
     57            qmax = 1.0 # Max q value to integrate up to when calculating IDF 
     58 
     59            # Units of x axis depend on qmax (for some reason?). This scales 
     60            # the xgamma array appropriately, since qmax was set to 0.6 in 
     61            # the original fortran code. 
     62            x_scale = qmax / 0.6 
     63 
     64            xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale) 
     65            idf = np.zeros(len(xgamma)) 
     66 
     67            # nth moment = integral(q^n * I(q), q=0, q=inf) 
     68            moment = np.zeros(5) 
     69            for n in range(5): 
     70                integrand = qs**n * (iqs-background) 
     71                moment[n] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
     72                if self.check_if_cancelled(): return 
     73 
     74            # idf(x) = -integral(q^4 * I(q)*cos(qx), q=0, q=inf) / 2nd moment 
     75            # => idf(0) = -integral(q^4 * I(q), 0, inf) / (2nd moment) 
     76            #  = -(4th moment)/(2nd moment) 
     77            idf[0] = -moment[4] / moment[2] 
     78            for i in range(1, len(xgamma)): 
     79                d = xgamma[i] 
     80 
     81                integrand = -qs**4 * (iqs-background) * np.cos(d*qs) 
     82                idf[i] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
     83                idf[i] /= moment[2] 
     84                if self.check_if_cancelled(): return 
     85 
     86            xgamma *= x_scale 
     87 
    4188        except Exception as e: 
    4289            import logging 
     
    53100        transform1 = Data1D(xs, gamma1) 
    54101        transform3 = Data1D(xs[xs <= 200], gamma3) 
     102        idf = Data1D(xgamma, idf) 
    55103 
    56         transforms = (transform1, transform3) 
     104        transforms = (transform1, transform3, idf) 
    57105 
    58106        self.complete(transforms=transforms) 
Note: See TracChangeset for help on using the changeset viewer.