Changeset 9b90bf8 in sasview for src/sas/sascalc/corfunc


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

File:
1 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) 
Note: See TracChangeset for help on using the changeset viewer.