source: sasmodels/explore/beta/sasfit_compare_new.py @ d299327

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d299327 was 3f818b2, checked in by richardh, 6 years ago

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