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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since b1a0f3e was 707cbdb, checked in by Suczewski <ges3@…>, 6 years ago

beta approximation comparison with sasfit
  • Property mode set to 100644
File size: 41.6 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)
7import os
8from matplotlib import pyplot as plt
9import numpy as np
10from numpy import pi, sin, cos, sqrt, fabs
11from numpy.polynomial.legendre import leggauss
12from scipy.special import j1 as J1
13from numpy import inf
14from scipy.special import gammaln  # type: ignore
15
16# THE FOLLOWING 7 BOOLEANS TAILOR WHICH GRAPHS ARE PRINTED WHEN THE PROGRAM RUNS
17#RICHARD, YOU WILL WANT SASVIEW, SASFIT, SCHULZ, ELLIPSOID, AND SPHERE ALL TRUE.
18ELLIPSOID = False
19SPHERE = True
20
21GAUSSIAN = True
22SCHULZ = False
23
24SASVIEW=True
25SASFIT=True
26YUN = True
27
28def data_file(name):
29    return os.path.join(BETA_DIR + '\\data_files', name)
30
31#Used to calculate F(q) for the cylinder, sphere, ellipsoid models
32def sas_sinx_x(x):
33    with np.errstate(all='ignore'):
34        retvalue = sin(x)/x
35    retvalue[x == 0.] = 1.
36    return retvalue
37
38def sas_2J1x_x(x):
39    with np.errstate(all='ignore'):
40        retvalue = 2*J1(x)/x
41    retvalue[x == 0] = 1.
42    return retvalue
43
44def sas_3j1x_x(x):
45    """return 3*j1(x)/x"""
46    retvalue = np.empty_like(x)
47    with np.errstate(all='ignore'):
48        # GSL bessel_j1 taylor expansion
49        index = (x < 0.25) 
50        y = x[index]**2
51        c1 = -1.0/10.0
52        c2 =  1.0/280.0
53        c3 = -1.0/15120.0
54        c4 =  1.0/1330560.0
55        c5 = -1.0/172972800.0
56        retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5))))
57        index = ~index
58        y = x[index]
59        retvalue[index] = 3*(sin(y) - y*cos(y))/y**3
60    retvalue[x == 0.] = 1.
61    return retvalue
62
63#Used to cross check my models with sasview models
64def build_model(model_name, q, **pars):
65    from sasmodels.core import load_model_info, build_model as build_sasmodel
66    from sasmodels.data import empty_data1D
67    from sasmodels.direct_model import DirectModel
68    model_info = load_model_info(model_name)
69    model = build_sasmodel(model_info, dtype='double!')
70    data = empty_data1D(q)
71    calculator = DirectModel(data, model,cutoff=0)
72    calculator.pars = pars.copy()
73    calculator.pars.setdefault('background', 1e-3)
74    return calculator
75
76#gives the hardsphere structure factor that sasview uses
77def hardsphere_simple(q, radius_effective, volfraction): 
78    CUTOFFHS=0.05 
79    if fabs(radius_effective) < 1.E-12:
80        HARDSPH=1.0
81        return HARDSPH
82    X = 1.0/( 1.0 -volfraction)
83    D= X*X
84    A= (1.+2.*volfraction)*D
85    A *=A
86    X=fabs(q*radius_effective*2.0)
87    if X < 5.E-06:
88        HARDSPH=1./A
89        return HARDSPH
90    X2 =X*X
91    B = (1.0 +0.5*volfraction)*D
92    B *= B
93    B *= -6.*volfraction
94    G=0.5*volfraction*A
95    if X < CUTOFFHS:
96        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
97        HARDSPH= 1./(1. + volfraction*FF )
98        return HARDSPH   
99    X4=X2*X2
100    S, C = sin(X), cos(X)
101    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 ;
102    HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
103    return HARDSPH
104
105#Used in gaussian quadrature for polydispersity
106#returns values and the probability of those values based on gaussian distribution
107def gaussian_distribution(center, sigma,lb,ub):
108    #3 standard deviations covers approx. 99.7%
109    if sigma != 0:
110        nsigmas=3
111        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=35)
112        x= x[(x >= lb) & (x <= ub)]
113        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
114        return x, px
115    else:
116        return np.array([center]), np.array([1])
117
118def schulz_distribution(center, sigma, lb, ub):
119    if sigma != 0:
120        nsigmas=8
121        x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=80)
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,volfraction=0,effective_radius=0):
163    #creates values z and corresponding probabilities w from legendre-gauss quadrature
164    z, w = leggauss(76)
165    F1 = np.zeros_like(q)
166    F2 = np.zeros_like(q)
167    #use a u subsition(u=cos) and then u=(z+1)/2 to change integration from
168    #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use
169    #legendre-gauss quadrature
170    for k, qk in enumerate(q):
171        r = sqrt(radius_equatorial**2*(1-((z+1)/2)**2)+radius_polar**2*((z+1)/2)**2)
172        F2i = ((sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r))**2
173        F2[k] = np.sum(w*F2i)
174        F1i = (sld-sld_solvent)*ellipsoid_volume(radius_polar,radius_equatorial)*sas_3j1x_x(qk*r)
175        F1[k] = np.sum(w*F1i)
176    #the 1/2 comes from the change of variables mentioned above
177    F2 = F2/2.0
178    F1 = F1/2.0
179    if effective_radius == 0:
180        effective_radius = ER_ellipsoid(radius_polar,radius_equatorial)
181    else:
182        effective_radius = effective_radius
183    SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])   
184    SQ_EFF = 1 + F1**2/F2*(SQ - 1) 
185    IQM = 1e-4/ellipsoid_volume(radius_polar,radius_equatorial)*F2
186    IQSM = volfraction*IQM*SQ
187    IQBM = volfraction*IQM*SQ_EFF
188    return F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF
189
190#IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc.
191#IQBD HAS NOT BEEN CROSS CHECKED AT ALL
192def ellipsoid_pe(q, Rp, Re, sld, sld_solvent, volfraction=0,effective_radius=0,radius_polar_pd=0.1,radius_equatorial_pd=0.1,distribution='gaussian'):
193    if distribution == 'gaussian':
194        Rp_val, Rp_prob = gaussian_distribution(Rp, radius_polar_pd*Rp, 0, inf)
195        Re_val, Re_prob = gaussian_distribution(Re, radius_equatorial_pd*Re, 0, inf)
196    elif distribution == 'schulz':
197        Rp_val, Rp_prob = schulz_distribution(Rp, radius_polar_pd*Rp, 0, inf)
198        Re_val, Re_prob = schulz_distribution(Re, radius_equatorial_pd*Re, 0, inf)
199    Normalization = 0
200    total_weight = 0
201    PQ = np.zeros_like(q)
202    F12,F21 = np.zeros_like(q), np.zeros_like(q)
203    radius_eff = 0
204    for k, Rpk in enumerate(Rp_val):
205        for i, Rei in enumerate(Re_val):
206            F1i, F2i, PQi, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,Rpk,Rei,sld,sld_solvent)
207            radius_eff += Rp_prob[k]*Re_prob[i]*ER_ellipsoid(Rpk,Rei)
208            total_weight +=  Rp_prob[k]*Re_prob[i]
209            Normalization += Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei)
210            PQ += PQi*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei)
211            F12 += F1i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei)
212            F21 += F2i*Rp_prob[k]*Re_prob[i]*ellipsoid_volume(Rpk, Rei)
213    F12 = (F12/Normalization)**2
214    F21 = F21/Normalization
215    if effective_radius == 0:
216        effective_radius = radius_eff/total_weight
217    else:
218        effective_radius = effective_radius
219    SQ = np.array([hardsphere_simple(qk, effective_radius, volfraction) for qk in q])   
220    SQ_EFF = 1 + F12/F21*(SQ - 1)     
221    IQD = PQ/Normalization
222    IQSD = volfraction*IQD*SQ
223    IQBD = volfraction*IQD*SQ_EFF
224    return IQD, IQSD, IQBD, SQ, SQ_EFF
225
226#polydispersity for sphere
227def sphere_r(q,radius,sld,sld_solvent,volfraction=0,radius_pd=0.1,distribution='gaussian',norm_type='volfraction'):
228    if distribution == 'gaussian':   
229        radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf)
230    elif distribution == 'schulz':
231        radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf)
232    Normalization=0
233    F2 = np.zeros_like(q)
234    F1 = np.zeros_like(q)
235    if norm_type == 'numdensity':
236        for k, rk in enumerate(radius_val):
237            volume = 4./3.*pi*rk**3
238            Normalization += radius_prob[k]*volume
239            F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2/volume
240            F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)/volume
241           
242            F2 += radius_prob[k]*volume*F2i
243            F1 += radius_prob[k]*volume*F1i
244        F21 = F2/Normalization
245        F12 = (F1/Normalization)**2
246        SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q])
247        SQ_EFF = 1 + F12/F21*(SQ - 1) 
248        IQD = 1e-4*F21
249        IQSD = volfraction*IQD*SQ
250        IQBD = volfraction*IQD*SQ_EFF
251    elif norm_type == 'volfraction':
252        for k, rk in enumerate(radius_val):
253            volume = 4./3.*pi*rk**3
254            Normalization += radius_prob[k]
255            F2i = ((sld-sld_solvent)*volume*sas_3j1x_x(q*rk))**2
256            F1i = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk)
257            F2 += radius_prob[k]*F2i
258            F1 += radius_prob[k]*F1i
259   
260        F21 = F2/Normalization
261        F12 = (F1/Normalization)**2
262        SQ = np.array([hardsphere_simple(qk, radius, volfraction) for qk in q])
263        SQ_EFF = 1 + F12/F21*(SQ - 1) 
264        IQD = 1e-4/(4./3.*pi*radius**3)*F21
265        IQSD = volfraction*IQD*SQ
266        IQBD = volfraction*IQD*SQ_EFF
267    return IQD, IQSD, IQBD, SQ, SQ_EFF
268
269###############################################################################
270###############################################################################
271###############################################################################
272##################                                           ##################
273##################                   TESTS                   ##################
274##################                                           ##################
275###############################################################################
276###############################################################################
277###############################################################################
278# SASVIEW   
279if (ELLIPSOID == True) & (GAUSSIAN == True) & (SASVIEW == True):
280    q = np.logspace(-5, 0, 50)
281    volfraction = 0.2
282    model = build_model("ellipsoid",q)
283    Sq = model(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,\
284               background=0,radius_polar_pd=.1,radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35)
285    IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q,20,400,4,1,volfraction)
286    Sq2 = model(radius_polar=20,radius_equatorial=10,sld=2,sld_solvent=1,background=0)
287    IQM = ellipsoid_Theta(q,20,10,2,1)[2]
288    model2 = build_model("ellipsoid@hardsphere", q)
289    Sq3 = model2(radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\
290                 radius_polar_pd_n=35,radius_equatorial_pd=.1,radius_equatorial_pd_n=35)
291    Sqp = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\
292                   background=0,radius_polar_pd=0.1,radius_polar_pd_n=35,radius_equatorial_pd=0,radius_equatorial_pd_n=35)
293    IQDp = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0] 
294    Sqp2 = model(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,\
295                   background=0,radius_polar_pd=0,radius_polar_pd_n=35,radius_equatorial_pd=0.1,radius_equatorial_pd_n=35)
296    IQDe = ellipsoid_pe(q,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]
297    print('\n\n\n\n ELLIPSOID COMPARISON WITH SASVIEW WITHOUT V')
298    plt.subplot(323)
299    plt.loglog(q, Sq,'r')
300    plt.loglog(q,IQD,'b')
301    plt.title('IQD')
302    plt.subplot(324)
303    plt.semilogx(q,Sq/IQD-1)
304    plt.subplot(321)
305    plt.title('IQM')
306    plt.loglog(q, Sq2,'r')
307    plt.loglog(q,IQM,'b')
308    plt.subplot(322)
309    plt.semilogx(q,IQM/Sq2-1)
310    plt.subplot(325)
311    plt.title('IQSD')
312    plt.loglog(q, Sq3, '-b')
313    plt.loglog(q, IQSD, '-r')
314    plt.subplot(326)
315    plt.semilogx(q, Sq3/IQSD-1, 'b')
316    plt.tight_layout()
317    plt.show()
318    plt.subplot(221)
319    plt.loglog(q,IQDp,'r')
320    plt.loglog(q,Sqp,'b')
321    plt.title('IQD Polar')
322    plt.subplot(222)
323    plt.plot(q,IQDp/Sqp-1)
324    plt.subplot(223)
325    plt.loglog(q,IQDe,'r')
326    plt.loglog(q,Sqp2,'b')
327    plt.title('IQD Equatorial')
328    plt.subplot(224)
329    plt.plot(q,IQDe/Sqp2-1)
330    plt.tight_layout()
331    plt.show()
332#comparing monodisperse Beta approximation to version without
333    q=np.linspace(1e-5,1)
334    F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF = ellipsoid_Theta(q,20,10,2,1,0.15,12.59921049894873)
335    plt.subplot(321)
336    plt.loglog(q,IQBM,'g',label='IQBM')
337    plt.loglog(q,IQSM,'r',label= 'IQSM')
338    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
339    plt.subplot(323)
340    plt.plot(q,IQBM,'g',label='IQBM')
341    plt.plot(q,IQSM,'r',label= 'IQSM')
342    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
343    plt.subplot(325)
344    plt.plot(q,IQBM/IQSM,'r',label= 'IQBM/IQSM')
345    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
346    plt.tight_layout()
347    plt.show()
348#comparing polydispersed Beta approximation to version without
349    IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(q, 20,10, 2,1, 0.15,12.59921049894873)
350    plt.subplot(321)
351    plt.loglog(q, SQ)
352    plt.loglog(q, SQ_EFF,label='SQ, SQ_EFF')
353    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
354    plt.subplot(323)
355    plt.loglog(q,IQBD)
356    plt.loglog(q,IQSD,label='IQBD,IQSD')
357    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
358    plt.subplot(325)
359    plt.plot(q,IQBD)
360    plt.plot(q,IQSD,label='IQBD,IQSD')
361    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
362    plt.tight_layout()
363    plt.show()
364
365# YUN
366if (ELLIPSOID == True) & (GAUSSIAN == True) & (YUN == True):
367    #Yun's data for Beta Approximation
368    data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T   
369    print('\n\n ELLIPSOID COMPARISON WITH YUN WITHOUT V')
370    Q=data[0]
371    F1,F2,IQM,IQSM,IQBM,SQ,SQ_EFF = ellipsoid_Theta(Q,20,10,2,1,0.15,12.59921049894873)
372    plt.subplot(211)
373    plt.loglog(Q,0.15*data[3]*data[5]*ellipsoid_volume(20,10)*1e-4)
374    plt.loglog(Q,IQSM,'-g',label='IQSM')
375    plt.loglog(data[0], data[3]*ellipsoid_volume(20,10)*1e-4,'-k')
376    plt.loglog(Q,IQM,'r',label = 'IQM')
377    plt.ylim(1e-5,4)
378    plt.legend(loc="upper left", bbox_to_anchor=(1,1))
379    plt.show()
380    plt.subplot(221)
381    plt.loglog(data[0],SQ)
382    plt.loglog(data[0],data[5])
383    plt.title('SQ')
384    plt.subplot(222)
385    plt.semilogx(data[0],SQ/data[5]-1)
386    plt.subplot(223)
387    plt.title('SQ_EFF')
388    plt.loglog(data[0],SQ_EFF)
389    plt.loglog(data[0],data[6])
390    plt.subplot(224)
391    plt.semilogx(data[0],(SQ_EFF/data[6])-1,label='Sq effective ratio')
392    plt.tight_layout()
393    plt.show()
394
395#SASFIT
396if ELLIPSOID and GAUSSIAN and SASFIT:
397    print('\n\n ELLIPSOID COMPARISON WITH SASFIT WITHOUT V')
398    #Rp=20,Re=10,eta_core=4,eta_solv=1
399    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQM.txt'),dtype=str,delimiter=';').T
400    sasqdata=map(float,sasdata[0])
401    sasvaluedata=map(float,sasdata[1])
402    plt.subplot(321)
403    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0],'b')
404    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata),'g')
405    plt.title('IQM')
406    plt.subplot(322)
407    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0)[0]-1)
408    plt.tight_layout()
409    plt.show()
410    print('========================================================================================')
411    print('========================================================================================')
412#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
413    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq.txt'),dtype=str,delimiter=';').T
414    sasqdata=map(float,sasdata[0])
415    sasvaluedata=map(float,sasdata[1])
416    plt.subplot(321)
417    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3])
418    plt.loglog(sasqdata,np.array(sasvaluedata))
419    plt.title('SQ poly 10% 0.15')
420    plt.subplot(322)
421    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1)
422#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
423    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq2.txt'),dtype=str,delimiter=';').T
424    sasqdata=map(float,sasdata[0])
425    sasvaluedata=map(float,sasdata[1])
426    plt.subplot(323)
427    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3])
428    plt.loglog(sasqdata,np.array(sasvaluedata))
429    plt.title('SQ poly 30% .15')
430    plt.subplot(324)
431    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1)
432#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
433    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq3.txt'),dtype=str,delimiter=';').T
434    sasqdata=map(float,sasdata[0])
435    sasvaluedata=map(float,sasdata[1])
436    plt.subplot(325)
437    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3])
438    plt.loglog(sasqdata,np.array(sasvaluedata))
439    plt.title('SQ poly 60% .15')
440    plt.subplot(326)
441    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1)
442    plt.tight_layout()
443    plt.show() 
444    print('========================================================================================')
445    print('========================================================================================')
446#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
447    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq4.txt'),dtype=str,delimiter=';').T
448    sasqdata=map(float,sasdata[0])
449    sasvaluedata=map(float,sasdata[1])
450    plt.subplot(321)
451    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3])
452    plt.loglog(sasqdata,np.array(sasvaluedata))
453    plt.title('SQ poly 10% .3')
454    plt.subplot(322)
455    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1)
456#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
457    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq5.txt'),dtype=str,delimiter=';').T
458    sasqdata=map(float,sasdata[0])
459    sasvaluedata=map(float,sasdata[1])
460    plt.subplot(323)
461    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3])
462    plt.loglog(sasqdata,np.array(sasvaluedata))
463    plt.title('SQ poly 30% .3')
464    plt.subplot(324)
465    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1)
466#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
467    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq6.txt'),dtype=str,delimiter=';').T
468    sasqdata=map(float,sasdata[0])
469    sasvaluedata=map(float,sasdata[1])
470    plt.subplot(325)
471    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3])
472    plt.loglog(sasqdata,np.array(sasvaluedata))
473    plt.title('SQ poly 60% .3')
474    plt.subplot(326)
475    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1)
476    plt.tight_layout()
477    plt.show()
478    print('========================================================================================')
479    print('========================================================================================')
480#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
481    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq7.txt'),dtype=str,delimiter=';').T
482    sasqdata=map(float,sasdata[0])
483    sasvaluedata=map(float,sasdata[1])
484    plt.subplot(321)
485    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3])
486    plt.loglog(sasqdata,np.array(sasvaluedata))
487    plt.title('SQ poly 10% .6')
488    plt.subplot(322)
489    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[3]-1)
490#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
491    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq8.txt'),dtype=str,delimiter=';').T
492    sasqdata=map(float,sasdata[0])
493    sasvaluedata=map(float,sasdata[1])
494    plt.subplot(323)
495    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3])
496    plt.loglog(sasqdata,np.array(sasvaluedata))
497    plt.title('SQ poly 30% .6')
498    plt.subplot(324)
499    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[3]-1)
500#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
501    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sq9.txt'),dtype=str,delimiter=';').T
502    sasqdata=map(float,sasdata[0])
503    sasvaluedata=map(float,sasdata[1])
504    plt.subplot(325)
505    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3])
506    plt.loglog(sasqdata,np.array(sasvaluedata))
507    plt.title('SQ poly 60% .6')
508    plt.subplot(326)
509    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[3]-1)
510    plt.tight_layout()
511    plt.show()
512    print('========================================================================================')
513    print('========================================================================================')
514#checks beta structre factor
515#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15
516    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt'),dtype=str,delimiter=';').T
517    sasqdata=map(float,sasdata[0])
518    sasvaluedata=map(float,sasdata[1])
519    plt.subplot(321)
520    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4])
521    plt.loglog(sasqdata,np.array(sasvaluedata))
522    plt.title('SQ_EFF poly 10% 0.15')
523    plt.subplot(322)
524    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1)
525#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15
526    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt'),dtype=str,delimiter=';').T
527    sasqdata=map(float,sasdata[0])
528    sasvaluedata=map(float,sasdata[1])
529    plt.subplot(323)
530    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4])
531    plt.loglog(sasqdata,np.array(sasvaluedata))
532    plt.title('SQ_EFF poly 30% .15')
533    plt.subplot(324)
534    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1)
535#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15
536    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt'),dtype=str,delimiter=';').T
537    sasqdata=map(float,sasdata[0])
538    sasvaluedata=map(float,sasdata[1])
539    plt.subplot(325)
540    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4])
541    plt.loglog(sasqdata,np.array(sasvaluedata))
542    plt.title('SQ_EFF poly 60% .15')
543    plt.subplot(326)
544    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,.15,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1)
545    plt.tight_layout()
546    plt.show() 
547    print('========================================================================================')
548    print('========================================================================================')
549#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3
550    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt'),dtype=str,delimiter=';').T
551    sasqdata=map(float,sasdata[0])
552    sasvaluedata=map(float,sasdata[1])
553    plt.subplot(321)
554    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4])
555    plt.loglog(sasqdata,np.array(sasvaluedata))
556    plt.title('SQ_EFF poly 10% .3')
557    plt.subplot(322)
558    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1)
559#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3
560    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt'),dtype=str,delimiter=';').T
561    sasqdata=map(float,sasdata[0])
562    sasvaluedata=map(float,sasdata[1])
563    plt.subplot(323)
564    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4])
565    plt.loglog(sasqdata,np.array(sasvaluedata))
566    plt.title('SQ_EFF poly 30% .3')
567    plt.subplot(324)
568    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1)
569#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3
570    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt'),dtype=str,delimiter=';').T
571    sasqdata=map(float,sasdata[0])
572    sasvaluedata=map(float,sasdata[1])
573    plt.subplot(325)
574    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4])
575    plt.loglog(sasqdata,np.array(sasvaluedata))
576    plt.title('SQ_EFF poly 60% .3')
577    plt.subplot(326)
578    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.3,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1)
579    plt.tight_layout()
580    plt.show()
581    print('========================================================================================')
582    print('========================================================================================')
583#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6
584    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt'),dtype=str,delimiter=';').T
585    sasqdata=map(float,sasdata[0])
586    sasvaluedata=map(float,sasdata[1])
587    plt.subplot(321)
588    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4])
589    plt.loglog(sasqdata,np.array(sasvaluedata))
590    plt.title('SQ_EFF poly 10% .6')
591    plt.subplot(322)
592    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.1,radius_equatorial_pd=0)[4]-1)
593#N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6
594    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt'),dtype=str,delimiter=';').T
595    sasqdata=map(float,sasdata[0])
596    sasvaluedata=map(float,sasdata[1])
597    plt.subplot(323)
598    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4])
599    plt.loglog(sasqdata,np.array(sasvaluedata))
600    plt.title('SQ_EFF poly 30% .6')
601    plt.subplot(324)
602    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.3,radius_equatorial_pd=0)[4]-1)
603#N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6
604    sasdata = np.loadtxt(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt'),dtype=str,delimiter=';').T
605    sasqdata=map(float,sasdata[0])
606    sasvaluedata=map(float,sasdata[1])
607    plt.subplot(325)
608    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4])
609    plt.loglog(sasqdata,np.array(sasvaluedata))
610    plt.title('SQ_EFF poly 60% .6')
611    plt.subplot(326)
612    plt.plot(sasqdata,np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,0.6,radius_polar_pd=0.6,radius_equatorial_pd=0)[4]-1)
613    plt.tight_layout()
614    plt.show()
615    print('========================================================================================')
616    print('========================================================================================')
617#checks polydispersion for single parameter and varying pd with sasfit
618#N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
619    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD.txt'),dtype=str,delimiter=';').T
620    sasqdata=map(float,sasdata[0])
621    sasvaluedata=map(float,sasdata[1])
622    plt.subplot(221)
623    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0])
624    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
625    plt.title('IQD polar 10%')
626    plt.subplot(222)
627    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.1,radius_equatorial_pd=0)[0]-1)
628#N=1,s=1,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
629    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD2.txt'),dtype=str,delimiter=';').T
630    sasqdata=map(float,sasdata[0])
631    sasvaluedata=map(float,sasdata[1])
632    plt.subplot(223)
633    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0])
634    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
635    plt.title('IQD equatorial 10%')
636    plt.subplot(224)
637    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.1)[0]-1)
638    plt.tight_layout()
639    plt.show()
640    print('========================================================================================')
641    print('========================================================================================')
642    #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
643    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD3.txt'),dtype=str,delimiter=';').T
644    sasqdata=map(float,sasdata[0])
645    sasvaluedata=map(float,sasdata[1])
646    plt.subplot(221)
647    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0])
648    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
649    plt.title('IQD polar 30%')
650    plt.subplot(222)
651    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.3,radius_equatorial_pd=0)[0]-1)
652#N=1,s=3,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
653    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD4.txt'),dtype=str,delimiter=';').T
654    sasqdata=map(float,sasdata[0])
655    sasvaluedata=map(float,sasdata[1])
656    plt.subplot(223)
657    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0])
658    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
659    plt.title('IQD equatorial 30%')
660    plt.subplot(224)
661    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.3)[0]-1)
662    plt.tight_layout()
663    plt.show()
664    print('========================================================================================')
665    print('========================================================================================')
666    #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
667    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD5.txt'),dtype=str,delimiter=';').T
668    sasqdata=map(float,sasdata[0])
669    sasvaluedata=map(float,sasdata[1])
670    plt.subplot(221)
671    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0])
672    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
673    plt.title('IQD polar 60%')
674    plt.subplot(222)
675    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0.6,radius_equatorial_pd=0)[0]-1)
676#N=1,s=6,X0=10,distr equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly
677    sasdata = np.loadtxt(data_file('sasfit_ellipsoid_IQD6.txt'),dtype=str,delimiter=';').T
678    sasqdata=map(float,sasdata[0])
679    sasvaluedata=map(float,sasdata[1])
680    plt.subplot(223)
681    plt.loglog(sasqdata, ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0])
682    plt.loglog(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata))
683    plt.title('IQD equatorial 60%')
684    plt.subplot(224)
685    plt.plot(sasqdata,1e-4/ellipsoid_volume(20,10)*np.array(sasvaluedata)/ellipsoid_pe(sasqdata,20,10,4,1,radius_polar_pd=0,radius_equatorial_pd=0.6)[0]-1)
686    plt.tight_layout()
687    plt.show()
688    print('========================================================================================')
689    print('========================================================================================')
690
691# TEST FOR RICHARD
692    #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1
693    #We have scaled the output from sasfit by 1e-4*volume*volfraction
694    #0.10050378152592121
695if (SPHERE == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True):
696    print('                   *****SPHERE MODEL*****')
697    sasdata = np.loadtxt(data_file('richard_test.txt'),dtype=str,delimiter=';').T
698    sasqdata=map(float,sasdata[0])
699    model=build_model('sphere',sasqdata)
700    Iq = model(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8)
701    model2=build_model('sphere@hardsphere',sasqdata)
702    IQS = model2(radius=20,radius_pd=0.1,radius_pd_n=80,sld=4,sld_solvent=1,background=0,radius_pd_type='schulz',radius_pd_nsigma=8,volfraction=0.3)
703   
704    sasvaluedata=map(float,sasdata[1])
705    sasdata2 = np.loadtxt(data_file('richard_test2.txt'),dtype=str,delimiter=';').T
706    sasqdata2=map(float,sasdata2[0])
707    sasvaluedata2=map(float,sasdata2[1])
708    sasdata3 = np.loadtxt(data_file('richard_test3.txt'),dtype=str,delimiter=';').T
709    sasqdata3=map(float,sasdata3[0])
710    sasvaluedata3=map(float,sasdata3[1])
711    IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3,0.1,distribution='schulz')
712    plt.grid(True)
713    plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)')
714    plt.loglog(sasqdata,Iq,'b')
715    plt.loglog(sasqdata,IQD,'r')
716    plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*np.array(sasvaluedata),'g')
717    plt.show()
718    plt.grid(True)
719    plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)')
720    plt.loglog(sasqdata,IQSD,'b')
721    plt.loglog(sasqdata,IQS,'r')
722    plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata2),'g')
723    plt.show()
724    plt.grid(True)
725    plt.title('IQBD(blue) vs. SASFIT(green)')
726    plt.loglog(sasqdata,IQBD,'b')
727    plt.loglog(sasqdata,1e-4/(4./3.*pi*20**3)*0.3*np.array(sasvaluedata3),'g')
728    plt.show()
729
730# TEST FOR RICHARD
731if (ELLIPSOID == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True):   
732    #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1
733         #Effective radius =13.1353356684
734        #We have scaled the output from sasfit by 1e-4*volume*volfraction
735        #0.10050378152592121
736    print('                   *****ELLIPSOID MODEL*****')
737    sasdata = np.loadtxt(data_file('richard_test4.txt'),dtype=str,delimiter=';').T
738    sasqdata=map(float,sasdata[0])
739    model=build_model('ellipsoid',sasqdata)
740    Iq = model(radius_polar=20,radius_polar_pd=0.1,radius_polar_pd_n=80, radius_equatorial=10, radius_equatorial_pd=0, sld=4,sld_solvent=1,background=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8)
741    model2=build_model('ellipsoid@hardsphere',sasqdata)
742    IQS =  model2(radius_polar=20,radius_equatorial=10,sld=4,sld_solvent=1,background=0,radius_polar_pd=.1,\
743                     radius_polar_pd_n=80,radius_equatorial_pd=0,radius_polar_pd_type='schulz',radius_polar_pd_nsigma=8,volfraction=0.3)
744    sasvaluedata=map(float,sasdata[1])
745    sasdata2 = np.loadtxt(data_file('richard_test5.txt'),dtype=str,delimiter=';').T
746    sasqdata2=map(float,sasdata2[0])
747    sasvaluedata2=map(float,sasdata2[1])
748    sasdata3 = np.loadtxt(data_file('richard_test6.txt'),dtype=str,delimiter=';').T
749    sasqdata3=map(float,sasdata3[0])
750    sasvaluedata3=map(float,sasdata3[1])
751    IQD, IQSD, IQBD, SQ, SQ_EFF = ellipsoid_pe(np.array(sasqdata),20,10,4,1,.3,radius_polar_pd=0.1,radius_equatorial_pd=0,distribution='schulz')   
752    plt.grid(True)
753    plt.title('IQD(blue) vs. SASVIEW(red) vs. SASFIT(green)')
754    plt.loglog(sasqdata,Iq,'b')
755    plt.loglog(sasqdata,IQD,'r')
756    plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*np.array(sasvaluedata),'g')
757    plt.show()
758    plt.grid(True)
759    plt.title('IQSD(blue) vs. SASVIEW(red) vs. SASFIT(green)')
760    plt.loglog(sasqdata,IQSD,'b')
761    plt.loglog(sasqdata,IQS,'r')
762    plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata2),'g')
763    plt.show()
764    plt.grid(True)
765    plt.title('IQBD(blue) vs. SASFIT(green)')
766    plt.loglog(sasqdata,IQBD,'b')
767    plt.loglog(sasqdata,1e-4/(4./3.*pi*20*10**2)*0.3*np.array(sasvaluedata3),'g')
768    plt.show()   
769   
770# SASVIEW GAUSS
771if (SPHERE == True) & (GAUSSIAN == True) & (SASVIEW == True):
772    q = np.logspace(-5, 0, 200)
773    IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,.3)
774    model = build_model("sphere", q)
775    Sq = model(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0)
776    model2 = build_model("sphere@hardsphere", q)
777    S=build_model("hardsphere",q)(radius_effective=20,volfraction=.3,background=0)
778    Sq2 = model2(radius=20,radius_pd=0.1,radius_pd_n=35,sld=4,sld_solvent=1,background=0,volfraction=.3)
779    print('\n\n SPHERE COMPARISON WITH SASVIEW WITHOUT V')
780    plt.subplot(221)
781    plt.title('IQD')
782    plt.loglog(q, IQD, '-b')
783    plt.loglog(q, Sq, '-r')
784    plt.subplot(222)
785    plt.semilogx(q, Sq/IQD-1, '-g')
786    plt.tight_layout()
787    plt.show()
788    plt.subplot(221)
789    plt.title('SQ')
790    plt.plot(q, SQ, '-r')
791    plt.plot(q,S,'-k')
792    plt.subplot(222)
793    plt.plot(SQ/S-1)
794    plt.tight_layout()
795    plt.show()
796    plt.subplot(221)
797    plt.title('IQSD')
798    plt.loglog(q, IQSD, '-b')
799    plt.loglog(q, Sq2, '-r',label='IQSD')
800    plt.subplot(222)
801    plt.semilogx(q, Sq2/IQSD-1, '-g',label='IQSD ratio')
802    plt.tight_layout()
803    plt.show()
804    IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(q,20,4,1,0.15)
805    plt.subplot(211)
806    plt.title('SQ vs SQ_EFF')
807    plt.plot(q, SQ)
808    plt.plot(q, SQ_EFF)
809    plt.tight_layout()
810    plt.show() 
811    plt.subplot(221)
812    plt.title('IQSD vs IQBD')
813    plt.loglog(q,IQBD)
814    plt.loglog(q,IQSD,label='IQBD,IQSD')
815    plt.subplot(222)
816    plt.plot(q,IQBD)
817    plt.plot(q,IQSD,)
818    plt.tight_layout()
819    plt.show() 
820
821# SASFIT GAUSS
822if (SPHERE == True) & (GAUSSIAN == True) & (SASFIT == True):
823    print('\n\n SPHERE COMPARISON WITH SASFIT WITHOUT V')
824#N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3
825    sasdata = np.loadtxt(data_file('sasfit_sphere_IQD.txt'),dtype=str,delimiter=';').T
826    sasqdata=map(float,sasdata[0])
827    sasvaluedata=map(float,sasdata[1])
828    IQD, IQSD, IQBD, SQ, SQ_EFF = sphere_r(np.array(sasqdata),20,4,1,.3)
829    plt.subplot(221)
830    plt.title('IQD')
831    plt.loglog(sasqdata,IQD )
832    plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata))
833    plt.subplot(222)
834    plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*np.array(sasvaluedata)/IQD-1)
835    plt.tight_layout()
836    plt.show()
837    sasdata = np.loadtxt(data_file('sasfit_sphere_IQSD.txt'),dtype=str,delimiter=';').T
838    sasqdata=map(float,sasdata[0])
839    sasvaluedata=map(float,sasdata[1])
840    plt.subplot(221)
841    plt.title('IQSD')
842    plt.loglog(sasqdata, IQSD)
843    plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata))
844    plt.subplot(222)
845    plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQSD-1)
846    plt.tight_layout()
847    plt.show()
848    sasdata = np.loadtxt(data_file('sasfit_sphere_IQBD.txt'),dtype=str,delimiter=';').T
849    sasqdata=map(float,sasdata[0])
850    sasvaluedata=map(float,sasdata[1])
851    plt.subplot(221)
852    plt.title('IQBD')
853    plt.loglog(sasqdata,IQBD)
854    plt.loglog(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata))
855    plt.subplot(222)
856    plt.semilogx(sasqdata,1e-4/(4./3*pi*20**3)*0.3*np.array(sasvaluedata)/IQBD-1)
857    plt.tight_layout()
858    plt.show()
859    sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sq.txt'),dtype=str,delimiter=';').T
860    sasqdata=map(float,sasdata[0])
861    sasvaluedata=map(float,sasdata[1])
862    plt.subplot(221)
863    plt.title('SQ')
864    plt.loglog(sasqdata, SQ)
865    plt.loglog(sasqdata,np.array(sasvaluedata))
866    plt.subplot(222)
867    plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ-1)
868    plt.tight_layout()
869    plt.show()
870    sasdata = np.loadtxt(data_file('sasfit_polydisperse_sphere_sqeff.txt'),dtype=str,delimiter=';').T
871    sasqdata=map(float,sasdata[0])
872    sasvaluedata=map(float,sasdata[1])
873    plt.subplot(221)
874    plt.title('SQ_EFF')
875    plt.loglog(sasqdata,SQ_EFF)
876    plt.loglog(sasqdata,np.array(sasvaluedata))
877    plt.subplot(222)
878    plt.semilogx(sasqdata,np.array(sasvaluedata)/SQ_EFF-1)
879    plt.tight_layout()
880    plt.show()
Note: See TracBrowser for help on using the repository browser.