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

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

Update Yun's code to produce <F2> and P columns

  • 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    F2 = data[2]
449    P = data[3]
450    S = data[5]
451    Seff = data[6]
452    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff)
453    actual = sphere_r(Q, norm='yun', **pars)
454    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
455    compare(title, target, actual)
456    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'),skiprows=2).T
457    pars.update(radius_pd=0.15)
458    Q = data[0]
459    F1 = data[1]
460    F2 = data[2]
461    P = data[3]
462    S = data[5]
463    Seff = data[6]
464    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff)
465    actual = sphere_r(Q, norm='yun', **pars)
466    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
467    compare(title, target, actual)
468COMPARISON[('yun','sphere','gaussian')] = compare_yun_sphere_gauss
469
470
471def compare_sasfit_sphere_gauss():
472    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
473    pars = {
474        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
475        'sld': 4, 'sld_solvent': 1,
476        'volfraction': 0.3,
477    }
478
479    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
480    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
481    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
482    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
483    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
484    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
485    actual = sphere_r(Q, norm="sasfit", **pars)
486    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
487    compare(title, target, actual)
488    #compare(title, target, actual, fields="P")
489COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss
490
491def compare_sasfit_sphere_schulz():
492    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
493    #We have scaled the output from sasfit by 1e-4*volume*volfraction
494    #0.10050378152592121
495    pars = {
496        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
497        'sld': 4, 'sld_solvent': 1,
498        'volfraction': 0.3,
499    }
500
501    Q, IQ = load_sasfit(data_file('richard_test.txt'))
502    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
503    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
504    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
505    actual = sphere_r(Q, norm="sasfit", **pars)
506    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
507    compare(title, target, actual)
508COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz
509
510def compare_sasfit_ellipsoid_schulz():
511    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
512    #Effective radius =13.1353356684
513    #We have scaled the output from sasfit by 1e-4*volume*volfraction
514    #0.10050378152592121
515    pars = {
516        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
517        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
518        'sld': 4, 'sld_solvent': 1,
519        'volfraction': 0.3, 'radius_effective': 13.1353356684,
520    }
521
522    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
523    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
524    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
525    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
526    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
527    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
528    compare(title, target, actual)
529COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz
530
531
532def compare_sasfit_ellipsoid_gaussian():
533    pars = {
534        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
535        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
536        'sld': 4, 'sld_solvent': 1,
537        'volfraction': 0, 'radius_effective': None,
538    }
539
540    #Rp=20,Re=10,eta_core=4,eta_solv=1
541    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
542    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
543    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
544    target = Theory(Q=Q, P=PQ0)
545    compare("sasfit ellipsoid no poly", target, actual); plt.show()
546
547    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
548    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
549    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
550    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
551    target = Theory(Q=Q, P=PQ_Rp10)
552    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
553    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
554    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
555    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
556    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
557    target = Theory(Q=Q, P=PQ_Re10)
558    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
559    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
560    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
561    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
562    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
563    target = Theory(Q=Q, P=PQ_Rp30)
564    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
565    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
566    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
567    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
568    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
569    target = Theory(Q=Q, P=PQ_Re30)
570    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
571    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
572    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
573    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
574    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
575    target = Theory(Q=Q, P=PQ_Rp60)
576    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
577    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
578    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
579    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
580    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
581    target = Theory(Q=Q, P=PQ_Re60)
582    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
583
584    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
585    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
586    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
587    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
588    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
589    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
590    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
591    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
592    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
593    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
594    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
595    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
596    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
597    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
598    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
599    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
600    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
601    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
602    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
603    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
604    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
605
606    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
607    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
608    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
609    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
610    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
611    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
612    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
613    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
614    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
615    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
616    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
617    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
618    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
619    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
620    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
621    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
622    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
623    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
624    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
625    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
626    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
627
628    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
629    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
630    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
631    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
632    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
633    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
634    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
635    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
636    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
637    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
638    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
639    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
640    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
641    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
642    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
643    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
644    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
645    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
646    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
647    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
648    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
649COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian
650
651def main():
652    key = tuple(sys.argv[1:])
653    if key not in COMPARISON:
654        print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]")
655        return
656    comparison = COMPARISON[key]
657    comparison()
658
659if __name__ == "__main__":
660    main()
Note: See TracBrowser for help on using the repository browser.