source: sasmodels/explore/beta/sasfit_compare.py @ aa44a6a

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

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