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

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since f691fac was 0b8a1fc, checked in by richardh, 6 years ago

progressing towards unit test with S(Q)
  • Property mode set to 100644
File size: 41.4 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# There were issues here from python3 (i.e. qt5 version of sasview),fixed by Paul K (else 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        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    print("in sphere_r : radius_effective  ",radius_effective,"  volfraction",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("in vesicle_pe : 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  Vtot/Vshell = (R+T)^3/ ( (R+T)^3 -R^3 )
351    vfff=volfraction*( 1.0/(1.0 - ( (radius/(radius+thickness))**3 )))
352    print("in vesicle_pe : radius_effective, fudged volfraction 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, length=80, 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", length_pd_type="gaussian", 
363        thickness_pd_type="gaussian", 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        L_val, L_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("in hollow_cylinder : 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("in hollow_cylinder : 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=20  # equivalent average sphere radius
594#       NOTE sasview computes its own radius_effective in "target" (the print(r) at end of sasmodels_theory will reveal its value),
595#       this one is only used locally for "actual"
596        )
597    target = sasmodels_theory(q, model, effective_radius_type=0, structure_factor_mode=1, **pars)
598    actual = sphere_r(q, norm='sasview', **pars)
599    title = " ".join(("sasmodels", model, pd_type))
600    compare(title, target, actual)
601COMPARISON[('sasview', 'sphere', 'gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')
602COMPARISON[('sasview', 'sphere', 'schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')
603
604
605def compare_sasview_vesicle(pd_type='gaussian'):
606    q = np.logspace(-5, 0, 250)
607    model = 'vesicle'
608    print("F2 seems OK, but big issues with S(Q), so work in progress")
609# NOTE parameters for vesicle model are soon to be changed to remove "volfraction" (though still there in 5.0)
610    pars = dict(
611        radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03,
612        background=0,
613        radius_pd=0.01, thickness_pd=0.01, radius_pd_type=pd_type, thickness_pd_type=pd_type,
614        radius_effective=30. ) 
615        # equivalent average sphere radius for local "actual" to match what sasview uses, use None to compute average outer radius here,
616
617    target = sasmodels_theory(q, model, effective_radius_type=0, structure_factor_mode=1, **pars)
618    actual = vesicle_pe(q, norm='sasview', **pars)
619    title = " ".join(("sasmodels", model, pd_type))
620    compare(title, target, actual)
621COMPARISON[('sasview', 'vesicle', 'gaussian')] = lambda: compare_sasview_vesicle(pd_type='gaussian')
622COMPARISON[('sasview', 'vesicle', 'schulz')] = lambda: compare_sasview_vesicle(pd_type='schulz')
623
624def compare_sasview_hollow_cylinder(pd_type='gaussian'):
625    q = np.logspace(-5, 0, 250)
626    model = 'hollow_cylinder'
627    print(model)
628    print("just about works for monodisperse, polydisperse needs work")
629# NOTE parameters for hollow_cylinder model are soon to be changed to remove "volfraction"
630# setting all three pd to zero is OK
631    pars = dict(
632        radius=20, thickness=10, length=80, sld=4, sld_solvent=1,
633        background=0,
634        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,
635        radius_effective=40.687) 
636        # equivalent average sphere radius for local "actual" to match what sasview uses
637    target = sasmodels_theory(q, model, effective_radius_type=0, structure_factor_mode=1, **pars)
638    actual = hollow_cylinder_pe(q, norm='sasview', **pars)
639# RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 )
640    title = " ".join(("sasmodels", model, pd_type))
641    compare(title, target, actual)
642COMPARISON[('sasview', 'hollow_cylinder', 'gaussian')] = lambda: compare_sasview_hollow_cylinder(pd_type='gaussian')
643COMPARISON[('sasview', 'hollow_cylinder', 'schulz')] = lambda: compare_sasview_hollow_cylinder(pd_type='schulz')
644
645
646def compare_sasview_ellipsoid(pd_type='gaussian'):
647    q = np.logspace(-5, 0, 50)
648    model = 'ellipsoid'
649    pars = dict(
650        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1,
651        background=0,
652        radius_polar_pd=0.1, radius_polar_pd_type=pd_type,
653        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type,
654        volfraction=0.15,
655        radius_effective=270.7543927018,
656# if change radius_effective to some other value, the S(Q) from sasview does not agree
657        )
658    target = sasmodels_theory(q, model, effective_radius_type=0, structure_factor_mode=1, **pars)
659    actual = ellipsoid_pe(q, norm='sasview', **pars)
660# RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.)
661    title = " ".join(("sasmodels", model, pd_type))
662    compare(title, target, actual)
663COMPARISON[('sasview', 'ellipsoid', 'gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')
664COMPARISON[('sasview', 'ellipsoid', 'schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')
665
666def compare_yun_ellipsoid_mono():
667    pars = {
668        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
669        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
670        'sld': 2, 'sld_solvent': 1,
671        'volfraction': 0.15,
672        # Yun uses radius for same volume sphere for effective radius
673        # whereas sasview uses the average curvature.
674        'radius_effective': 12.59921049894873,
675    }
676
677    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T
678    Q = data[0]
679    F1 = data[1]
680    P = data[3]*pars['volfraction']
681    S = data[5]
682    Seff = data[6]
683    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
684    actual = ellipsoid_pe(Q, norm='yun', **pars)
685    title = " ".join(("yun", "ellipsoid", "no pd"))
686    #compare(title, target, actual, fields="P S I Seff Ibeta")
687    compare(title, target, actual)
688COMPARISON[('yun', 'ellipsoid', 'gaussian')] = compare_yun_ellipsoid_mono
689COMPARISON[('yun', 'ellipsoid', 'schulz')] = compare_yun_ellipsoid_mono
690
691def compare_yun_sphere_gauss():
692    # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100
693    # With pd = 0.1, that's 14 sigma and 1400 points.
694    pars = {
695        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
696        'sld': 6, 'sld_solvent': 0,
697        'volfraction': 0.1,
698    }
699
700    data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T
701    Q = data[0]
702    F1 = data[1]
703    F2 = data[2]
704    P = data[3]
705    S = data[5]
706    Seff = data[6]
707    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
708    actual = sphere_r(Q, norm='yun', **pars)
709    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf"))
710    compare(title, target, actual)
711    data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T
712    pars.update(radius_pd=0.15)
713    Q = data[0]
714    F1 = data[1]
715    F2 = data[2]
716    P = data[3]
717    S = data[5]
718    Seff = data[6]
719    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff)
720    actual = sphere_r(Q, norm='yun', **pars)
721    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf"))
722    compare(title, target, actual)
723COMPARISON[('yun', 'sphere', 'gaussian')] = compare_yun_sphere_gauss
724
725
726def compare_sasfit_sphere_gauss():
727    #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
728    pars = {
729        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian',
730        'sld': 4, 'sld_solvent': 1,
731        'volfraction': 0.3,
732    }
733
734    Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt'))
735    Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt'))
736    Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt'))
737    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt'))
738    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt'))
739    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD)
740    actual = sphere_r(Q, norm="sasfit", **pars)
741    title = " ".join(("sasfit", "sphere", "pd=10% gaussian"))
742    compare(title, target, actual)
743    #compare(title, target, actual, fields="P")
744COMPARISON[('sasfit', 'sphere', 'gaussian')] = compare_sasfit_sphere_gauss
745
746def compare_sasfit_sphere_schulz():
747    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
748    #We have scaled the output from sasfit by 1e-4*volume*volfraction
749    #0.10050378152592121
750    pars = {
751        'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz',
752        'sld': 4, 'sld_solvent': 1,
753        'volfraction': 0.3,
754    }
755
756    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt'))
757    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt'))
758    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt'))
759    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
760    actual = sphere_r(Q, norm="sasfit", **pars)
761    title = " ".join(("sasfit", "sphere", "pd=10% schulz"))
762    compare(title, target, actual)
763COMPARISON[('sasfit', 'sphere', 'schulz')] = compare_sasfit_sphere_schulz
764
765def compare_sasfit_ellipsoid_schulz():
766    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
767    #Effective radius =13.1353356684
768    #We have scaled the output from sasfit by 1e-4*volume*volfraction
769    #0.10050378152592121
770    pars = {
771        'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz',
772        'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz',
773        'sld': 4, 'sld_solvent': 1,
774        'volfraction': 0.3, 'radius_effective': 13.1353356684,
775    }
776
777    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQD.txt'))
778    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQSD.txt'))
779    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQBD.txt'))
780    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
781    actual = ellipsoid_pe(Q, norm="sasfit", **pars)
782    title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz"))
783    compare(title, target, actual)
784COMPARISON[('sasfit', 'ellipsoid', 'schulz')] = compare_sasfit_ellipsoid_schulz
785
786
787def compare_sasfit_ellipsoid_gaussian():
788    pars = {
789        'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian',
790        'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian',
791        'sld': 4, 'sld_solvent': 1,
792        'volfraction': 0, 'radius_effective': None,
793    }
794
795    #Rp=20,Re=10,eta_core=4,eta_solv=1
796    Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt'))
797    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None)
798    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
799    target = Theory(Q=Q, P=PQ0)
800    compare("sasfit ellipsoid no poly", target, actual); plt.show()
801
802    #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
803    Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt'))
804    pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None)
805    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
806    target = Theory(Q=Q, P=PQ_Rp10)
807    compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show()
808    #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
809    Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt'))
810    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None)
811    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
812    target = Theory(Q=Q, P=PQ_Re10)
813    compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show()
814    #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
815    Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt'))
816    pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None)
817    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
818    target = Theory(Q=Q, P=PQ_Rp30)
819    compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show()
820    #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
821    Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt'))
822    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None)
823    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
824    target = Theory(Q=Q, P=PQ_Re30)
825    compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show()
826    #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
827    Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt'))
828    pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None)
829    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
830    target = Theory(Q=Q, P=PQ_Rp60)
831    compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show()
832    #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
833    Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt'))
834    pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None)
835    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
836    target = Theory(Q=Q, P=PQ_Re60)
837    compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show()
838
839    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
840    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt'))
841    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'))
842    pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
843    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
844    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
845    compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show()
846    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
847    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'))
848    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'))
849    pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
850    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
851    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
852    compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show()
853    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
854    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'))
855    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'))
856    pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
857    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
858    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
859    compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show()
860
861    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
862    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'))
863    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'))
864    pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
865    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
866    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
867    compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show()
868    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
869    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'))
870    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'))
871    pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
872    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
873    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
874    compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show()
875    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
876    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'))
877    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'))
878    pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
879    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
880    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
881    compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show()
882
883    #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
884    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'))
885    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'))
886    pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254)
887    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
888    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
889    compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show()
890    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
891    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'))
892    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'))
893    pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149)
894    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
895    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
896    compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show()
897    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
898    Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'))
899    Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'))
900    pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917)
901    actual = ellipsoid_pe(Q, norm='sasfit', **pars)
902    target = Theory(Q=Q, S=SQ, Seff=SQ_EFF)
903    compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show()
904COMPARISON[('sasfit', 'ellipsoid', 'gaussian')] = compare_sasfit_ellipsoid_gaussian
905
906def main():
907    key = tuple(sys.argv[1:])
908    if key not in COMPARISON:
909        print("Usage: python sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid|vesicle|hollow_cylinder] [gaussian|schulz]")
910        print("But not for [yun sphere schulz], and vesicle & hollow_cylinder only with sasview.")
911        print("Note this compares chosen source against internal python code here, with adjustment for known scale etc issues.")
912        return
913    comparison = COMPARISON[key]
914    comparison()
915
916if __name__ == "__main__":
917    main()
Note: See TracBrowser for help on using the repository browser.