Changeset 9b90bf8 in sasview


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

Location:
src/sas
Files:
3 edited

Legend:

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

    rc728295 r9b90bf8  
    3737            # ----- 1D Correlation Function ----- 
    3838            gamma1 = dct((iqs-background)*qs**2) 
    39             gamma1 = gamma1 / gamma1.max() 
     39            Q = gamma1.max() 
     40            gamma1 /= Q 
    4041 
    4142            if self.check_if_cancelled(): return 
     
    5253 
    5354            # ----- Interface Distribution function ----- 
    54             dmax = 200.0 # Max real space value to calculate IDF up to 
    55             dstep = 0.5 
    56             qmax = 1.0 # Max q value to integrate up to when calculating IDF 
     55            # dmax = 200.0 # Max real space value to calculate IDF up to 
     56            # dstep = 0.5 # Evaluate the IDF in steps of dstep along the real axis 
     57            # qmax = 1.0 # Max q value to integrate up to when calculating IDF 
    5758 
    5859            # Units of x axis depend on qmax (for some reason?). This scales 
    5960            # the xgamma array appropriately, since qmax was set to 0.6 in 
    6061            # the original fortran code. 
    61             x_scale = qmax / 0.6 
     62            # x_scale = qmax / 0.6 
    6263 
    63             xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale) 
    64             idf = np.zeros(len(xgamma)) 
     64            # xgamma = np.arange(0, dmax/x_scale, step=dstep/x_scale) 
     65            # idf = np.zeros(len(xgamma)) 
     66            idf = dct(-qs**4 * (iqs-background)) 
     67            idf[0] = trapz(-qs**4 * (iqs-background), qs) 
     68            idf /= Q 
    6569 
    6670            # nth moment = integral(q^n * I(q), q=0, q=inf) 
    67             moment = np.zeros(5) 
    68             for n in range(5): 
    69                 integrand = qs**n * (iqs-background) 
    70                 moment[n] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
    71                 if self.check_if_cancelled(): return 
    72  
    73             # idf(x) = -integral(q^4 * I(q)*cos(qx), q=0, q=inf) / 2nd moment 
    74             # => idf(0) = -integral(q^4 * I(q), 0, inf) / (2nd moment) 
    75             # = -(4th moment)/(2nd moment) 
    76             idf[0] = -moment[4] / moment[2] 
    77             for i in range(1, len(xgamma)): 
    78                 d = xgamma[i] 
    79  
    80                 integrand = -qs**4 * (iqs-background) * np.cos(d*qs) 
    81                 idf[i] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
    82                 idf[i] /= moment[2] 
    83                 if self.check_if_cancelled(): return 
    84  
    85             xgamma *= x_scale 
     71            # moment = np.zeros(5) 
     72            # for n in range(5): 
     73            #     integrand = qs**n * (iqs-background) 
     74            #     moment[n] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
     75            #     if self.check_if_cancelled(): return 
     76            # 
     77            # # idf(x) = -integral(q^4 * I(q)*cos(qx), q=0, q=inf) / 2nd moment 
     78            # # => idf(0) = -integral(q^4 * I(q), 0, inf) / (2nd moment) 
     79            # # = -(4th moment)/(2nd moment) 
     80            # idf[0] = -moment[4] / moment[2] 
     81            # for i in range(1, len(xgamma)): 
     82            #     d = xgamma[i] 
     83            # 
     84            #     integrand = -qs**4 * (iqs-background) * np.cos(d*qs) 
     85            #     idf[i] = trapz(integrand[qs < qmax], qs[qs < qmax]) 
     86            #     idf[i] /= moment[2] 
     87            #     if self.check_if_cancelled(): return 
     88            # 
     89            # xgamma *= x_scale 
    8690 
    8791        except Exception as e: 
     
    99103        transform1 = Data1D(xs, gamma1) 
    100104        transform3 = Data1D(xs[xs <= 200], gamma3) 
    101         idf = Data1D(xgamma, idf) 
     105        idf = Data1D(xs, idf) 
    102106 
    103107        transforms = (transform1, transform3, idf) 
  • src/sas/sasgui/perspectives/corfunc/corfunc.py

    r7dda833 r9b90bf8  
    192192            new_plot.xaxis("{x}", 'A') 
    193193            new_plot.yaxis("{g_1}", '') 
    194  
     194            # Linear scale 
    195195            new_plot.xtransform = 'x' 
    196196            new_plot.ytransform = 'y' 
    197197            group_id = GROUP_ID_IDF 
     198            # Show IDF as a curve instead of points 
    198199            new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM 
    199200        new_plot.id = label 
  • src/sas/sasgui/perspectives/corfunc/corfunc_panel.py

    rc728295 r9b90bf8  
    277277        self._transformed_data = transforms 
    278278        (transform1, transform3, idf) = transforms 
    279         plot_x = transform1.x[np.where(transform1.x <= 200)] 
    280         plot_y = transform1.y[np.where(transform1.x <= 200)] 
     279        plot_x = transform1.x[transform1.x <= 200] 
     280        plot_y = transform1.y[transform1.x <= 200] 
    281281        self._manager.show_data(Data1D(plot_x, plot_y), TRANSFORM_LABEL1) 
    282282        # No need to shorten gamma3 as it's only calculated up to x=200 
    283283        self._manager.show_data(transform3, TRANSFORM_LABEL3) 
    284         self._manager.show_data(idf, IDF_LABEL) 
     284 
     285        plot_x = idf.x[idf.x <= 200] 
     286        plot_y = idf.y[idf.x <= 200] 
     287        self._manager.show_data(Data1D(plot_x, plot_y), IDF_LABEL) 
    285288 
    286289        # Only enable extract params button if a fourier trans. has been done 
Note: See TracChangeset for help on using the changeset viewer.