source: sasmodels/explore/beta/sasfit_compare.py @ 2cefd79

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

restructure beta approx tests to allow command line options

  • Property mode set to 100644
File size: 28.1 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        F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(qk*r))**2
173        F2[k] = np.sum(w*F2i)
174        F1i = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r)
175        F1[k] = np.sum(w*F1i)
176    #the 1/2 comes from the change of variables mentioned above
177    F2 = F2/2.0
178    F1 = F1/2.0
179    if radius_effective is None:
180        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial)
181    SQ = hardsphere_simple(q, radius_effective, volfraction)
182    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
183    IQM = 1e-4*F2/volume
184    IQSM = volfraction*IQM*SQ
185    IQBM = volfraction*IQM*SQ_EFF
186    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
187
188#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
189#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
190def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent,
191                 radius_polar_pd=0.1, radius_equatorial_pd=0.1,
192                 radius_polar_pd_type='gaussian',
193                 radius_equatorial_pd_type='gaussian',
194                 volfraction=0, radius_effective=None,
195                 background=0, scale=1,
196                 norm='sasview'):
197    if norm not in ['sasview', 'sasfit', 'yun']:
198        raise TypeError("unknown norm "+norm)
199    if radius_polar_pd_type == 'gaussian':
200        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
201    elif radius_polar_pd_type == 'schulz':
202        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
203    if radius_equatorial_pd_type == 'gaussian':
204        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
205    elif radius_equatorial_pd_type == 'schulz':
206        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
207    Normalization = 0
208    F1,F2 = np.zeros_like(q), np.zeros_like(q)
209    radius_eff = total_weight = 0
210    for k, Rpk in enumerate(Rp_val):
211        for i, Rei in enumerate(Re_val):
212            theory = ellipsoid_theta(q,Rpk,Rei,sld,sld_solvent)
213            volume = ellipsoid_volume(Rpk, Rei)
214            if norm == 'sasfit':
215                Normalization += Rp_prob[k]*Re_prob[i]
216            elif norm == 'sasview' or norm == 'yun':
217                Normalization += Rp_prob[k]*Re_prob[i]*volume
218            F1 += theory.F1*Rp_prob[k]*Re_prob[i]
219            F2 += theory.F2*Rp_prob[k]*Re_prob[i]
220            radius_eff += Rp_prob[k]*Re_prob[i]*ER_ellipsoid(Rpk,Rei)
221            total_weight += Rp_prob[k]*Re_prob[i]
222    F1 = F1/Normalization
223    F2 = F2/Normalization
224    if radius_effective is None:
225        radius_effective = radius_eff/total_weight
226    SQ = hardsphere_simple(q, radius_effective, volfraction)
227    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
228    volume = ellipsoid_volume(radius_polar, radius_equatorial)
229    if norm == 'sasfit':
230        IQD = F2
231        IQSD = IQD*SQ
232        IQBD = IQD*SQ_EFF
233    elif norm == 'sasview':
234        IQD = F2*1e-4*volfraction
235        IQSD = IQD*SQ
236        IQBD = IQD*SQ_EFF
237    elif norm == 'yun':
238        SQ_EFF = 1 + Normalization*F1**2/F2*(SQ - 1)
239        F2 = F2/volume
240        IQD = F2
241        IQSD = IQD*SQ
242        IQBD = IQD*SQ_EFF
243    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
244
245#polydispersity for sphere
246def sphere_r(q,radius,sld,sld_solvent,
247             radius_pd=0.1, radius_pd_type='gaussian',
248             volfraction=0, radius_effective=None,
249             background=0, scale=1,
250             norm='sasview'):
251    if norm not in ['sasview', 'sasfit', 'yun']:
252        raise TypeError("unknown norm "+norm)
253    if radius_pd_type == 'gaussian':
254        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
255    elif radius_pd_type == 'schulz':
256        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
257    Normalization=0
258    F1 = np.zeros_like(q)
259    F2 = np.zeros_like(q)
260    for k, rk in enumerate(radius_val):
261        volume = 4./3.*pi*rk**3
262        if norm == 'sasfit':
263            Normalization += radius_prob[k]
264        elif norm == 'sasview' or norm == 'yun':
265            Normalization += radius_prob[k]*volume
266        F2k = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2
267        F1k = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
268        F2 += radius_prob[k]*F2k
269        F1 += radius_prob[k]*F1k
270
271    F2 = F2/Normalization
272    F1 = F1/Normalization
273    if radius_effective is None:
274        radius_effective = radius
275    SQ = hardsphere_simple(q, radius_effective, volfraction)
276    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
277    volume = 4./3.*pi*radius**3
278    if norm == 'sasfit':
279        IQD = F2
280        IQSD = IQD*SQ
281        IQBD = IQD*SQ_EFF
282    elif norm == 'sasview':
283        IQD = F2*1e-4*volfraction
284        IQSD = IQD*SQ
285        IQBD = IQD*SQ_EFF
286    elif norm == 'yun':
287        SQ_EFF = 1 + Normalization*F1**2/F2*(SQ - 1)
288        F2 = F2/volume
289        IQD = F2
290        IQSD = IQD*SQ
291        IQBD = IQD*SQ_EFF
292    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
293
294###############################################################################
295###############################################################################
296###############################################################################
297##################                                           ##################
298##################                   TESTS                   ##################
299##################                                           ##################
300###############################################################################
301###############################################################################
302###############################################################################
303
304def popn(d, keys):
305    """
306    Splits a dict into two, with any key of *d* which is in *keys* removed
307    from *d* and added to *b*. Returns *b*.
308    """
309    b = {}
310    for k in keys:
311        try:
312            b[k] = d.pop(k)
313        except KeyError:
314            pass
315    return b
316
317def sasmodels_theory(q, Pname, **pars):
318    """
319    Call sasmodels to compute the model with and without a hard sphere
320    structure factor.
321    """
322    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()}
323    Ppars = pars.copy()
324    Spars = popn(Ppars, ['radius_effective', 'volfraction'])
325    Ipars = pars.copy()
326
327    # Autofill npts and nsigmas for the given pd type
328    for k, v in pars.items():
329        if k.endswith("_pd_type"):
330            if v == "gaussian":
331                n, nsigmas = N_GAUSS, NSIGMA_GAUSS
332            elif v == "schulz":
333                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ
334            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n)
335            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
336            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n)
337            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
338
339    #Ppars['scale'] = Spars.get('volfraction', 1)
340    P = build_model(Pname, q)
341    S = build_model("hardsphere", q)
342    I = build_model(Pname+"@hardsphere", q)
343    Pq = P(**Ppars)*pars.get('volfraction', 1)
344    #Sq = S(**Spars)
345    Iq = I(**Ipars)
346    #Iq = Pq*Sq*pars.get('volfraction', 1)
347    Sq = Iq/Pq
348    return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None)
349
350def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'):
351    """
352    Plot fields in common between target and actual, along with relative error.
353    """
354    available = [s for s in fields.split()
355                 if getattr(target, s) is not None and getattr(actual, s) is not None]
356    rows = len(available)
357    for row, field in enumerate(available):
358        Q = target.Q
359        I1, I2 = getattr(target, field), getattr(actual, field)
360        plt.subplot(rows, 2, 2*row+1)
361        plt.loglog(Q, abs(I1), label="target "+field)
362        plt.loglog(Q, abs(I2), label="value "+field)
363        #plt.legend(loc="upper left", bbox_to_anchor=(1,1))
364        plt.legend(loc='lower left')
365        plt.subplot(rows, 2, 2*row+2)
366        #plt.semilogx(Q, I2/I1 - 1, label="relative error")
367        plt.semilogx(Q, I1/I2 - 1, label="relative error")
368    plt.tight_layout()
369    plt.suptitle(title)
370    plt.show()
371
372def data_file(name):
373    return os.path.join(BETA_DIR, 'data_files', name)
374
375def load_sasfit(path):
376    data = np.loadtxt(path, dtype=str, delimiter=';').T
377    data = np.vstack((map(float, v) for v in data[0:2]))
378    return data
379
380COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None]
381
382def compare_sasview_sphere(pd_type='schulz'):
383    q = np.logspace(-5, 0, 250)
384    model = 'sphere'
385    pars = dict(
386        radius=20,sld=4,sld_solvent=1,
387        background=0,
388        radius_pd=.1, radius_pd_type=pd_type,
389        volfraction=0.15,
390        #radius_effective=12.59921049894873,  # equivalent average sphere radius
391        )
392    target = sasmodels_theory(q, model, **pars)
393    actual = sphere_r(q, norm='sasview', **pars)
394    title = " ".join(("sasmodels", model, pd_type))
395    compare(title, target, actual)
396COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
397COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
398
399def compare_sasview_ellipsoid(pd_type='gaussian'):
400    q = np.logspace(-5, 0, 50)
401    model = 'ellipsoid'
402    pars = dict(
403        radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,
404        background=0,
405        radius_polar_pd=.1, radius_polar_pd_type=pd_type,
406        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type,
407        volfraction=0.15,
408        #radius_effective=12.59921049894873,
409        )
410    target = sasmodels_theory(q, model, **pars)
411    actual = ellipsoid_pe(q, norm='sasview', **pars)
412    title = " ".join(("sasmodels", model, pd_type))
413    compare(title, target, actual)
414COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
415COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
416
417def compare_yun_ellipsoid_mono():
418    pars = {
419        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
420        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
421        'sld': 2, 'sld_solvent': 1,
422        'volfraction': 0.15,
423        # Yun uses radius for same volume sphere for effective radius
424        # whereas sasview uses the average curvature.
425        'radius_effective': 12.59921049894873,
426    }
427    volume = ellipsoid_volume(pars['radius_polar'], pars['radius_equatorial'])
428
429    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
430    Q = data[0]
431    F1 = data[1]
432    F2 = data[3]
433    S = data[5]
434    Seff = data[6]
435    P = F2
436    I = P*S
437    Ibeta = P*Seff
438    P = I = Ibeta = None
439    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, I=I, Seff=Seff, Ibeta=Ibeta)
440    actual = ellipsoid_pe(Q, norm='yun', **pars)
441    title = " ".join(("yun", "ellipsoid", "no pd"))
442    #compare(title, target, actual, fields="P S I Seff Ibeta")
443    compare(title, target, actual)
444COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono
445COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono
446
447def compare_sasfit_sphere_gauss():
448    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
449    pars = {
450        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
451        'sld': 4, 'sld_solvent': 1,
452        'volfraction': 0.3,
453    }
454    volume = 4./3.*pi*pars['radius']**3
455    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
456    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
457    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
458    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
459    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
460    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
461    actual = sphere_r(Q, norm="sasfit", **pars)
462    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
463    compare(title, target, actual)
464    #compare(title, target, actual, fields="P")
465COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss
466
467def compare_sasfit_sphere_schulz():
468    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
469    #We have scaled the output from sasfit by 1e-4*volume*volfraction
470    #0.10050378152592121
471    pars = {
472        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
473        'sld': 4, 'sld_solvent': 1,
474        'volfraction': 0.3,
475    }
476    volume = 4./3.*pi*pars['radius']**3
477
478    Q, IQ = load_sasfit(data_file('richard_test.txt'))
479    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
480    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
481    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
482    actual = sphere_r(Q, norm="sasfit", **pars)
483    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
484    compare(title, target, actual)
485COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz
486
487def compare_sasfit_ellipsoid_schulz():
488    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
489    #Effective radius =13.1353356684
490    #We have scaled the output from sasfit by 1e-4*volume*volfraction
491    #0.10050378152592121
492    pars = {
493        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
494        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
495        'sld': 4, 'sld_solvent': 1,
496        'volfraction': 0.3, 'radius_effective': 13.1353356684,
497    }
498    volume = ellipsoid_volume(pars['radius_polar'], pars['radius_equatorial'])
499    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
500    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
501    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
502    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
503    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
504    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
505    compare(title, target, actual)
506COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz
507
508
509def compare_sasfit_ellipsoid_gaussian():
510    pars = {
511        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
512        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
513        'sld': 4, 'sld_solvent': 1,
514        'volfraction': 0, 'radius_effective': None,
515    }
516
517    #Rp=20,Re=10,eta_core=4,eta_solv=1
518    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
519    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
520    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
521    target = Theory(Q=Q, P=PQ0)
522    compare("sasfit ellipsoid no poly", target, actual); plt.show()
523
524    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
525    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
526    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
527    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
528    target = Theory(Q=Q, P=PQ_Rp10)
529    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
530    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
531    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
532    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
533    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
534    target = Theory(Q=Q, P=PQ_Re10)
535    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
536    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
537    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
538    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
539    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
540    target = Theory(Q=Q, P=PQ_Rp30)
541    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
542    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
543    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
544    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
545    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
546    target = Theory(Q=Q, P=PQ_Re30)
547    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
548    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
549    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
550    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
551    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
552    target = Theory(Q=Q, P=PQ_Rp60)
553    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
554    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
555    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
556    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
557    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
558    target = Theory(Q=Q, P=PQ_Re60)
559    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
560
561    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
562    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
563    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
564    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
565    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
566    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
567    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
568    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
569    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
570    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
571    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
572    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
573    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
574    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
575    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
576    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
577    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
578    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
579    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
580    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
581    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
582
583    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
584    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
585    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
586    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
587    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
588    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
589    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
590    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
591    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
592    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
593    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
594    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
595    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
596    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
597    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
598    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
599    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
600    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
601    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
602    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
603    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
604
605    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
606    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
607    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
608    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
609    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
610    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
611    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
612    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
613    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
614    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
615    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
616    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
617    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
618    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
619    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
620    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
621    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
622    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
623    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
624    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
625    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
626COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_schulz
627
628def main():
629    key = tuple(sys.argv[1:])
630    if key not in COMPARISON:
631        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]")
632        return
633    comparison = COMPARISON[key]
634    comparison()
635
636if __name__ == "__main__":
637    main()
Note: See TracBrowser for help on using the repository browser.