Changeset 01c8d9e in sasmodels for explore/beta/sasfit_compare.py
- Timestamp:
- Aug 7, 2018 12:32:18 PM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- c11d09f
- Parents:
- 707cbdb
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/beta/sasfit_compare.py
r707cbdb r01c8d9e 3 3 import sys,os 4 4 BETA_DIR = os.path.dirname(os.path.realpath(__file__)) 5 SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 5 #SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 6 SASMODELS_DIR = r"C:\Source\sasmodels" 6 7 sys.path.insert(0, SASMODELS_DIR) 7 8 import os 9 from collections import namedtuple 10 8 11 from matplotlib import pyplot as plt 9 12 import numpy as np … … 14 17 from scipy.special import gammaln # type: ignore 15 18 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. 18 ELLIPSOID = False 19 SPHERE = True 20 21 GAUSSIAN = True 22 SCHULZ = False 23 24 SASVIEW=True 25 SASFIT=True 26 YUN = True 27 28 def data_file(name): 29 return os.path.join(BETA_DIR + '\\data_files', name) 19 Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta') 20 Theory.__new__.__defaults__ = (None,) * len(Theory._fields) 30 21 31 22 #Used to calculate F(q) for the cylinder, sphere, ellipsoid models … … 47 38 with np.errstate(all='ignore'): 48 39 # GSL bessel_j1 taylor expansion 49 index = (x < 0.25) 40 index = (x < 0.25) 50 41 y = x[index]**2 51 42 c1 = -1.0/10.0 … … 71 62 calculator = DirectModel(data, model,cutoff=0) 72 63 calculator.pars = pars.copy() 73 calculator.pars.setdefault('background', 1e-3)64 calculator.pars.setdefault('background', 0) 74 65 return calculator 75 66 76 67 #gives the hardsphere structure factor that sasview uses 77 def hardsphere_simple(q, radius_effective, volfraction):78 CUTOFFHS=0.05 68 def _hardsphere_simple(q, radius_effective, volfraction): 69 CUTOFFHS=0.05 79 70 if fabs(radius_effective) < 1.E-12: 80 71 HARDSPH=1.0 … … 96 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 97 88 HARDSPH= 1./(1. + volfraction*FF ) 98 return HARDSPH 89 return HARDSPH 99 90 X4=X2*X2 100 91 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 ) ;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 ) 103 94 return HARDSPH 95 96 def hardsphere_simple(q, radius_effective, volfraction): 97 SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q] 98 return np.array(SQ) 104 99 105 100 #Used in gaussian quadrature for polydispersity 106 101 #returns values and the probability of those values based on gaussian distribution 107 def gaussian_distribution(center, sigma,lb,ub): 108 #3 standard deviations covers approx. 99.7% 102 N_GAUSS = 35 103 NSIGMA_GAUSS = 3 104 def gaussian_distribution(center, sigma, lb, ub): 105 #3 standard deviations covers approx. 99.7% 109 106 if sigma != 0: 110 nsigmas =3111 x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num= 35)107 nsigmas = NSIGMA_GAUSS 108 x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS) 112 109 x= x[(x >= lb) & (x <= ub)] 113 110 px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) … … 116 113 return np.array([center]), np.array([1]) 117 114 115 N_SCHULZ = 80 116 NSIGMA_SCHULZ = 8 118 117 def schulz_distribution(center, sigma, lb, ub): 119 118 if sigma != 0: 120 nsigmas =8121 x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num= 80)119 nsigmas = NSIGMA_SCHULZ 120 x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ) 122 121 x= x[(x >= lb) & (x <= ub)] 123 122 R = x/center … … 160 159 #SQ is monodisperse approach for structure factor 161 160 #SQ_EFF is the effective structure factor from beta approx 162 def ellipsoid_Theta(q, radius_polar, radius_equatorial, sld, sld_solvent,volfraction=0,effective_radius=0): 161 def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent, 162 volfraction=0, radius_effective=None): 163 163 #creates values z and corresponding probabilities w from legendre-gauss quadrature 164 volume = ellipsoid_volume(radius_polar, radius_equatorial) 164 165 z, w = leggauss(76) 165 166 F1 = np.zeros_like(q) 166 167 F2 = np.zeros_like(q) 167 168 #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 #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use 169 170 #legendre-gauss quadrature 170 171 for k, qk in enumerate(q): 171 172 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) 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 176 #the 1/2 comes from the change of variables mentioned above 177 177 F2 = F2/2.0 178 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 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 186 184 IQSM = volfraction*IQM*SQ 187 185 IQBM = volfraction*IQM*SQ_EFF 188 return F1, F2, IQM, IQSM, IQBM, SQ, SQ_EFF189 190 #IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc. 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. 191 189 #IQBD HAS NOT BEEN CROSS CHECKED AT ALL 192 def 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) 190 def 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 203 208 radius_eff = 0 209 F1, F2 = np.zeros_like(q), np.zeros_like(q) 204 210 for k, Rpk in enumerate(Rp_val): 205 211 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 212 theory = ellipsoid_theta(q,Rpk,Rei,sld,sld_solvent) 213 volume = ellipsoid_volume(Rpk, Rei) 214 weight = Rp_prob[k]*Re_prob[i] 215 total_weight += weight 216 total_volume += weight*volume 217 F1 += theory.F1*weight 218 F2 += theory.F2*weight 219 radius_eff += weight*ER_ellipsoid(Rpk,Rei) 220 F1 /= total_weight 221 F2 /= total_weight 222 average_volume = total_volume/total_weight 223 if radius_effective is None: 224 radius_effective = radius_eff/total_weight 225 print("this is the effective radius for pure python",radius_effective) 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 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 print("\n\n\n\n\n this is F1" + str(F1)) 248 249 print("\n\n\n\n\n this is F2" + str(F2)) 250 print("\n\n\n\n\n this is SQ" + str(SQ)) 251 return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 225 252 226 253 #polydispersity for sphere 227 def sphere_r(q,radius,sld,sld_solvent,volfraction=0,radius_pd=0.1,distribution='gaussian',norm_type='volfraction'): 228 if distribution == 'gaussian': 254 def sphere_r(q,radius,sld,sld_solvent, 255 radius_pd=0.1, radius_pd_type='gaussian', 256 volfraction=0, radius_effective=None, 257 background=0, scale=1, 258 norm='sasview'): 259 if norm not in ['sasview', 'sasfit', 'yun']: 260 raise TypeError("unknown norm "+norm) 261 if radius_pd_type == 'gaussian': 229 262 radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf) 230 elif distribution== 'schulz':263 elif radius_pd_type == 'schulz': 231 264 radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf) 232 Normalization=0 265 total_weight = total_volume = 0 266 F1 = np.zeros_like(q) 233 267 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 for k, rk in enumerate(radius_val): 269 volume = 4./3.*pi*rk**3 270 total_weight += radius_prob[k] 271 total_volume += radius_prob[k]*volume 272 form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) 273 F2 += radius_prob[k]*form**2 274 F1 += radius_prob[k]*form 275 F1 /= total_weight 276 F2 /= total_weight 277 average_volume = total_volume/total_weight 278 279 if radius_effective is None: 280 radius_effective = radius 281 average_volume = total_volume/total_weight 282 if norm == 'sasfit': 283 IQD = F2 284 elif norm == 'sasview': 285 IQD = F2/average_volume*1e-4*volfraction 286 elif norm == 'yun': 287 F1 *= 1e-6 # Yun is using sld in 1/A^2, not 1e-6/A^2 288 F2 *= 1e-12 289 IQD = F2/average_volume*1e8*volfraction 290 SQ = hardsphere_simple(q, radius_effective, volfraction) 291 beta = F1**2/F2 292 SQ_EFF = 1 + beta*(SQ - 1) 293 IQSD = IQD*SQ 294 IQBD = IQD*SQ_EFF 295 return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 268 296 269 297 ############################################################################### … … 276 304 ############################################################################### 277 305 ############################################################################### 278 # SASVIEW 279 if (ELLIPSOID == True) & (GAUSSIAN == True) & (SASVIEW == True): 306 307 def popn(d, keys): 308 """ 309 Splits a dict into two, with any key of *d* which is in *keys* removed 310 from *d* and added to *b*. Returns *b*. 311 """ 312 b = {} 313 for k in keys: 314 try: 315 b[k] = d.pop(k) 316 except KeyError: 317 pass 318 return b 319 320 def sasmodels_theory(q, Pname, **pars): 321 """ 322 Call sasmodels to compute the model with and without a hard sphere 323 structure factor. 324 """ 325 #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()} 326 Ppars = pars.copy() 327 Spars = popn(Ppars, ['radius_effective', 'volfraction']) 328 Ipars = pars.copy() 329 330 # Autofill npts and nsigmas for the given pd type 331 for k, v in pars.items(): 332 if k.endswith("_pd_type"): 333 if v == "gaussian": 334 n, nsigmas = N_GAUSS, NSIGMA_GAUSS 335 elif v == "schulz": 336 n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ 337 Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n) 338 Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 339 Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n) 340 Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) 341 342 343 #Ppars['scale'] = Spars.get('volfraction', 1) 344 P = build_model(Pname, q) 345 S = build_model("hardsphere", q) 346 I = build_model(Pname+"@hardsphere", q) 347 Pq = P(**Ppars)*pars.get('volfraction', 1) 348 Sq = S(**Spars) 349 Iq = I(**Ipars) 350 #Iq = Pq*Sq*pars.get('volfraction', 1) 351 #Sq = Iq/Pq 352 #Iq = None#= Sq = None 353 r=I._kernel.results 354 return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 355 356 def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 357 """ 358 Plot fields in common between target and actual, along with relative error. 359 """ 360 available = [s for s in fields.split() 361 if getattr(target, s) is not None and getattr(actual, s) is not None] 362 rows = len(available) 363 for row, field in enumerate(available): 364 Q = target.Q 365 I1, I2 = getattr(target, field), getattr(actual, field) 366 plt.subplot(rows, 2, 2*row+1) 367 plt.loglog(Q, abs(I1), label="target "+field) 368 plt.loglog(Q, abs(I2), label="value "+field) 369 #plt.legend(loc="upper left", bbox_to_anchor=(1,1)) 370 plt.legend(loc='lower left') 371 plt.subplot(rows, 2, 2*row+2) 372 plt.semilogx(Q, I2/I1 - 1, label="relative error") 373 #plt.semilogx(Q, I1/I2 - 1, label="relative error") 374 plt.tight_layout() 375 plt.suptitle(title) 376 plt.show() 377 378 def data_file(name): 379 return os.path.join(BETA_DIR, 'data_files', name) 380 381 def load_sasfit(path): 382 data = np.loadtxt(path, dtype=str, delimiter=';').T 383 data = np.vstack((map(float, v) for v in data[0:2])) 384 return data 385 386 COMPARISON = {} # Type: Dict[(str,str,str)] -> Callable[(), None] 387 388 def compare_sasview_sphere(pd_type='schulz'): 389 q = np.logspace(-5, 0, 250) 390 model = 'sphere' 391 pars = dict( 392 radius=20,sld=4,sld_solvent=1, 393 background=0, 394 radius_pd=.1, radius_pd_type=pd_type, 395 volfraction=0.15, 396 #radius_effective=12.59921049894873, # equivalent average sphere radius 397 ) 398 target = sasmodels_theory(q, model, **pars) 399 actual = sphere_r(q, norm='sasview', **pars) 400 title = " ".join(("sasmodels", model, pd_type)) 401 compare(title, target, actual) 402 COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') 403 COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') 404 405 def compare_sasview_ellipsoid(pd_type='gaussian'): 280 406 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 366 if (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 396 if 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 407 model = 'ellipsoid' 408 pars = dict( 409 radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1, 410 background=0, 411 radius_polar_pd=.1, radius_polar_pd_type=pd_type, 412 radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 413 volfraction=0.15, 414 radius_effective=270.7543927018, 415 #radius_effective=12.59921049894873, 416 ) 417 target = sasmodels_theory(q, model, beta_mode=1, **pars) 418 actual = ellipsoid_pe(q, norm='sasview', **pars) 419 print(actual) 420 title = " ".join(("sasmodels", model, pd_type)) 421 compare(title, target, actual) 422 COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian') 423 COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz') 424 425 def compare_yun_ellipsoid_mono(): 426 pars = { 427 'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 428 'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 429 'sld': 2, 'sld_solvent': 1, 430 'volfraction': 0.15, 431 # Yun uses radius for same volume sphere for effective radius 432 # whereas sasview uses the average curvature. 433 'radius_effective': 12.59921049894873, 434 } 435 436 data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T 437 Q = data[0] 438 F1 = data[1] 439 P = data[3]*pars['volfraction'] 440 S = data[5] 441 Seff = data[6] 442 target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 443 actual = ellipsoid_pe(Q, norm='yun', **pars) 444 title = " ".join(("yun", "ellipsoid", "no pd")) 445 #compare(title, target, actual, fields="P S I Seff Ibeta") 446 compare(title, target, actual) 447 COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono 448 COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono 449 450 def compare_yun_sphere_gauss(): 451 # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100 452 # With pd = 0.1, that's 14 sigma and 1400 points. 453 pars = { 454 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', 455 'sld': 6, 'sld_solvent': 0, 456 'volfraction': 0.1, 457 } 458 459 data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'),skiprows=2).T 460 Q = data[0] 461 F1 = data[1] 462 P = data[3] 463 S = data[5] 464 Seff = data[6] 465 target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 466 actual = sphere_r(Q, norm='yun', **pars) 467 title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) 468 compare(title, target, actual) 469 data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'),skiprows=2).T 470 pars.update(radius_pd=0.15) 471 Q = data[0] 472 F1 = data[1] 473 P = data[3] 474 S = data[5] 475 Seff = data[6] 476 target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 477 actual = sphere_r(Q, norm='yun', **pars) 478 title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) 479 compare(title, target, actual) 480 COMPARISON[('yun','sphere','gaussian')] = compare_yun_sphere_gauss 481 482 483 def compare_sasfit_sphere_gauss(): 484 #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 485 pars = { 486 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', 487 'sld': 4, 'sld_solvent': 1, 488 'volfraction': 0.3, 489 } 490 491 Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt')) 492 Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt')) 493 Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt')) 494 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt')) 495 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt')) 496 target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) 497 actual = sphere_r(Q, norm="sasfit", **pars) 498 title = " ".join(("sasfit", "sphere", "pd=10% gaussian")) 499 compare(title, target, actual) 500 #compare(title, target, actual, fields="P") 501 COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss 502 503 def compare_sasfit_sphere_schulz(): 692 504 #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1 693 505 #We have scaled the output from sasfit by 1e-4*volume*volfraction 694 506 #0.10050378152592121 695 if (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 731 if (ELLIPSOID == True) & (SCHULZ == True) & (SASVIEW == True) & (SASFIT == True): 507 pars = { 508 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz', 509 'sld': 4, 'sld_solvent': 1, 510 'volfraction': 0.3, 511 } 512 513 Q, IQ = load_sasfit(data_file('richard_test.txt')) 514 Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 515 Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 516 target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 517 actual = sphere_r(Q, norm="sasfit", **pars) 518 title = " ".join(("sasfit", "sphere", "pd=10% schulz")) 519 compare(title, target, actual) 520 COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz 521 522 def compare_sasfit_ellipsoid_schulz(): 732 523 #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 771 if (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 822 if (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() 524 #Effective radius =13.1353356684 525 #We have scaled the output from sasfit by 1e-4*volume*volfraction 526 #0.10050378152592121 527 pars = { 528 'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz', 529 'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz', 530 'sld': 4, 'sld_solvent': 1, 531 'volfraction': 0.3, 'radius_effective': 13.1353356684, 532 } 533 534 Q, IQ = load_sasfit(data_file('richard_test4.txt')) 535 Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 536 Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 537 target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 538 actual = ellipsoid_pe(Q, norm="sasfit", **pars) 539 title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz")) 540 compare(title, target, actual) 541 COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz 542 543 544 def compare_sasfit_ellipsoid_gaussian(): 545 pars = { 546 'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', 547 'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', 548 'sld': 4, 'sld_solvent': 1, 549 'volfraction': 0, 'radius_effective': None, 550 } 551 552 #Rp=20,Re=10,eta_core=4,eta_solv=1 553 Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt')) 554 pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None) 555 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 556 target = Theory(Q=Q, P=PQ0) 557 compare("sasfit ellipsoid no poly", target, actual); plt.show() 558 559 #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 560 Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt')) 561 pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None) 562 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 563 target = Theory(Q=Q, P=PQ_Rp10) 564 compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show() 565 #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 566 Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt')) 567 pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None) 568 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 569 target = Theory(Q=Q, P=PQ_Re10) 570 compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show() 571 #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 572 Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt')) 573 pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None) 574 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 575 target = Theory(Q=Q, P=PQ_Rp30) 576 compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show() 577 #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 578 Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt')) 579 pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None) 580 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 581 target = Theory(Q=Q, P=PQ_Re30) 582 compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show() 583 #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 584 Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt')) 585 pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None) 586 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 587 target = Theory(Q=Q, P=PQ_Rp60) 588 compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show() 589 #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly 590 Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt')) 591 pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None) 592 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 593 target = Theory(Q=Q, P=PQ_Re60) 594 compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show() 595 596 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 597 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt')) 598 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt')) 599 pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 600 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 601 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 602 compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show() 603 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 604 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt')) 605 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt')) 606 pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 607 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 608 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 609 compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show() 610 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 611 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt')) 612 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt')) 613 pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 614 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 615 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 616 compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show() 617 618 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 619 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt')) 620 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt')) 621 pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 622 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 623 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 624 compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show() 625 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 626 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt')) 627 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt')) 628 pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 629 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 630 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 631 compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show() 632 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 633 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt')) 634 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt')) 635 pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 636 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 637 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 638 compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show() 639 640 #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 641 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt')) 642 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt')) 643 pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) 644 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 645 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 646 compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show() 647 #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 648 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt')) 649 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt')) 650 pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) 651 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 652 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 653 compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show() 654 #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 655 Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt')) 656 Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt')) 657 pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) 658 actual = ellipsoid_pe(Q, norm='sasfit', **pars) 659 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 660 compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show() 661 COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian 662 663 def main(): 664 key = tuple(sys.argv[1:]) 665 if key not in COMPARISON: 666 print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]") 667 return 668 comparison = COMPARISON[key] 669 comparison() 670 671 if __name__ == "__main__": 672 main()
Note: See TracChangeset
for help on using the changeset viewer.