source: sasmodels/explore/beta/sasfit_compare.py @ 9f2216f

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 9f2216f was e262dd6, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

futher cleanup; include Yun's hardsphere and sphere models from matlab
  • Property mode set to 100644
File size: 28.8 KB
Line 
1from __future__ import division, print_function
2# Make sasmodels available on the path
3import sys, os
4BETA_DIR = os.path.dirname(os.path.realpath(__file__))
5SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR))
6sys.path.insert(0, SASMODELS_DIR)
7
8from collections import namedtuple
9
10from matplotlib import pyplot as plt
11import numpy as np
12from numpy import pi, sin, cos, sqrt, fabs
13from numpy.polynomial.legendre import leggauss
14from scipy.special import j1 as J1
15from numpy import inf
16from scipy.special import gammaln  # type: ignore
17
18Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta')
19Theory.__new__.__defaults__ = (None,) * len(Theory._fields)
20
21#Used to calculate F(q) for the cylinder, sphere, ellipsoid models
22def sas_sinx_x(x):
23    with np.errstate(all='ignore'):
24        retvalue = sin(x)/x
25    retvalue[x == 0.] = 1.
26    return retvalue
27
28def sas_2J1x_x(x):
29    with np.errstate(all='ignore'):
30        retvalue = 2*J1(x)/x
31    retvalue[x == 0] = 1.
32    return retvalue
33
34def sas_3j1x_x(x):
35    """return 3*j1(x)/x"""
36    retvalue = np.empty_like(x)
37    with np.errstate(all='ignore'):
38        # GSL bessel_j1 taylor expansion
39        index = (x < 0.25)
40        y = x[index]**2
41        c1 = -1.0/10.0
42        c2 =  1.0/280.0
43        c3 = -1.0/15120.0
44        c4 =  1.0/1330560.0
45        c5 = -1.0/172972800.0
46        retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))))
47        index = ~index
48        y = x[index]
49        retvalue[index] = 3*(sin(y) - y*cos(y))/y**3
50    retvalue[x == 0.] = 1.
51    return retvalue
52
53#Used to cross check my models with sasview models
54def build_model(model_name, q, **pars):
55    from sasmodels.core import load_model_info, build_model as build_sasmodel
56    from sasmodels.data import empty_data1D
57    from sasmodels.direct_model import DirectModel
58    model_info = load_model_info(model_name)
59    model = build_sasmodel(model_info, dtype='double!')
60    data = empty_data1D(q)
61    calculator = DirectModel(data, model,cutoff=0)
62    calculator.pars = pars.copy()
63    calculator.pars.setdefault('background', 1e-3)
64    return calculator
65
66#gives the hardsphere structure factor that sasview uses
67def _hardsphere_simple(q, radius_effective, volfraction):
68    CUTOFFHS=0.05
69    if fabs(radius_effective) < 1.E-12:
70        HARDSPH=1.0
71        return HARDSPH
72    X = 1.0/( 1.0 -volfraction)
73    D= X*X
74    A= (1.+2.*volfraction)*D
75    A *=A
76    X=fabs(q*radius_effective*2.0)
77    if X < 5.E-06:
78        HARDSPH=1./A
79        return HARDSPH
80    X2 =X*X
81    B = (1.0 +0.5*volfraction)*D
82    B *= B
83    B *= -6.*volfraction
84    G=0.5*volfraction*A
85    if X < CUTOFFHS:
86        FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2
87        HARDSPH= 1./(1. + volfraction*FF )
88        return HARDSPH
89    X4=X2*X2
90    S, C = sin(X), cos(X)
91    FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X
92    HARDSPH= 1./(1. + 24.*volfraction*FF/X2 )
93    return HARDSPH
94
95def hardsphere_simple(q, radius_effective, volfraction):
96    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q]
97    return np.array(SQ)
98
99#Used in gaussian quadrature for polydispersity
100#returns values and the probability of those values based on gaussian distribution
101N_GAUSS = 35
102NSIGMA_GAUSS = 3
103def gaussian_distribution(center, sigma, lb, ub):
104    #3 standard deviations covers approx. 99.7%
105    if sigma != 0:
106        nsigmas = NSIGMA_GAUSS
107        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS)
108        x= x[(x >= lb) & (x <= ub)]
109        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
110        return x, px
111    else:
112        return np.array([center]), np.array([1])
113
114N_SCHULZ = 80
115NSIGMA_SCHULZ = 8
116def schulz_distribution(center, sigma, lb, ub):
117    if sigma != 0:
118        nsigmas = NSIGMA_SCHULZ
119        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ)
120        x= x[(x >= lb) & (x <= ub)]
121        R = x/center
122        z = (center/sigma)**2
123        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
124        px = np.exp(arg)
125        return x, px
126    else:
127        return np.array([center]), np.array([1])
128
129#returns the effective radius used in sasview
130def ER_ellipsoid(radius_polar, radius_equatorial):
131    ee = np.empty_like(radius_polar)
132    if radius_polar > radius_equatorial:
133        ee = (radius_polar**2 - radius_equatorial**2)/radius_polar**2
134    elif radius_polar < radius_equatorial:
135        ee = (radius_equatorial**2 - radius_polar**2) / radius_equatorial**2
136    else:
137        ee = 2*radius_polar
138    if (radius_polar * radius_equatorial != 0):
139        bd = 1.0 - ee
140        e1 = np.sqrt(ee)
141        b1 = 1.0 + np.arcsin(e1) / (e1*np.sqrt(bd))
142        bL = (1.0 + e1) / (1.0 - e1)
143        b2 = 1.0 + bd / 2 / e1 * np.log(bL)
144        delta = 0.75 * b1 * b2
145    ddd = np.zeros_like(radius_polar)
146    ddd = 2.0*(delta + 1.0)*radius_polar*radius_equatorial**2
147    return 0.5*ddd**(1.0 / 3.0)
148
149def ellipsoid_volume(radius_polar,radius_equatorial):
150    volume = (4./3.)*pi*radius_polar*radius_equatorial**2
151    return volume
152
153# F1 is F(q)
154# F2 is F(g)^2
155#IQM is I(q) with monodispersity
156#IQSM is I(q) with structure factor S(q) and monodispersity
157#IQBM is I(q) with Beta Approximation and monodispersity
158#SQ is monodisperse approach for structure factor
159#SQ_EFF is the effective structure factor from beta approx
160def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent,
161                    volfraction=0, radius_effective=None):
162    #creates values z and corresponding probabilities w from legendre-gauss quadrature
163    volume = ellipsoid_volume(radius_polar, radius_equatorial)
164    z, w = leggauss(76)
165    F1 = np.zeros_like(q)
166    F2 = np.zeros_like(q)
167    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
168    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
169    #legendre-gauss quadrature
170    for k, qk in enumerate(q):
171        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2)
172        form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r)
173        F2[k] = np.sum(w*form**2)
174        F1[k] = np.sum(w*form)
175    #the 1/2 comes from the change of variables mentioned above
176    F2 = F2/2.0
177    F1 = F1/2.0
178    if radius_effective is None:
179        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial)
180    SQ = hardsphere_simple(q, radius_effective, volfraction)
181    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
182    IQM = 1e-4*F2/volume
183    IQSM = volfraction*IQM*SQ
184    IQBM = volfraction*IQM*SQ_EFF
185    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
186
187#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
188#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
189def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent,
190                 radius_polar_pd=0.1, radius_equatorial_pd=0.1,
191                 radius_polar_pd_type='gaussian',
192                 radius_equatorial_pd_type='gaussian',
193                 volfraction=0, radius_effective=None,
194                 background=0, scale=1,
195                 norm='sasview'):
196    if norm not in ['sasview', 'sasfit', 'yun']:
197        raise TypeError("unknown norm "+norm)
198    if radius_polar_pd_type == 'gaussian':
199        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
200    elif radius_polar_pd_type == 'schulz':
201        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
202    if radius_equatorial_pd_type == 'gaussian':
203        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
204    elif radius_equatorial_pd_type == 'schulz':
205        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
206    total_weight = total_volume = 0
207    radius_eff = 0
208    F1, F2 = np.zeros_like(q), np.zeros_like(q)
209    for k, Rpk in enumerate(Rp_val):
210        for i, Rei in enumerate(Re_val):
211            theory = ellipsoid_theta(q,Rpk,Rei,sld,sld_solvent)
212            volume = ellipsoid_volume(Rpk, Rei)
213            weight = Rp_prob[k]*Re_prob[i]
214            total_weight += weight
215            total_volume += weight*volume
216            F1 += theory.F1*weight
217            F2 += theory.F2*weight
218            radius_eff += weight*ER_ellipsoid(Rpk,Rei)
219    F1 /= total_weight
220    F2 /= total_weight
221    average_volume = total_volume/total_weight
222    if radius_effective is None:
223        radius_effective = radius_eff/total_weight
224    if norm == 'sasfit':
225        IQD = F2
226    elif norm == 'sasview':
227        # Note: internally, sasview uses F2/total_volume because:
228        #   average_volume = total_volume/total_weight
229        #   F2/total_weight / average_volume
230        #     = F2/total_weight / total_volume/total_weight
231        #     = F2/total_volume
232        IQD = F2/average_volume*1e-4*volfraction
233    elif norm == 'yun':
234        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
235        F2 *= 1e-12
236        IQD = F2/average_volume*1e8*volfraction
237    SQ = hardsphere_simple(q, radius_effective, volfraction)
238    beta = F1**2/F2
239    SQ_EFF = 1 + beta*(SQ - 1)
240    IQSD = IQD*SQ
241    IQBD = IQD*SQ_EFF
242    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
243
244#polydispersity for sphere
245def sphere_r(q,radius,sld,sld_solvent,
246             radius_pd=0.1, radius_pd_type='gaussian',
247             volfraction=0, radius_effective=None,
248             background=0, scale=1,
249             norm='sasview'):
250    if norm not in ['sasview', 'sasfit', 'yun']:
251        raise TypeError("unknown norm "+norm)
252    if radius_pd_type == 'gaussian':
253        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
254    elif radius_pd_type == 'schulz':
255        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
256    total_weight = total_volume = 0
257    F1 = np.zeros_like(q)
258    F2 = np.zeros_like(q)
259    for k, rk in enumerate(radius_val):
260        volume = 4./3.*pi*rk**3
261        total_weight += radius_prob[k]
262        total_volume += radius_prob[k]*volume
263        form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
264        F2 += radius_prob[k]*form**2
265        F1 += radius_prob[k]*form
266    F1 /= total_weight
267    F2 /= total_weight
268    average_volume = total_volume/total_weight
269
270    if radius_effective is None:
271        radius_effective = radius
272    average_volume = total_volume/total_weight
273    if norm == 'sasfit':
274        IQD = F2
275    elif norm == 'sasview':
276        IQD = F2/average_volume*1e-4*volfraction
277    elif norm == 'yun':
278        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
279        F2 *= 1e-12
280        IQD = F2/average_volume*1e8*volfraction
281    SQ = hardsphere_simple(q, radius_effective, volfraction)
282    beta = F1**2/F2
283    SQ_EFF = 1 + beta*(SQ - 1)
284    IQSD = IQD*SQ
285    IQBD = IQD*SQ_EFF
286    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
287
288###############################################################################
289###############################################################################
290###############################################################################
291##################                                           ##################
292##################                   TESTS                   ##################
293##################                                           ##################
294###############################################################################
295###############################################################################
296###############################################################################
297
298def popn(d, keys):
299    """
300    Splits a dict into two, with any key of *d* which is in *keys* removed
301    from *d* and added to *b*. Returns *b*.
302    """
303    b = {}
304    for k in keys:
305        try:
306            b[k] = d.pop(k)
307        except KeyError:
308            pass
309    return b
310
311def sasmodels_theory(q, Pname, **pars):
312    """
313    Call sasmodels to compute the model with and without a hard sphere
314    structure factor.
315    """
316    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()}
317    Ppars = pars.copy()
318    Spars = popn(Ppars, ['radius_effective', 'volfraction'])
319    Ipars = pars.copy()
320
321    # Autofill npts and nsigmas for the given pd type
322    for k, v in pars.items():
323        if k.endswith("_pd_type"):
324            if v == "gaussian":
325                n, nsigmas = N_GAUSS, NSIGMA_GAUSS
326            elif v == "schulz":
327                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ
328            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n)
329            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
330            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n)
331            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
332
333    #Ppars['scale'] = Spars.get('volfraction', 1)
334    P = build_model(Pname, q)
335    S = build_model("hardsphere", q)
336    I = build_model(Pname+"@hardsphere", q)
337    Pq = P(**Ppars)*pars.get('volfraction', 1)
338    #Sq = S(**Spars)
339    Iq = I(**Ipars)
340    #Iq = Pq*Sq*pars.get('volfraction', 1)
341    Sq = Iq/Pq
342    return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None)
343
344def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'):
345    """
346    Plot fields in common between target and actual, along with relative error.
347    """
348    available = [s for s in fields.split()
349                 if getattr(target, s) is not None and getattr(actual, s) is not None]
350    rows = len(available)
351    for row, field in enumerate(available):
352        Q = target.Q
353        I1, I2 = getattr(target, field), getattr(actual, field)
354        plt.subplot(rows, 2, 2*row+1)
355        plt.loglog(Q, abs(I1), label="target "+field)
356        plt.loglog(Q, abs(I2), label="value "+field)
357        #plt.legend(loc="upper left", bbox_to_anchor=(1,1))
358        plt.legend(loc='lower left')
359        plt.subplot(rows, 2, 2*row+2)
360        plt.semilogx(Q, I2/I1 - 1, label="relative error")
361        #plt.semilogx(Q, I1/I2 - 1, label="relative error")
362    plt.tight_layout()
363    plt.suptitle(title)
364    plt.show()
365
366def data_file(name):
367    return os.path.join(BETA_DIR, 'data_files', name)
368
369def load_sasfit(path):
370    data = np.loadtxt(path, dtype=str, delimiter=';').T
371    data = np.vstack((map(float, v) for v in data[0:2]))
372    return data
373
374COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None]
375
376def compare_sasview_sphere(pd_type='schulz'):
377    q = np.logspace(-5, 0, 250)
378    model = 'sphere'
379    pars = dict(
380        radius=20,sld=4,sld_solvent=1,
381        background=0,
382        radius_pd=.1, radius_pd_type=pd_type,
383        volfraction=0.15,
384        #radius_effective=12.59921049894873,  # equivalent average sphere radius
385        )
386    target = sasmodels_theory(q, model, **pars)
387    actual = sphere_r(q, norm='sasview', **pars)
388    title = " ".join(("sasmodels", model, pd_type))
389    compare(title, target, actual)
390COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
391COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
392
393def compare_sasview_ellipsoid(pd_type='gaussian'):
394    q = np.logspace(-5, 0, 50)
395    model = 'ellipsoid'
396    pars = dict(
397        radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,
398        background=0,
399        radius_polar_pd=.1, radius_polar_pd_type=pd_type,
400        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type,
401        volfraction=0.15,
402        #radius_effective=12.59921049894873,
403        )
404    target = sasmodels_theory(q, model, **pars)
405    actual = ellipsoid_pe(q, norm='sasview', **pars)
406    title = " ".join(("sasmodels", model, pd_type))
407    compare(title, target, actual)
408COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')
409COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')
410
411def compare_yun_ellipsoid_mono():
412    pars = {
413        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
414        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
415        'sld': 2, 'sld_solvent': 1,
416        'volfraction': 0.15,
417        # Yun uses radius for same volume sphere for effective radius
418        # whereas sasview uses the average curvature.
419        'radius_effective': 12.59921049894873,
420    }
421
422    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
423    Q = data[0]
424    F1 = data[1]
425    P = data[3]*pars['volfraction']
426    S = data[5]
427    Seff = data[6]
428    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
429    actual = ellipsoid_pe(Q, norm='yun', **pars)
430    title = " ".join(("yun", "ellipsoid", "no pd"))
431    #compare(title, target, actual, fields="P S I Seff Ibeta")
432    compare(title, target, actual)
433COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono
434COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono
435
436def compare_yun_sphere_gauss():
437    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100
438    # With pd = 0.1, that's 14 sigma and 1400 points.
439    pars = {
440        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
441        'sld': 6, 'sld_solvent': 0,
442        'volfraction': 0.1,
443    }
444
445    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'),skiprows=2).T
446    Q = data[0]
447    F1 = data[1]
448    P = data[3]
449    S = data[5]
450    Seff = data[6]
451    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
452    actual = sphere_r(Q, norm='yun', **pars)
453    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
454    compare(title, target, actual)
455    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'),skiprows=2).T
456    pars.update(radius_pd=0.15)
457    Q = data[0]
458    F1 = data[1]
459    P = data[3]
460    S = data[5]
461    Seff = data[6]
462    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
463    actual = sphere_r(Q, norm='yun', **pars)
464    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
465    compare(title, target, actual)
466COMPARISON[('yun','sphere','gaussian')] = compare_yun_sphere_gauss
467
468
469def compare_sasfit_sphere_gauss():
470    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
471    pars = {
472        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
473        'sld': 4, 'sld_solvent': 1,
474        'volfraction': 0.3,
475    }
476
477    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
478    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
479    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
480    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
481    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
482    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
483    actual = sphere_r(Q, norm="sasfit", **pars)
484    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
485    compare(title, target, actual)
486    #compare(title, target, actual, fields="P")
487COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss
488
489def compare_sasfit_sphere_schulz():
490    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
491    #We have scaled the output from sasfit by 1e-4*volume*volfraction
492    #0.10050378152592121
493    pars = {
494        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
495        'sld': 4, 'sld_solvent': 1,
496        'volfraction': 0.3,
497    }
498
499    Q, IQ = load_sasfit(data_file('richard_test.txt'))
500    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
501    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
502    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
503    actual = sphere_r(Q, norm="sasfit", **pars)
504    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
505    compare(title, target, actual)
506COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz
507
508def compare_sasfit_ellipsoid_schulz():
509    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
510    #Effective radius =13.1353356684
511    #We have scaled the output from sasfit by 1e-4*volume*volfraction
512    #0.10050378152592121
513    pars = {
514        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
515        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
516        'sld': 4, 'sld_solvent': 1,
517        'volfraction': 0.3, 'radius_effective': 13.1353356684,
518    }
519
520    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
521    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
522    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
523    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
524    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
525    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
526    compare(title, target, actual)
527COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz
528
529
530def compare_sasfit_ellipsoid_gaussian():
531    pars = {
532        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
533        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
534        'sld': 4, 'sld_solvent': 1,
535        'volfraction': 0, 'radius_effective': None,
536    }
537
538    #Rp=20,Re=10,eta_core=4,eta_solv=1
539    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
540    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
541    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
542    target = Theory(Q=Q, P=PQ0)
543    compare("sasfit ellipsoid no poly", target, actual); plt.show()
544
545    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
546    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
547    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
548    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
549    target = Theory(Q=Q, P=PQ_Rp10)
550    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
551    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
552    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
553    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
554    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
555    target = Theory(Q=Q, P=PQ_Re10)
556    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
557    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
558    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
559    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
560    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
561    target = Theory(Q=Q, P=PQ_Rp30)
562    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
563    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
564    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
565    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
566    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
567    target = Theory(Q=Q, P=PQ_Re30)
568    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
569    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
570    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
571    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
572    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
573    target = Theory(Q=Q, P=PQ_Rp60)
574    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
575    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
576    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
577    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
578    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
579    target = Theory(Q=Q, P=PQ_Re60)
580    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
581
582    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
583    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
584    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
585    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
586    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
587    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
588    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
589    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
590    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
591    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
592    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
593    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
594    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
595    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
596    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
597    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
598    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
599    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
600    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
601    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
602    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
603
604    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
605    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
606    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
607    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
608    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
609    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
610    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
611    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
612    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
613    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
614    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
615    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
616    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
617    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
618    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
619    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
620    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
621    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
622    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
623    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
624    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
625
626    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
627    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
628    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
629    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
630    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
631    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
632    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
633    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
634    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
635    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
636    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
637    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
638    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
639    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
640    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
641    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
642    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
643    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
644    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
645    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
646    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
647COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian
648
649def main():
650    key = tuple(sys.argv[1:])
651    if key not in COMPARISON:
652        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]")
653        return
654    comparison = COMPARISON[key]
655    comparison()
656
657if __name__ == "__main__":
658    main()
Note: See TracBrowser for help on using the repository browser.