source: sasmodels/explore/beta/sasfit_compare_new.py @ 2a64722

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

tweak explore/beta for python 3
  • Property mode set to 100644
File size: 40.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
22# RKH adding vesicle and hollow_cylinder to test sasview special cases of ER and VR
23# BEWARE there are issues here if you run this from python3 (i.e. qt5 version of sasview), so do "source activate sasview"
24def sas_sinx_x(x):
25    with np.errstate(all='ignore'):
26        retvalue = sin(x)/x
27    retvalue[x == 0.] = 1.
28    return retvalue
29
30def sas_2J1x_x(x):
31    with np.errstate(all='ignore'):
32        retvalue = 2*J1(x)/x
33    retvalue[x == 0] = 1.
34    return retvalue
35
36def sas_3j1x_x(x):
37    """return 3*j1(x)/x"""
38    retvalue = np.empty_like(x)
39    with np.errstate(all='ignore'):
40        # GSL bessel_j1 taylor expansion
41        index = (x < 0.25)
42        y = x[index]**2
43        c1 = -1.0/10.0
44        c2 = +1.0/280.0
45        c3 = -1.0/15120.0
46        c4 = +1.0/1330560.0
47        c5 = -1.0/172972800.0
48        retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))))
49        index = ~index
50        y = x[index]
51        retvalue[index] = 3*(sin(y) - y*cos(y))/y**3
52    retvalue[x == 0.] = 1.
53    return retvalue
54
55#Used to cross check my models with sasview models
56def build_model(model_name, q, **pars):
57    from sasmodels.core import load_model_info, build_model as build_sasmodel
58    from sasmodels.data import empty_data1D
59    from sasmodels.direct_model import DirectModel
60    model_info = load_model_info(model_name)
61    model = build_sasmodel(model_info, dtype='double!')
62    data = empty_data1D(q)
63    calculator = DirectModel(data, model,cutoff=0)
64    calculator.pars = pars.copy()
65    calculator.pars.setdefault('background', 0)
66    return calculator
67
68#gives the hardsphere structure factor that sasview uses
69def _hardsphere_simple(q, radius_effective, volfraction):
70    CUTOFFHS = 0.05
71    if fabs(radius_effective) < 1.E-12:
72        HARDSPH = 1.0
73        return HARDSPH
74    X = 1.0/(1.0 -volfraction)
75    D = X*X
76    A = (1.+2.*volfraction)*D
77    A *= A
78    X = fabs(q*radius_effective*2.0)
79    if X < 5.E-06:
80        HARDSPH = 1./A
81        return HARDSPH
82    X2 = X*X
83    B = (1.0 +0.5*volfraction)*D
84    B *= B
85    B *= -6.*volfraction
86    G = 0.5*volfraction*A
87    if X < CUTOFFHS:
88        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
89        HARDSPH = 1./(1. + volfraction*FF )
90        return HARDSPH
91    X4 = X2*X2
92    S, C = sin(X), cos(X)
93    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
94    HARDSPH = 1./(1. + 24.*volfraction*FF/X2)
95    return HARDSPH
96
97def hardsphere_simple(q, radius_effective, volfraction):
98    SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q]
99    return np.array(SQ)
100
101#Used in gaussian quadrature for polydispersity
102#returns values and the probability of those values based on gaussian distribution
103N_GAUSS = 35
104NSIGMA_GAUSS = 3
105def gaussian_distribution(center, sigma, lb, ub):
106    #3 standard deviations covers approx. 99.7%
107    if sigma != 0:
108        nsigmas = NSIGMA_GAUSS
109        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS)
110        x = x[(x >= lb) & (x <= ub)]
111        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
112        return x, px
113    else:
114        return np.array([center]), np.array([1])
115
116N_SCHULZ = 80
117NSIGMA_SCHULZ = 8
118def schulz_distribution(center, sigma, lb, ub):
119    if sigma != 0:
120        nsigmas = NSIGMA_SCHULZ
121        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ)
122        x = x[(x >= lb) & (x <= ub)]
123        R = x/center
124        z = (center/sigma)**2
125        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z)
126        px = np.exp(arg)
127        return x, px
128    else:
129        return np.array([center]), np.array([1])
130
131#returns the effective radius used in sasview
132def ER_ellipsoid(radius_polar, radius_equatorial):
133    ee = np.empty_like(radius_polar)
134    if radius_polar > radius_equatorial:
135        ee = (radius_polar**2 - radius_equatorial**2)/radius_polar**2
136    elif radius_polar < radius_equatorial:
137        ee = (radius_equatorial**2 - radius_polar**2) / radius_equatorial**2
138    else:
139        ee = 2*radius_polar
140    if radius_polar * radius_equatorial != 0:
141        bd = 1.0 - ee
142        e1 = np.sqrt(ee)
143        b1 = 1.0 + np.arcsin(e1) / (e1*np.sqrt(bd))
144        bL = (1.0 + e1) / (1.0 - e1)
145        b2 = 1.0 + bd / 2 / e1 * np.log(bL)
146        delta = 0.75 * b1 * b2
147    ddd = np.zeros_like(radius_polar)
148    ddd = 2.0*(delta + 1.0)*radius_polar*radius_equatorial**2
149    return 0.5*ddd**(1.0 / 3.0)
150
151def ellipsoid_volume(radius_polar, radius_equatorial):
152    volume = (4./3.)*pi*radius_polar*radius_equatorial**2
153    return volume
154
155# F1 is F(q)
156# F2 is F(g)^2
157#IQM is I(q) with monodispersity
158#IQSM is I(q) with structure factor S(q) and monodispersity
159#IQBM is I(q) with Beta Approximation and monodispersity
160#SQ is monodisperse approach for structure factor
161#SQ_EFF is the effective structure factor from beta approx
162def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent,
163                    volfraction=0, radius_effective=None):
164    #creates values z and corresponding probabilities w from legendre-gauss quadrature
165    volume = ellipsoid_volume(radius_polar, radius_equatorial)
166    z, w = leggauss(76)
167    F1 = np.zeros_like(q)
168    F2 = np.zeros_like(q)
169    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
170    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
171    #legendre-gauss quadrature
172    for k, qk in enumerate(q):
173        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2)
174        form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r)
175        F2[k] = np.sum(w*form**2)
176        F1[k] = np.sum(w*form)
177    #the 1/2 comes from the change of variables mentioned above
178    F2 = F2/2.0
179    F1 = F1/2.0
180    if radius_effective is None:
181        radius_effective = ER_ellipsoid(radius_polar,radius_equatorial)
182    SQ = hardsphere_simple(q, radius_effective, volfraction)
183    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
184    IQM = 1e-4*F2/volume
185    IQSM = volfraction*IQM*SQ
186    IQBM = volfraction*IQM*SQ_EFF
187    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
188
189#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
190#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
191def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent,
192                 radius_polar_pd=0.1, radius_equatorial_pd=0.1,
193                 radius_polar_pd_type='gaussian',
194                 radius_equatorial_pd_type='gaussian',
195                 volfraction=0, radius_effective=None,
196                 background=0, scale=1,
197                 norm='sasview'):
198    if norm not in ['sasview', 'sasfit', 'yun']:
199        raise TypeError("unknown norm "+norm)
200    if radius_polar_pd_type == 'gaussian':
201        Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
202    elif radius_polar_pd_type == 'schulz':
203        Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf)
204    if radius_equatorial_pd_type == 'gaussian':
205        Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
206    elif radius_equatorial_pd_type == 'schulz':
207        Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf)
208    total_weight = total_volume = 0
209    radius_eff = 0
210    F1, F2 = np.zeros_like(q), np.zeros_like(q)
211    for k, Rpk in enumerate(Rp_val):
212        print("ellipsoid cycle", k, "of", len(Rp_val))
213        for i, Rei in enumerate(Re_val):
214            theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent)
215            volume = ellipsoid_volume(Rpk, Rei)
216            weight = Rp_prob[k]*Re_prob[i]
217            total_weight += weight
218            total_volume += weight*volume
219            F1 += theory.F1*weight
220            F2 += theory.F2*weight
221            radius_eff += weight*ER_ellipsoid(Rpk, Rei)
222    F1 /= total_weight
223    F2 /= total_weight
224    average_volume = total_volume/total_weight
225    if radius_effective is None:
226        radius_effective = radius_eff/total_weight
227    if norm == 'sasfit':
228        IQD = F2
229    elif norm == 'sasview':
230        # Note: internally, sasview uses F2/total_volume because:
231        #   average_volume = total_volume/total_weight
232        #   F2/total_weight / average_volume
233        #     = F2/total_weight / total_volume/total_weight
234        #     = F2/total_volume
235        IQD = F2/average_volume*1e-4*volfraction
236        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
237        F2 *= 1e-4
238    elif norm == 'yun':
239        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
240        F2 *= 1e-12
241        IQD = F2/average_volume*1e8*volfraction
242    SQ = hardsphere_simple(q, radius_effective, volfraction)
243    beta = F1**2/F2
244    SQ_EFF = 1 + beta*(SQ - 1)
245    IQSD = IQD*SQ
246    IQBD = IQD*SQ_EFF
247    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
248
249#polydispersity for sphere
250def sphere_r(q,radius,sld,sld_solvent,
251             radius_pd=0.1, radius_pd_type='gaussian',
252             volfraction=0, radius_effective=None,
253             background=0, scale=1,
254             norm='sasview'):
255    if norm not in ['sasview', 'sasfit', 'yun']:
256        raise TypeError("unknown norm "+norm)
257    if radius_pd_type == 'gaussian':
258        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
259    elif radius_pd_type == 'schulz':
260        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
261    total_weight = total_volume = 0
262    F1 = np.zeros_like(q)
263    F2 = np.zeros_like(q)
264    for k, rk in enumerate(radius_val):
265        volume = 4./3.*pi*rk**3
266        total_weight += radius_prob[k]
267        total_volume += radius_prob[k]*volume
268        form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
269        F2 += radius_prob[k]*form**2
270        F1 += radius_prob[k]*form
271    F1 /= total_weight
272    F2 /= total_weight
273    average_volume = total_volume/total_weight
274
275    if radius_effective is None:
276        radius_effective = radius
277        print("radius_effective  ",radius_effective)
278        average_volume = total_volume/total_weight
279    if norm == 'sasfit':
280        IQD = F2
281    elif norm == 'sasview':
282        IQD = F2/average_volume*1e-4*volfraction
283    elif norm == 'yun':
284        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
285        F2 *= 1e-12
286        IQD = F2/average_volume*1e8*volfraction
287    SQ = hardsphere_simple(q, radius_effective, volfraction)
288    beta = F1**2/F2
289    SQ_EFF = 1 + beta*(SQ - 1)
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#polydispersity for vesicle
295def vesicle_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
296        radius_pd=0.1, thickness_pd=0.2, radius_pd_type="gaussian", thickness_pd_type="gaussian",
297        radius_effective=None, background=0, scale=1,
298                 norm='sasview'):
299    if norm not in ['sasview', 'sasfit', 'yun']:
300        raise TypeError("unknown norm "+norm)
301    if radius_pd_type == 'gaussian':
302        R_val, R_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
303    elif radius_pd_type == 'schulz':
304        R_val, R_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
305    if thickness_pd_type == 'gaussian':
306        T_val, T_prob = gaussian_distribution(thickness, thickness_pd*thickness, 0, inf)
307    elif thickness_pd_type == 'schulz':
308        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf)
309    total_weight = total_volume = total_shell = 0
310    radius_eff = 0
311    F1, F2 = np.zeros_like(q), np.zeros_like(q)
312    for k, Rk in enumerate(R_val):
313        print("vesicle cycle", k, "of", len(R_val)," Rk ",Rk)
314        volume_in = 4./3.*pi*Rk**3
315        form_in = (sld-sld_solvent)*volume_in*sas_3j1x_x(q*Rk)
316        for i, Ti in enumerate(T_val):
317            Rout = Rk + Ti
318            volume_out = 4./3.*pi*Rout**3
319            form_out = (sld-sld_solvent)*volume_out*sas_3j1x_x(q*Rout)
320            form = form_out - form_in
321            volume = volume_out - volume_in
322            weight = R_prob[k]*T_prob[i]
323            total_weight += weight
324            total_shell += weight*volume
325            total_volume += weight*volume_out
326            F1 += form*weight
327            F2 += form**2*weight
328            radius_eff += weight*Rout
329    F1 /= total_weight
330    F2 /= total_weight
331    average_volume = total_volume/total_weight
332    if radius_effective is None:
333        radius_effective = radius_eff/total_weight
334        print("radius_effective  ",radius_effective)
335    if norm == 'sasfit':
336        IQD = F2
337    elif norm == 'sasview':
338        # Note: internally, sasview uses F2/total_volume because:
339        #   average_volume = total_volume/total_weight
340        #   F2/total_weight / average_volume
341        #     = F2/total_weight / total_volume/total_weight
342        #     = F2/total_volume
343        IQD = F2/average_volume*1e-4*volfraction
344        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
345        F2 *= 1e-4
346    elif norm == 'yun':
347        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
348        F2 *= 1e-12
349        IQD = F2/average_volume*1e8*volfraction
350    # RKH SORT THIS OUT WHEN HAVE NEW MODEL
351    vfff=volfraction*(30./20.)**3
352    print("radius_effective, fudged volfaction for S(Q) ",radius_effective,vfff)
353    SQ = hardsphere_simple(q, radius_effective, vfff)
354    beta = F1**2/F2
355    SQ_EFF = 1 + beta*(SQ - 1)
356    IQSD = IQD*SQ
357    IQBD = IQD*SQ_EFF
358    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
359#
360#polydispersity for hollow_cylinder
361def hollow_cylinder_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.15,
362        radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", thickness_pd_type="gaussian",
363        radius_effective=None, background=0, scale=1,
364                 norm='sasview'):
365    if norm not in ['sasview', 'sasfit', 'yun']:
366        raise TypeError("unknown norm "+norm)
367    if radius_pd_type == 'gaussian':
368        R_val, R_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
369    elif radius_pd_type == 'schulz':
370        R_val, R_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
371    if thickness_pd_type == 'gaussian':
372        T_val, T_prob = gaussian_distribution(thickness, thickness_pd*thickness, 0, inf)
373    elif thickness_pd_type == 'schulz':
374        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf)
375    if length_pd_type == 'gaussian':
376        T_val, T_prob = gaussian_distribution(length, length_pd*length, 0, inf)
377    elif length_pd_type == 'schulz':
378        L_val, L_prob = schulz_distribution(length, length_pd*length, 0, inf)
379    total_weight = total_volume = total_shell = 0
380    radius_eff = 0
381    F1, F2 = np.zeros_like(q), np.zeros_like(q)
382    for k, Rk in enumerate(R_val):
383        print("hollow_cylinder cycle", k, "of", len(R_val)," Rk ",Rk)
384        for i, Ti in enumerate(T_val):
385            for j, Lj in enumerate(L_val):
386                Rout = Rk + Ti
387                volume_out = pi*Rout**2*Lj
388                volume_in = pi*Rk**2*Lj
389                volume = volume_out - volume_in
390                theory = hollow_cylinder_theta(q, Rk, Ti, Lj, sld, sld_solvent,
391                    volfraction=0, radius_effective=None)
392                weight = R_prob[k]*T_prob[i]*L_prob[j]
393                total_weight += weight
394                total_shell += weight*volume
395                total_volume += weight*volume_out
396                F1 += theory.F1*weight
397                F2 += theory.F2*weight
398                radius_eff += weight*ER_hollow_cylinder(radius, thickness, length)
399    F1 /= total_weight
400    F2 /= total_weight
401    average_volume = total_volume/total_weight
402    if radius_effective is None:
403        radius_effective = radius_eff/total_weight
404        print("radius_effective  ",radius_effective)
405    if norm == 'sasfit':
406        IQD = F2
407    elif norm == 'sasview':
408        # Note: internally, sasview uses F2/total_volume because:
409        #   average_volume = total_volume/total_weight
410        #   F2/total_weight / average_volume
411        #     = F2/total_weight / total_volume/total_weight
412        #     = F2/total_volume
413        IQD = F2/average_volume*1e-4*volfraction
414        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2
415        F2 *= 1e-4
416    elif norm == 'yun':
417        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2
418        F2 *= 1e-12
419        IQD = F2/average_volume*1e8*volfraction
420# RKH SORT THIS OUT WHEN HAVE NEW MODEL
421    vfff=volfraction
422    print("radius_effective, volfaction for S(Q) ",radius_effective,vfff)
423    SQ = hardsphere_simple(q, radius_effective, vfff)
424    beta = F1**2/F2
425    SQ_EFF = 1 + beta*(SQ - 1)
426    IQSD = IQD*SQ
427    IQBD = IQD*SQ_EFF
428    return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
429
430# RKH copying Paul K's ellipsoid here, this computes everything for monodisperse case
431# and supplies partial results for polydisperse case.
432def hollow_cylinder_theta(q, radius, thickness, length, sld, sld_solvent,
433                    volfraction=0, radius_effective=None):
434    #creates values z and corresponding probabilities w from legendre-gauss quadrature
435    z, w = leggauss(76)
436    F1 = np.zeros_like(q)
437    F2 = np.zeros_like(q)
438    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
439    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
440    #legendre-gauss quadrature
441    lower = 0.0
442    upper = 1.0
443    gamma_sq = (radius/(radius+thickness))**2
444    for k, qk in enumerate(q):
445        for i, zt in enumerate(z):
446            # quadrature loop
447            cos_theta = 0.5*( zt* (upper-lower) + lower + upper  )
448            sin_theta = sqrt(1.0 - cos_theta*cos_theta)
449            aaa=(radius+thickness)*qk*sin_theta
450            ## sas_J1 expects an array, so try direct use of J1
451            lam1 = J1(aaa)/aaa
452            aaa=radius*qk*sin_theta
453            lam2 = J1(aaa)/aaa
454            psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq)
455            #Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
456            #Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
457            aaa=0.5*length*qk*cos_theta
458            t2 = sin(aaa)/aaa
459#            t2 = sas_sinx_x(0.5*length*qk*cos_theta)
460            form = psi*t2
461            F1[i] += w[i] * form
462            F2[i] += w[i] * form * form
463    volume = pi*((radius + thickness)**2 - radius**2)*length
464    s = (sld - sld_solvent) * volume
465    F2 = F2*s*(upper - lower)/2.0
466    F1 = F1*s*s*(upper - lower)/2.0
467    if radius_effective is None:
468        radius_effective = ER_hollow_cylinder(radius, thickness, length)
469    SQ = hardsphere_simple(q, radius_effective, volfraction)
470    SQ_EFF = 1 + F1**2/F2*(SQ - 1)
471    IQM = 1e-4*F2/volume
472    IQSM = volfraction*IQM*SQ
473    IQBM = volfraction*IQM*SQ_EFF
474    return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM)
475#
476def ER_hollow_cylinder(radius, thickness, length):
477    """
478    :param radius:      Cylinder core radius
479    :param thickness:   Cylinder wall thickness
480    :param length:      Cylinder length
481    :return:            Effective radius
482    """
483    router = radius + thickness
484    if router == 0 or length == 0:
485        return 0.0
486    len1 = router
487    len2 = length/2.0
488    term1 = len1*len1*2.0*len2/2.0
489    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0)
490    ddd = 3.0*term1*term2
491    diam = pow(ddd, (1.0/3.0))
492    return diam
493
494###############################################################################
495###############################################################################
496###############################################################################
497##################                                           ##################
498##################                   TESTS                   ##################
499##################                                           ##################
500###############################################################################
501###############################################################################
502###############################################################################
503
504def popn(d, keys):
505    """
506    Splits a dict into two, with any key of *d* which is in *keys* removed
507    from *d* and added to *b*. Returns *b*.
508    """
509    b = {}
510    for k in keys:
511        try:
512            b[k] = d.pop(k)
513        except KeyError:
514            pass
515    return b
516
517def sasmodels_theory(q, Pname, **pars):
518    """
519    Call sasmodels to compute the model with and without a hard sphere
520    structure factor.
521    """
522    #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()}
523    Ppars = pars.copy()
524    Spars = popn(Ppars, ['radius_effective', 'volfraction'])
525    Ipars = pars.copy()
526
527    # Autofill npts and nsigmas for the given pd type
528    for k, v in pars.items():
529        if k.endswith("_pd_type"):
530            if v == "gaussian":
531                n, nsigmas = N_GAUSS, NSIGMA_GAUSS
532            elif v == "schulz":
533                n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ
534            Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n)
535            Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
536            Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n)
537            Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas)
538
539    #Ppars['scale'] = Spars.get('volfraction', 1)
540    P = build_model(Pname, q)
541    S = build_model("hardsphere", q)
542    I = build_model(Pname+"@hardsphere", q)
543    Pq = P(**Ppars)*pars.get('volfraction', 1)
544    Sq = S(**Spars)
545    Iq = I(**Ipars)
546    #Iq = Pq*Sq*pars.get('volfraction', 1)
547    #Sq = Iq/Pq
548    #Iq = None#= Sq = None
549    r = dict(I._kernel.results())
550    #print(r)
551    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq)
552
553def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'):
554    """
555    Plot fields in common between target and actual, along with relative error.
556    """
557    available = [s for s in fields.split()
558                 if getattr(target, s) is not None and getattr(actual, s) is not None]
559    rows = len(available)
560    for row, field in enumerate(available):
561        Q = target.Q
562        I1, I2 = getattr(target, field), getattr(actual, field)
563        plt.subplot(rows, 2, 2*row+1)
564        plt.loglog(Q, abs(I1), label="target "+field)
565        plt.loglog(Q, abs(I2), label="value "+field)
566        #plt.legend(loc="upper left", bbox_to_anchor=(1,1))
567        plt.legend(loc='lower left')
568        plt.subplot(rows, 2, 2*row+2)
569        plt.semilogx(Q, I2/I1 - 1, label="relative error")
570        #plt.semilogx(Q, I1/I2 - 1, label="relative error")
571    plt.tight_layout()
572    plt.suptitle(title)
573    plt.show()
574
575def data_file(name):
576    return os.path.join(BETA_DIR, 'data_files', name)
577
578def load_sasfit(path):
579    data = np.loadtxt(path, dtype=str, delimiter=';').T
580    data = np.vstack((list(map(float, v)) for v in data[0:2]))
581    return data
582
583COMPARISON = {}  # Type: Dict[(str,str,str)] -> Callable[(), None]
584
585def compare_sasview_sphere(pd_type='schulz'):
586    q = np.logspace(-5, 0, 250)
587    model = 'sphere'
588    pars = dict(
589        radius=20, sld=4, sld_solvent=1,
590        background=0,
591        radius_pd=.1, radius_pd_type=pd_type,
592        volfraction=0.15,
593        radius_effective=12.59921049894873,  # equivalent average sphere radius
594        )
595    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
596    actual = sphere_r(q, norm='sasview', **pars)
597    title = " ".join(("sasmodels", model, pd_type))
598    compare(title, target, actual)
599COMPARISON[('sasview', 'sphere', 'gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
600COMPARISON[('sasview', 'sphere', 'schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
601
602
603def compare_sasview_vesicle(pd_type='gaussian'):
604    q = np.logspace(-5, 0, 250)
605    model = 'vesicle'
606    print("F2 seems OK, but big issues with S(Q), wo work in progress")
607# NOTE parameers for vesicle model are soon to be changed to remove "volfraction"
608    pars = dict(
609        radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
610        background=0,
611        radius_pd=0.1, thickness_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type
612        #,radius_effective=12.59921049894873,  # equivalent average sphere radius
613        )
614    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
615    actual = vesicle_pe(q, norm='sasview', **pars)
616    title = " ".join(("sasmodels", model, pd_type))
617    compare(title, target, actual)
618COMPARISON[('sasview', 'vesicle', 'gaussian')] = lambda: compare_sasview_vesicle(pd_type='gaussian')
619COMPARISON[('sasview', 'vesicle', 'schulz')] = lambda: compare_sasview_vesicle(pd_type='schulz')
620
621def compare_sasview_hollow_cylinder(pd_type='gaussian'):
622    q = np.logspace(-5, 0, 250)
623    model = 'hollow_cylinder'
624    print(model)
625    print("just about works for monodisperse, polydisperse needs work")
626# NOTE parameters for hollow_cylinder model are soon to be changed to remove "volfraction"
627# setting all three pd to zero is OK
628    pars = dict(
629        radius=20, thickness=10, length=80, sld=4, sld_solvent=1,
630        background=0,
631        radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type
632        #,radius_effective=12.59921049894873,  # equivalent average sphere radius
633        )
634    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
635    actual = hollow_cylinder_pe(q, norm='sasview', **pars)
636# RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 )
637    title = " ".join(("sasmodels", model, pd_type))
638    compare(title, target, actual)
639COMPARISON[('sasview', 'hollow_cylinder', 'gaussian')] = lambda: compare_sasview_hollow_cylinder(pd_type='gaussian')
640COMPARISON[('sasview', 'hollow_cylinder', 'schulz')] = lambda: compare_sasview_hollow_cylinder(pd_type='schulz')
641
642
643def compare_sasview_ellipsoid(pd_type='gaussian'):
644    q = np.logspace(-5, 0, 50)
645    model = 'ellipsoid'
646    pars = dict(
647        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1,
648        background=0,
649        radius_polar_pd=0.1, radius_polar_pd_type=pd_type,
650        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type,
651        volfraction=0.15,
652        radius_effective=270.7543927018,
653        #radius_effective=12.59921049894873,
654        )
655    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars)
656    actual = ellipsoid_pe(q, norm='sasview', **pars)
657# RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.)
658    title = " ".join(("sasmodels", model, pd_type))
659    compare(title, target, actual)
660COMPARISON[('sasview', 'ellipsoid', 'gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')
661COMPARISON[('sasview', 'ellipsoid', 'schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')
662
663def compare_yun_ellipsoid_mono():
664    pars = {
665        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
666        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
667        'sld': 2, 'sld_solvent': 1,
668        'volfraction': 0.15,
669        # Yun uses radius for same volume sphere for effective radius
670        # whereas sasview uses the average curvature.
671        'radius_effective': 12.59921049894873,
672    }
673
674    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
675    Q = data[0]
676    F1 = data[1]
677    P = data[3]*pars['volfraction']
678    S = data[5]
679    Seff = data[6]
680    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
681    actual = ellipsoid_pe(Q, norm='yun', **pars)
682    title = " ".join(("yun", "ellipsoid", "no pd"))
683    #compare(title, target, actual, fields="P S I Seff Ibeta")
684    compare(title, target, actual)
685COMPARISON[('yun', 'ellipsoid', 'gaussian')] = compare_yun_ellipsoid_mono
686COMPARISON[('yun', 'ellipsoid', 'schulz')] = compare_yun_ellipsoid_mono
687
688def compare_yun_sphere_gauss():
689    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100
690    # With pd = 0.1, that's 14 sigma and 1400 points.
691    pars = {
692        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
693        'sld': 6, 'sld_solvent': 0,
694        'volfraction': 0.1,
695    }
696
697    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T
698    Q = data[0]
699    F1 = data[1]
700    F2 = data[2]
701    P = data[3]
702    S = data[5]
703    Seff = data[6]
704    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
705    actual = sphere_r(Q, norm='yun', **pars)
706    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
707    compare(title, target, actual)
708    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T
709    pars.update(radius_pd=0.15)
710    Q = data[0]
711    F1 = data[1]
712    F2 = data[2]
713    P = data[3]
714    S = data[5]
715    Seff = data[6]
716    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
717    actual = sphere_r(Q, norm='yun', **pars)
718    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
719    compare(title, target, actual)
720COMPARISON[('yun', 'sphere', 'gaussian')] = compare_yun_sphere_gauss
721
722
723def compare_sasfit_sphere_gauss():
724    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
725    pars = {
726        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
727        'sld': 4, 'sld_solvent': 1,
728        'volfraction': 0.3,
729    }
730
731    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
732    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
733    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
734    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
735    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
736    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
737    actual = sphere_r(Q, norm="sasfit", **pars)
738    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
739    compare(title, target, actual)
740    #compare(title, target, actual, fields="P")
741COMPARISON[('sasfit', 'sphere', 'gaussian')] = compare_sasfit_sphere_gauss
742
743def compare_sasfit_sphere_schulz():
744    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
745    #We have scaled the output from sasfit by 1e-4*volume*volfraction
746    #0.10050378152592121
747    pars = {
748        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
749        'sld': 4, 'sld_solvent': 1,
750        'volfraction': 0.3,
751    }
752
753    Q, IQ = load_sasfit(data_file('richard_test.txt'))
754    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
755    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
756    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
757    actual = sphere_r(Q, norm="sasfit", **pars)
758    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
759    compare(title, target, actual)
760COMPARISON[('sasfit', 'sphere', 'schulz')] = compare_sasfit_sphere_schulz
761
762def compare_sasfit_ellipsoid_schulz():
763    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
764    #Effective radius =13.1353356684
765    #We have scaled the output from sasfit by 1e-4*volume*volfraction
766    #0.10050378152592121
767    pars = {
768        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
769        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
770        'sld': 4, 'sld_solvent': 1,
771        'volfraction': 0.3, 'radius_effective': 13.1353356684,
772    }
773
774    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
775    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
776    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
777    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
778    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
779    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
780    compare(title, target, actual)
781COMPARISON[('sasfit', 'ellipsoid', 'schulz')] = compare_sasfit_ellipsoid_schulz
782
783
784def compare_sasfit_ellipsoid_gaussian():
785    pars = {
786        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
787        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
788        'sld': 4, 'sld_solvent': 1,
789        'volfraction': 0, 'radius_effective': None,
790    }
791
792    #Rp=20,Re=10,eta_core=4,eta_solv=1
793    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
794    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
795    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
796    target = Theory(Q=Q, P=PQ0)
797    compare("sasfit ellipsoid no poly", target, actual); plt.show()
798
799    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
800    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
801    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
802    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
803    target = Theory(Q=Q, P=PQ_Rp10)
804    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
805    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
806    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
807    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
808    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
809    target = Theory(Q=Q, P=PQ_Re10)
810    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
811    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
812    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
813    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
814    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
815    target = Theory(Q=Q, P=PQ_Rp30)
816    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
817    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
818    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
819    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
820    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
821    target = Theory(Q=Q, P=PQ_Re30)
822    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
823    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
824    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
825    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
826    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
827    target = Theory(Q=Q, P=PQ_Rp60)
828    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
829    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
830    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
831    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
832    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
833    target = Theory(Q=Q, P=PQ_Re60)
834    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
835
836    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
837    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
838    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
839    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
840    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
841    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
842    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
843    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
844    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
845    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
846    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
847    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
848    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
849    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
850    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
851    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
852    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
853    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
854    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
855    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
856    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
857
858    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
859    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
860    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
861    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
862    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
863    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
864    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
865    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
866    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
867    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
868    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
869    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
870    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
871    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
872    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
873    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
874    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
875    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
876    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
877    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
878    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
879
880    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
881    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
882    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
883    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
884    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
885    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
886    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
887    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
888    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
889    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
890    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
891    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
892    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
893    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
894    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
895    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
896    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
897    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
898    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
899    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
900    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
901COMPARISON[('sasfit', 'ellipsoid', 'gaussian')] = compare_sasfit_ellipsoid_gaussian
902
903def main():
904    key = tuple(sys.argv[1:])
905    if key not in COMPARISON:
906        print("Usage: python sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid|vesicle|hollow_cylinder] [gaussian|schulz]")
907        print("But not for [yun sphere schulz], and vesicle & hollow_cylinder only with sasview.")
908        print("Note this compares chosen source against internal python code here, with adjustment for known scale etc issues.")
909        return
910    comparison = COMPARISON[key]
911    comparison()
912
913if __name__ == "__main__":
914    main()
Note: See TracBrowser for help on using the repository browser.