[707cbdb] | 1 | from __future__ import division, print_function |
---|
| 2 | # Make sasmodels available on the path |
---|
[2cefd79] | 3 | import sys, os |
---|
[707cbdb] | 4 | BETA_DIR = os.path.dirname(os.path.realpath(__file__)) |
---|
| 5 | SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) |
---|
| 6 | sys.path.insert(0, SASMODELS_DIR) |
---|
[2cefd79] | 7 | |
---|
| 8 | from collections import namedtuple |
---|
| 9 | |
---|
[707cbdb] | 10 | from matplotlib import pyplot as plt |
---|
| 11 | import numpy as np |
---|
| 12 | from numpy import pi, sin, cos, sqrt, fabs |
---|
| 13 | from numpy.polynomial.legendre import leggauss |
---|
| 14 | from scipy.special import j1 as J1 |
---|
| 15 | from numpy import inf |
---|
| 16 | from scipy.special import gammaln # type: ignore |
---|
| 17 | |
---|
[2cefd79] | 18 | Theory = namedtuple('Theory', 'Q F1 F2 P S I Seff Ibeta') |
---|
| 19 | Theory.__new__.__defaults__ = (None,) * len(Theory._fields) |
---|
[707cbdb] | 20 | |
---|
| 21 | #Used to calculate F(q) for the cylinder, sphere, ellipsoid models |
---|
| 22 | def sas_sinx_x(x): |
---|
| 23 | with np.errstate(all='ignore'): |
---|
| 24 | retvalue = sin(x)/x |
---|
| 25 | retvalue[x == 0.] = 1. |
---|
| 26 | return retvalue |
---|
| 27 | |
---|
| 28 | def sas_2J1x_x(x): |
---|
| 29 | with np.errstate(all='ignore'): |
---|
| 30 | retvalue = 2*J1(x)/x |
---|
| 31 | retvalue[x == 0] = 1. |
---|
| 32 | return retvalue |
---|
| 33 | |
---|
| 34 | def sas_3j1x_x(x): |
---|
| 35 | """return 3*j1(x)/x""" |
---|
| 36 | retvalue = np.empty_like(x) |
---|
| 37 | with np.errstate(all='ignore'): |
---|
| 38 | # GSL bessel_j1 taylor expansion |
---|
[2cefd79] | 39 | index = (x < 0.25) |
---|
[707cbdb] | 40 | y = x[index]**2 |
---|
| 41 | c1 = -1.0/10.0 |
---|
| 42 | c2 = 1.0/280.0 |
---|
| 43 | c3 = -1.0/15120.0 |
---|
| 44 | c4 = 1.0/1330560.0 |
---|
| 45 | c5 = -1.0/172972800.0 |
---|
| 46 | retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5)))) |
---|
| 47 | index = ~index |
---|
| 48 | y = x[index] |
---|
| 49 | retvalue[index] = 3*(sin(y) - y*cos(y))/y**3 |
---|
| 50 | retvalue[x == 0.] = 1. |
---|
| 51 | return retvalue |
---|
| 52 | |
---|
| 53 | #Used to cross check my models with sasview models |
---|
| 54 | def build_model(model_name, q, **pars): |
---|
| 55 | from sasmodels.core import load_model_info, build_model as build_sasmodel |
---|
| 56 | from sasmodels.data import empty_data1D |
---|
| 57 | from sasmodels.direct_model import DirectModel |
---|
| 58 | model_info = load_model_info(model_name) |
---|
| 59 | model = build_sasmodel(model_info, dtype='double!') |
---|
| 60 | data = empty_data1D(q) |
---|
| 61 | calculator = DirectModel(data, model,cutoff=0) |
---|
| 62 | calculator.pars = pars.copy() |
---|
| 63 | calculator.pars.setdefault('background', 1e-3) |
---|
| 64 | return calculator |
---|
| 65 | |
---|
| 66 | #gives the hardsphere structure factor that sasview uses |
---|
[2cefd79] | 67 | def _hardsphere_simple(q, radius_effective, volfraction): |
---|
| 68 | CUTOFFHS=0.05 |
---|
[707cbdb] | 69 | if fabs(radius_effective) < 1.E-12: |
---|
| 70 | HARDSPH=1.0 |
---|
| 71 | return HARDSPH |
---|
| 72 | X = 1.0/( 1.0 -volfraction) |
---|
| 73 | D= X*X |
---|
| 74 | A= (1.+2.*volfraction)*D |
---|
| 75 | A *=A |
---|
| 76 | X=fabs(q*radius_effective*2.0) |
---|
| 77 | if X < 5.E-06: |
---|
| 78 | HARDSPH=1./A |
---|
| 79 | return HARDSPH |
---|
| 80 | X2 =X*X |
---|
| 81 | B = (1.0 +0.5*volfraction)*D |
---|
| 82 | B *= B |
---|
| 83 | B *= -6.*volfraction |
---|
| 84 | G=0.5*volfraction*A |
---|
| 85 | if X < CUTOFFHS: |
---|
| 86 | 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 |
---|
| 87 | HARDSPH= 1./(1. + volfraction*FF ) |
---|
[2cefd79] | 88 | return HARDSPH |
---|
[707cbdb] | 89 | X4=X2*X2 |
---|
| 90 | S, C = sin(X), cos(X) |
---|
[2cefd79] | 91 | 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 |
---|
| 92 | HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ) |
---|
[707cbdb] | 93 | return HARDSPH |
---|
| 94 | |
---|
[2cefd79] | 95 | def hardsphere_simple(q, radius_effective, volfraction): |
---|
| 96 | SQ = [_hardsphere_simple(qk, radius_effective, volfraction) for qk in q] |
---|
| 97 | return np.array(SQ) |
---|
| 98 | |
---|
[707cbdb] | 99 | #Used in gaussian quadrature for polydispersity |
---|
| 100 | #returns values and the probability of those values based on gaussian distribution |
---|
[2cefd79] | 101 | N_GAUSS = 35 |
---|
| 102 | NSIGMA_GAUSS = 3 |
---|
| 103 | def gaussian_distribution(center, sigma, lb, ub): |
---|
| 104 | #3 standard deviations covers approx. 99.7% |
---|
[707cbdb] | 105 | if sigma != 0: |
---|
[2cefd79] | 106 | nsigmas = NSIGMA_GAUSS |
---|
| 107 | x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS) |
---|
[707cbdb] | 108 | x= x[(x >= lb) & (x <= ub)] |
---|
| 109 | px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) |
---|
| 110 | return x, px |
---|
| 111 | else: |
---|
| 112 | return np.array([center]), np.array([1]) |
---|
| 113 | |
---|
[2cefd79] | 114 | N_SCHULZ = 80 |
---|
| 115 | NSIGMA_SCHULZ = 8 |
---|
[707cbdb] | 116 | def schulz_distribution(center, sigma, lb, ub): |
---|
| 117 | if sigma != 0: |
---|
[2cefd79] | 118 | nsigmas = NSIGMA_SCHULZ |
---|
| 119 | x = np.linspace(center-sigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ) |
---|
[707cbdb] | 120 | x= x[(x >= lb) & (x <= ub)] |
---|
| 121 | R = x/center |
---|
| 122 | z = (center/sigma)**2 |
---|
| 123 | arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) |
---|
| 124 | px = np.exp(arg) |
---|
| 125 | return x, px |
---|
| 126 | else: |
---|
| 127 | return np.array([center]), np.array([1]) |
---|
| 128 | |
---|
| 129 | #returns the effective radius used in sasview |
---|
| 130 | def ER_ellipsoid(radius_polar, radius_equatorial): |
---|
| 131 | ee = np.empty_like(radius_polar) |
---|
| 132 | if radius_polar > radius_equatorial: |
---|
| 133 | ee = (radius_polar**2 - radius_equatorial**2)/radius_polar**2 |
---|
| 134 | elif radius_polar < radius_equatorial: |
---|
| 135 | ee = (radius_equatorial**2 - radius_polar**2) / radius_equatorial**2 |
---|
| 136 | else: |
---|
| 137 | ee = 2*radius_polar |
---|
| 138 | if (radius_polar * radius_equatorial != 0): |
---|
| 139 | bd = 1.0 - ee |
---|
| 140 | e1 = np.sqrt(ee) |
---|
| 141 | b1 = 1.0 + np.arcsin(e1) / (e1*np.sqrt(bd)) |
---|
| 142 | bL = (1.0 + e1) / (1.0 - e1) |
---|
| 143 | b2 = 1.0 + bd / 2 / e1 * np.log(bL) |
---|
| 144 | delta = 0.75 * b1 * b2 |
---|
| 145 | ddd = np.zeros_like(radius_polar) |
---|
| 146 | ddd = 2.0*(delta + 1.0)*radius_polar*radius_equatorial**2 |
---|
| 147 | return 0.5*ddd**(1.0 / 3.0) |
---|
| 148 | |
---|
| 149 | def ellipsoid_volume(radius_polar,radius_equatorial): |
---|
| 150 | volume = (4./3.)*pi*radius_polar*radius_equatorial**2 |
---|
| 151 | return volume |
---|
| 152 | |
---|
| 153 | # F1 is F(q) |
---|
| 154 | # F2 is F(g)^2 |
---|
| 155 | #IQM is I(q) with monodispersity |
---|
| 156 | #IQSM is I(q) with structure factor S(q) and monodispersity |
---|
| 157 | #IQBM is I(q) with Beta Approximation and monodispersity |
---|
| 158 | #SQ is monodisperse approach for structure factor |
---|
| 159 | #SQ_EFF is the effective structure factor from beta approx |
---|
[2cefd79] | 160 | def ellipsoid_theta(q, radius_polar, radius_equatorial, sld, sld_solvent, |
---|
| 161 | volfraction=0, radius_effective=None): |
---|
[707cbdb] | 162 | #creates values z and corresponding probabilities w from legendre-gauss quadrature |
---|
[2cefd79] | 163 | volume = ellipsoid_volume(radius_polar, radius_equatorial) |
---|
[707cbdb] | 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 |
---|
[2cefd79] | 168 | #0->2pi with respect to alpha to -1->1 with respect to z, allowing us to use |
---|
[707cbdb] | 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) |
---|
[0076d6e] | 172 | form = (sld-sld_solvent)*volume*sas_3j1x_x(qk*r) |
---|
| 173 | F2[k] = np.sum(w*form**2) |
---|
| 174 | F1[k] = np.sum(w*form) |
---|
[707cbdb] | 175 | #the 1/2 comes from the change of variables mentioned above |
---|
| 176 | F2 = F2/2.0 |
---|
| 177 | F1 = F1/2.0 |
---|
[2cefd79] | 178 | if radius_effective is None: |
---|
| 179 | radius_effective = ER_ellipsoid(radius_polar,radius_equatorial) |
---|
| 180 | SQ = hardsphere_simple(q, radius_effective, volfraction) |
---|
| 181 | SQ_EFF = 1 + F1**2/F2*(SQ - 1) |
---|
| 182 | IQM = 1e-4*F2/volume |
---|
[707cbdb] | 183 | IQSM = volfraction*IQM*SQ |
---|
| 184 | IQBM = volfraction*IQM*SQ_EFF |
---|
[2cefd79] | 185 | return Theory(Q=q, F1=F1, F2=F2, P=IQM, S=SQ, I=IQSM, Seff=SQ_EFF, Ibeta=IQBM) |
---|
[707cbdb] | 186 | |
---|
[2cefd79] | 187 | #IQD is I(q) polydispursed, IQSD is I(q)S(q) polydispursed, etc. |
---|
[707cbdb] | 188 | #IQBD HAS NOT BEEN CROSS CHECKED AT ALL |
---|
[2cefd79] | 189 | def ellipsoid_pe(q, radius_polar, radius_equatorial, sld, sld_solvent, |
---|
| 190 | radius_polar_pd=0.1, radius_equatorial_pd=0.1, |
---|
| 191 | radius_polar_pd_type='gaussian', |
---|
| 192 | radius_equatorial_pd_type='gaussian', |
---|
| 193 | volfraction=0, radius_effective=None, |
---|
| 194 | background=0, scale=1, |
---|
| 195 | norm='sasview'): |
---|
| 196 | if norm not in ['sasview', 'sasfit', 'yun']: |
---|
| 197 | raise TypeError("unknown norm "+norm) |
---|
| 198 | if radius_polar_pd_type == 'gaussian': |
---|
| 199 | Rp_val, Rp_prob = gaussian_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) |
---|
| 200 | elif radius_polar_pd_type == 'schulz': |
---|
| 201 | Rp_val, Rp_prob = schulz_distribution(radius_polar, radius_polar_pd*radius_polar, 0, inf) |
---|
| 202 | if radius_equatorial_pd_type == 'gaussian': |
---|
| 203 | Re_val, Re_prob = gaussian_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) |
---|
| 204 | elif radius_equatorial_pd_type == 'schulz': |
---|
| 205 | Re_val, Re_prob = schulz_distribution(radius_equatorial, radius_equatorial_pd*radius_equatorial, 0, inf) |
---|
[0076d6e] | 206 | total_weight = total_volume = 0 |
---|
| 207 | radius_eff = 0 |
---|
| 208 | F1, F2 = np.zeros_like(q), np.zeros_like(q) |
---|
[707cbdb] | 209 | for k, Rpk in enumerate(Rp_val): |
---|
| 210 | for i, Rei in enumerate(Re_val): |
---|
[2cefd79] | 211 | theory = ellipsoid_theta(q,Rpk,Rei,sld,sld_solvent) |
---|
| 212 | volume = ellipsoid_volume(Rpk, Rei) |
---|
[0076d6e] | 213 | weight = Rp_prob[k]*Re_prob[i] |
---|
| 214 | total_weight += weight |
---|
| 215 | total_volume += weight*volume |
---|
| 216 | F1 += theory.F1*weight |
---|
| 217 | F2 += theory.F2*weight |
---|
| 218 | radius_eff += weight*ER_ellipsoid(Rpk,Rei) |
---|
| 219 | F1 /= total_weight |
---|
| 220 | F2 /= total_weight |
---|
| 221 | average_volume = total_volume/total_weight |
---|
[2cefd79] | 222 | if radius_effective is None: |
---|
| 223 | radius_effective = radius_eff/total_weight |
---|
| 224 | if norm == 'sasfit': |
---|
| 225 | IQD = F2 |
---|
| 226 | elif norm == 'sasview': |
---|
[0076d6e] | 227 | # Note: internally, sasview uses F2/total_volume because: |
---|
| 228 | # average_volume = total_volume/total_weight |
---|
| 229 | # F2/total_weight / average_volume |
---|
| 230 | # = F2/total_weight / total_volume/total_weight |
---|
| 231 | # = F2/total_volume |
---|
| 232 | IQD = F2/average_volume*1e-4*volfraction |
---|
[2cefd79] | 233 | elif norm == 'yun': |
---|
[0076d6e] | 234 | F1 *= 1e-6 # Yun is using sld in 1/A^2, not 1e-6/A^2 |
---|
| 235 | F2 *= 1e-12 |
---|
| 236 | IQD = F2/average_volume*1e8*volfraction |
---|
[e262dd6] | 237 | SQ = hardsphere_simple(q, radius_effective, volfraction) |
---|
| 238 | beta = F1**2/F2 |
---|
| 239 | SQ_EFF = 1 + beta*(SQ - 1) |
---|
| 240 | IQSD = IQD*SQ |
---|
| 241 | IQBD = IQD*SQ_EFF |
---|
[2cefd79] | 242 | return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) |
---|
[707cbdb] | 243 | |
---|
| 244 | #polydispersity for sphere |
---|
[2cefd79] | 245 | def sphere_r(q,radius,sld,sld_solvent, |
---|
| 246 | radius_pd=0.1, radius_pd_type='gaussian', |
---|
| 247 | volfraction=0, radius_effective=None, |
---|
| 248 | background=0, scale=1, |
---|
| 249 | norm='sasview'): |
---|
| 250 | if norm not in ['sasview', 'sasfit', 'yun']: |
---|
| 251 | raise TypeError("unknown norm "+norm) |
---|
| 252 | if radius_pd_type == 'gaussian': |
---|
[707cbdb] | 253 | radius_val, radius_prob = gaussian_distribution(radius, radius_pd*radius, 0, inf) |
---|
[2cefd79] | 254 | elif radius_pd_type == 'schulz': |
---|
[707cbdb] | 255 | radius_val, radius_prob = schulz_distribution(radius, radius_pd*radius, 0, inf) |
---|
[0076d6e] | 256 | total_weight = total_volume = 0 |
---|
[707cbdb] | 257 | F1 = np.zeros_like(q) |
---|
[2cefd79] | 258 | F2 = np.zeros_like(q) |
---|
| 259 | for k, rk in enumerate(radius_val): |
---|
| 260 | volume = 4./3.*pi*rk**3 |
---|
[0076d6e] | 261 | total_weight += radius_prob[k] |
---|
| 262 | total_volume += radius_prob[k]*volume |
---|
| 263 | form = (sld-sld_solvent)*volume*sas_3j1x_x(q*rk) |
---|
| 264 | F2 += radius_prob[k]*form**2 |
---|
| 265 | F1 += radius_prob[k]*form |
---|
| 266 | F1 /= total_weight |
---|
| 267 | F2 /= total_weight |
---|
| 268 | average_volume = total_volume/total_weight |
---|
| 269 | |
---|
[2cefd79] | 270 | if radius_effective is None: |
---|
| 271 | radius_effective = radius |
---|
[0076d6e] | 272 | average_volume = total_volume/total_weight |
---|
[2cefd79] | 273 | if norm == 'sasfit': |
---|
| 274 | IQD = F2 |
---|
| 275 | elif norm == 'sasview': |
---|
[0076d6e] | 276 | IQD = F2/average_volume*1e-4*volfraction |
---|
[2cefd79] | 277 | elif norm == 'yun': |
---|
[0076d6e] | 278 | F1 *= 1e-6 # Yun is using sld in 1/A^2, not 1e-6/A^2 |
---|
| 279 | F2 *= 1e-12 |
---|
| 280 | IQD = F2/average_volume*1e8*volfraction |
---|
[e262dd6] | 281 | SQ = hardsphere_simple(q, radius_effective, volfraction) |
---|
| 282 | beta = F1**2/F2 |
---|
| 283 | SQ_EFF = 1 + beta*(SQ - 1) |
---|
| 284 | IQSD = IQD*SQ |
---|
| 285 | IQBD = IQD*SQ_EFF |
---|
[2cefd79] | 286 | return Theory(Q=q, F1=F1, F2=F2, P=IQD, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) |
---|
[707cbdb] | 287 | |
---|
| 288 | ############################################################################### |
---|
| 289 | ############################################################################### |
---|
| 290 | ############################################################################### |
---|
| 291 | ################## ################## |
---|
| 292 | ################## TESTS ################## |
---|
| 293 | ################## ################## |
---|
| 294 | ############################################################################### |
---|
| 295 | ############################################################################### |
---|
| 296 | ############################################################################### |
---|
| 297 | |
---|
[2cefd79] | 298 | def popn(d, keys): |
---|
| 299 | """ |
---|
| 300 | Splits a dict into two, with any key of *d* which is in *keys* removed |
---|
| 301 | from *d* and added to *b*. Returns *b*. |
---|
| 302 | """ |
---|
| 303 | b = {} |
---|
| 304 | for k in keys: |
---|
| 305 | try: |
---|
| 306 | b[k] = d.pop(k) |
---|
| 307 | except KeyError: |
---|
| 308 | pass |
---|
| 309 | return b |
---|
[707cbdb] | 310 | |
---|
[2cefd79] | 311 | def sasmodels_theory(q, Pname, **pars): |
---|
| 312 | """ |
---|
| 313 | Call sasmodels to compute the model with and without a hard sphere |
---|
| 314 | structure factor. |
---|
| 315 | """ |
---|
| 316 | #mono_pars = {k: (0 if k.endswith('_pd') else v) for k, v in pars.items()} |
---|
| 317 | Ppars = pars.copy() |
---|
| 318 | Spars = popn(Ppars, ['radius_effective', 'volfraction']) |
---|
| 319 | Ipars = pars.copy() |
---|
| 320 | |
---|
| 321 | # Autofill npts and nsigmas for the given pd type |
---|
| 322 | for k, v in pars.items(): |
---|
| 323 | if k.endswith("_pd_type"): |
---|
| 324 | if v == "gaussian": |
---|
| 325 | n, nsigmas = N_GAUSS, NSIGMA_GAUSS |
---|
| 326 | elif v == "schulz": |
---|
| 327 | n, nsigmas = N_SCHULZ, NSIGMA_SCHULZ |
---|
| 328 | Ppars.setdefault(k.replace("_pd_type", "_pd_n"), n) |
---|
| 329 | Ppars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) |
---|
| 330 | Ipars.setdefault(k.replace("_pd_type", "_pd_n"), n) |
---|
| 331 | Ipars.setdefault(k.replace("_pd_type", "_pd_nsigma"), nsigmas) |
---|
| 332 | |
---|
| 333 | #Ppars['scale'] = Spars.get('volfraction', 1) |
---|
| 334 | P = build_model(Pname, q) |
---|
| 335 | S = build_model("hardsphere", q) |
---|
| 336 | I = build_model(Pname+"@hardsphere", q) |
---|
| 337 | Pq = P(**Ppars)*pars.get('volfraction', 1) |
---|
| 338 | #Sq = S(**Spars) |
---|
| 339 | Iq = I(**Ipars) |
---|
| 340 | #Iq = Pq*Sq*pars.get('volfraction', 1) |
---|
| 341 | Sq = Iq/Pq |
---|
| 342 | return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None) |
---|
| 343 | |
---|
| 344 | def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): |
---|
| 345 | """ |
---|
| 346 | Plot fields in common between target and actual, along with relative error. |
---|
| 347 | """ |
---|
| 348 | available = [s for s in fields.split() |
---|
| 349 | if getattr(target, s) is not None and getattr(actual, s) is not None] |
---|
| 350 | rows = len(available) |
---|
| 351 | for row, field in enumerate(available): |
---|
| 352 | Q = target.Q |
---|
| 353 | I1, I2 = getattr(target, field), getattr(actual, field) |
---|
| 354 | plt.subplot(rows, 2, 2*row+1) |
---|
| 355 | plt.loglog(Q, abs(I1), label="target "+field) |
---|
| 356 | plt.loglog(Q, abs(I2), label="value "+field) |
---|
| 357 | #plt.legend(loc="upper left", bbox_to_anchor=(1,1)) |
---|
| 358 | plt.legend(loc='lower left') |
---|
| 359 | plt.subplot(rows, 2, 2*row+2) |
---|
[0076d6e] | 360 | plt.semilogx(Q, I2/I1 - 1, label="relative error") |
---|
| 361 | #plt.semilogx(Q, I1/I2 - 1, label="relative error") |
---|
[707cbdb] | 362 | plt.tight_layout() |
---|
[2cefd79] | 363 | plt.suptitle(title) |
---|
[707cbdb] | 364 | plt.show() |
---|
| 365 | |
---|
[2cefd79] | 366 | def data_file(name): |
---|
| 367 | return os.path.join(BETA_DIR, 'data_files', name) |
---|
| 368 | |
---|
| 369 | def load_sasfit(path): |
---|
| 370 | data = np.loadtxt(path, dtype=str, delimiter=';').T |
---|
| 371 | data = np.vstack((map(float, v) for v in data[0:2])) |
---|
| 372 | return data |
---|
| 373 | |
---|
| 374 | COMPARISON = {} # Type: Dict[(str,str,str)] -> Callable[(), None] |
---|
| 375 | |
---|
| 376 | def compare_sasview_sphere(pd_type='schulz'): |
---|
| 377 | q = np.logspace(-5, 0, 250) |
---|
| 378 | model = 'sphere' |
---|
| 379 | pars = dict( |
---|
| 380 | radius=20,sld=4,sld_solvent=1, |
---|
| 381 | background=0, |
---|
| 382 | radius_pd=.1, radius_pd_type=pd_type, |
---|
| 383 | volfraction=0.15, |
---|
| 384 | #radius_effective=12.59921049894873, # equivalent average sphere radius |
---|
| 385 | ) |
---|
| 386 | target = sasmodels_theory(q, model, **pars) |
---|
| 387 | actual = sphere_r(q, norm='sasview', **pars) |
---|
| 388 | title = " ".join(("sasmodels", model, pd_type)) |
---|
| 389 | compare(title, target, actual) |
---|
| 390 | COMPARISON[('sasview','sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') |
---|
| 391 | COMPARISON[('sasview','sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') |
---|
| 392 | |
---|
| 393 | def compare_sasview_ellipsoid(pd_type='gaussian'): |
---|
| 394 | q = np.logspace(-5, 0, 50) |
---|
| 395 | model = 'ellipsoid' |
---|
| 396 | pars = dict( |
---|
| 397 | radius_polar=20,radius_equatorial=400,sld=4,sld_solvent=1, |
---|
| 398 | background=0, |
---|
| 399 | radius_polar_pd=.1, radius_polar_pd_type=pd_type, |
---|
| 400 | radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, |
---|
| 401 | volfraction=0.15, |
---|
| 402 | #radius_effective=12.59921049894873, |
---|
| 403 | ) |
---|
| 404 | target = sasmodels_theory(q, model, **pars) |
---|
| 405 | actual = ellipsoid_pe(q, norm='sasview', **pars) |
---|
| 406 | title = " ".join(("sasmodels", model, pd_type)) |
---|
| 407 | compare(title, target, actual) |
---|
[0076d6e] | 408 | COMPARISON[('sasview','ellipsoid','gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian') |
---|
| 409 | COMPARISON[('sasview','ellipsoid','schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz') |
---|
[2cefd79] | 410 | |
---|
| 411 | def compare_yun_ellipsoid_mono(): |
---|
| 412 | pars = { |
---|
| 413 | 'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', |
---|
| 414 | 'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', |
---|
| 415 | 'sld': 2, 'sld_solvent': 1, |
---|
| 416 | 'volfraction': 0.15, |
---|
| 417 | # Yun uses radius for same volume sphere for effective radius |
---|
| 418 | # whereas sasview uses the average curvature. |
---|
| 419 | 'radius_effective': 12.59921049894873, |
---|
| 420 | } |
---|
| 421 | |
---|
| 422 | data = np.loadtxt(data_file('yun_ellipsoid.dat'),skiprows=2).T |
---|
| 423 | Q = data[0] |
---|
| 424 | F1 = data[1] |
---|
[0076d6e] | 425 | P = data[3]*pars['volfraction'] |
---|
[2cefd79] | 426 | S = data[5] |
---|
| 427 | Seff = data[6] |
---|
[0076d6e] | 428 | target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) |
---|
[2cefd79] | 429 | actual = ellipsoid_pe(Q, norm='yun', **pars) |
---|
| 430 | title = " ".join(("yun", "ellipsoid", "no pd")) |
---|
| 431 | #compare(title, target, actual, fields="P S I Seff Ibeta") |
---|
| 432 | compare(title, target, actual) |
---|
| 433 | COMPARISON[('yun','ellipsoid','gaussian')] = compare_yun_ellipsoid_mono |
---|
| 434 | COMPARISON[('yun','ellipsoid','schulz')] = compare_yun_ellipsoid_mono |
---|
| 435 | |
---|
[0076d6e] | 436 | def compare_yun_sphere_gauss(): |
---|
[e262dd6] | 437 | # Note: yun uses gauss limits from R0/10 to R0 + 5 sigma steps sigma/100 |
---|
| 438 | # With pd = 0.1, that's 14 sigma and 1400 points. |
---|
[0076d6e] | 439 | pars = { |
---|
| 440 | 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', |
---|
| 441 | 'sld': 6, 'sld_solvent': 0, |
---|
| 442 | 'volfraction': 0.1, |
---|
| 443 | } |
---|
| 444 | |
---|
| 445 | data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'),skiprows=2).T |
---|
| 446 | Q = data[0] |
---|
| 447 | F1 = data[1] |
---|
[cdd676e] | 448 | F2 = data[2] |
---|
[0076d6e] | 449 | P = data[3] |
---|
| 450 | S = data[5] |
---|
| 451 | Seff = data[6] |
---|
[cdd676e] | 452 | target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) |
---|
[0076d6e] | 453 | actual = sphere_r(Q, norm='yun', **pars) |
---|
| 454 | title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) |
---|
| 455 | compare(title, target, actual) |
---|
| 456 | data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'),skiprows=2).T |
---|
| 457 | pars.update(radius_pd=0.15) |
---|
| 458 | Q = data[0] |
---|
| 459 | F1 = data[1] |
---|
[cdd676e] | 460 | F2 = data[2] |
---|
[0076d6e] | 461 | P = data[3] |
---|
| 462 | S = data[5] |
---|
| 463 | Seff = data[6] |
---|
[cdd676e] | 464 | target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) |
---|
[0076d6e] | 465 | actual = sphere_r(Q, norm='yun', **pars) |
---|
| 466 | title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) |
---|
| 467 | compare(title, target, actual) |
---|
| 468 | COMPARISON[('yun','sphere','gaussian')] = compare_yun_sphere_gauss |
---|
| 469 | |
---|
| 470 | |
---|
[2cefd79] | 471 | def compare_sasfit_sphere_gauss(): |
---|
| 472 | #N=1,s=2,X0=20,distr radius R=20,eta_core=4,eta_solv=1,.3 |
---|
| 473 | pars = { |
---|
| 474 | 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'gaussian', |
---|
| 475 | 'sld': 4, 'sld_solvent': 1, |
---|
| 476 | 'volfraction': 0.3, |
---|
| 477 | } |
---|
[0076d6e] | 478 | |
---|
[2cefd79] | 479 | Q, IQ = load_sasfit(data_file('sasfit_sphere_IQD.txt')) |
---|
| 480 | Q, IQSD = load_sasfit(data_file('sasfit_sphere_IQSD.txt')) |
---|
| 481 | Q, IQBD = load_sasfit(data_file('sasfit_sphere_IQBD.txt')) |
---|
| 482 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_sphere_sq.txt')) |
---|
| 483 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_sphere_sqeff.txt')) |
---|
| 484 | target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=SQ, I=IQSD, Seff=SQ_EFF, Ibeta=IQBD) |
---|
| 485 | actual = sphere_r(Q, norm="sasfit", **pars) |
---|
| 486 | title = " ".join(("sasfit", "sphere", "pd=10% gaussian")) |
---|
| 487 | compare(title, target, actual) |
---|
| 488 | #compare(title, target, actual, fields="P") |
---|
| 489 | COMPARISON[('sasfit','sphere','gaussian')] = compare_sasfit_sphere_gauss |
---|
| 490 | |
---|
| 491 | def compare_sasfit_sphere_schulz(): |
---|
[707cbdb] | 492 | #radius=20,sld=4,sld_solvent=1,volfraction=0.3,radius_pd=0.1 |
---|
| 493 | #We have scaled the output from sasfit by 1e-4*volume*volfraction |
---|
| 494 | #0.10050378152592121 |
---|
[2cefd79] | 495 | pars = { |
---|
| 496 | 'radius': 20, 'radius_pd': 0.1, 'radius_pd_type': 'schulz', |
---|
| 497 | 'sld': 4, 'sld_solvent': 1, |
---|
| 498 | 'volfraction': 0.3, |
---|
| 499 | } |
---|
| 500 | |
---|
| 501 | Q, IQ = load_sasfit(data_file('richard_test.txt')) |
---|
| 502 | Q, IQSD = load_sasfit(data_file('richard_test2.txt')) |
---|
| 503 | Q, IQBD = load_sasfit(data_file('richard_test3.txt')) |
---|
| 504 | target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) |
---|
| 505 | actual = sphere_r(Q, norm="sasfit", **pars) |
---|
| 506 | title = " ".join(("sasfit", "sphere", "pd=10% schulz")) |
---|
| 507 | compare(title, target, actual) |
---|
| 508 | COMPARISON[('sasfit','sphere','schulz')] = compare_sasfit_sphere_schulz |
---|
[707cbdb] | 509 | |
---|
[2cefd79] | 510 | def compare_sasfit_ellipsoid_schulz(): |
---|
[707cbdb] | 511 | #polarradius=20, equatorialradius=10, sld=4,sld_solvent=1,volfraction=0.3,radius_polar_pd=0.1 |
---|
[2cefd79] | 512 | #Effective radius =13.1353356684 |
---|
| 513 | #We have scaled the output from sasfit by 1e-4*volume*volfraction |
---|
| 514 | #0.10050378152592121 |
---|
| 515 | pars = { |
---|
| 516 | 'radius_polar': 20, 'radius_polar_pd': 0.1, 'radius_polar_pd_type': 'schulz', |
---|
| 517 | 'radius_equatorial': 10, 'radius_equatorial_pd': 0., 'radius_equatorial_pd_type': 'schulz', |
---|
| 518 | 'sld': 4, 'sld_solvent': 1, |
---|
| 519 | 'volfraction': 0.3, 'radius_effective': 13.1353356684, |
---|
| 520 | } |
---|
[0076d6e] | 521 | |
---|
[2cefd79] | 522 | Q, IQ = load_sasfit(data_file('richard_test4.txt')) |
---|
| 523 | Q, IQSD = load_sasfit(data_file('richard_test5.txt')) |
---|
| 524 | Q, IQBD = load_sasfit(data_file('richard_test6.txt')) |
---|
| 525 | target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) |
---|
| 526 | actual = ellipsoid_pe(Q, norm="sasfit", **pars) |
---|
| 527 | title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz")) |
---|
| 528 | compare(title, target, actual) |
---|
| 529 | COMPARISON[('sasfit','ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz |
---|
| 530 | |
---|
| 531 | |
---|
| 532 | def compare_sasfit_ellipsoid_gaussian(): |
---|
| 533 | pars = { |
---|
| 534 | 'radius_polar': 20, 'radius_polar_pd': 0, 'radius_polar_pd_type': 'gaussian', |
---|
| 535 | 'radius_equatorial': 10, 'radius_equatorial_pd': 0, 'radius_equatorial_pd_type': 'gaussian', |
---|
| 536 | 'sld': 4, 'sld_solvent': 1, |
---|
| 537 | 'volfraction': 0, 'radius_effective': None, |
---|
| 538 | } |
---|
| 539 | |
---|
| 540 | #Rp=20,Re=10,eta_core=4,eta_solv=1 |
---|
| 541 | Q, PQ0 = load_sasfit(data_file('sasfit_ellipsoid_IQM.txt')) |
---|
| 542 | pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0, radius_effective=None) |
---|
| 543 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 544 | target = Theory(Q=Q, P=PQ0) |
---|
| 545 | compare("sasfit ellipsoid no poly", target, actual); plt.show() |
---|
| 546 | |
---|
| 547 | #N=1,s=2,X0=20,distr 10% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 548 | Q, PQ_Rp10 = load_sasfit(data_file('sasfit_ellipsoid_IQD.txt')) |
---|
| 549 | pars.update(volfraction=0, radius_polar_pd=0.1, radius_equatorial_pd=0.0, radius_effective=None) |
---|
| 550 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 551 | target = Theory(Q=Q, P=PQ_Rp10) |
---|
| 552 | compare("sasfit ellipsoid P(Q) 10% Rp", target, actual); plt.show() |
---|
| 553 | #N=1,s=1,X0=10,distr 10% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 554 | Q, PQ_Re10 = load_sasfit(data_file('sasfit_ellipsoid_IQD2.txt')) |
---|
| 555 | pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.1, radius_effective=None) |
---|
| 556 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 557 | target = Theory(Q=Q, P=PQ_Re10) |
---|
| 558 | compare("sasfit ellipsoid P(Q) 10% Re", target, actual); plt.show() |
---|
| 559 | #N=1,s=6,X0=20,distr 30% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 560 | Q, PQ_Rp30 = load_sasfit(data_file('sasfit_ellipsoid_IQD3.txt')) |
---|
| 561 | pars.update(volfraction=0, radius_polar_pd=0.3, radius_equatorial_pd=0.0, radius_effective=None) |
---|
| 562 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 563 | target = Theory(Q=Q, P=PQ_Rp30) |
---|
| 564 | compare("sasfit ellipsoid P(Q) 30% Rp", target, actual); plt.show() |
---|
| 565 | #N=1,s=3,X0=10,distr 30% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 566 | Q, PQ_Re30 = load_sasfit(data_file('sasfit_ellipsoid_IQD4.txt')) |
---|
| 567 | pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.3, radius_effective=None) |
---|
| 568 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 569 | target = Theory(Q=Q, P=PQ_Re30) |
---|
| 570 | compare("sasfit ellipsoid P(Q) 30% Re", target, actual); plt.show() |
---|
| 571 | #N=1,s=12,X0=20,distr 60% polar Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 572 | Q, PQ_Rp60 = load_sasfit(data_file('sasfit_ellipsoid_IQD5.txt')) |
---|
| 573 | pars.update(volfraction=0, radius_polar_pd=0.6, radius_equatorial_pd=0.0, radius_effective=None) |
---|
| 574 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 575 | target = Theory(Q=Q, P=PQ_Rp60) |
---|
| 576 | compare("sasfit ellipsoid P(Q) 60% Rp", target, actual); plt.show() |
---|
| 577 | #N=1,s=6,X0=10,distr 60% equatorial Rp=20,Re=10,eta_core=4,eta_solv=1, no structure poly |
---|
| 578 | Q, PQ_Re60 = load_sasfit(data_file('sasfit_ellipsoid_IQD6.txt')) |
---|
| 579 | pars.update(volfraction=0, radius_polar_pd=0.0, radius_equatorial_pd=0.6, radius_effective=None) |
---|
| 580 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 581 | target = Theory(Q=Q, P=PQ_Re60) |
---|
| 582 | compare("sasfit ellipsoid P(Q) 60% Re", target, actual); plt.show() |
---|
| 583 | |
---|
| 584 | #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.15 |
---|
| 585 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq.txt')) |
---|
| 586 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff.txt')) |
---|
| 587 | pars.update(volfraction=0.15, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) |
---|
| 588 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 589 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 590 | compare("sasfit ellipsoid P(Q) 10% Rp 15% Vf", target, actual); plt.show() |
---|
| 591 | #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.15 |
---|
| 592 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq2.txt')) |
---|
| 593 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff2.txt')) |
---|
| 594 | pars.update(volfraction=0.15, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) |
---|
| 595 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 596 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 597 | compare("sasfit ellipsoid P(Q) 30% Rp 15% Vf", target, actual); plt.show() |
---|
| 598 | #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.15 |
---|
| 599 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq3.txt')) |
---|
| 600 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff3.txt')) |
---|
| 601 | pars.update(volfraction=0.15, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) |
---|
| 602 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 603 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 604 | compare("sasfit ellipsoid P(Q) 60% Rp 15% Vf", target, actual); plt.show() |
---|
| 605 | |
---|
| 606 | #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.3 |
---|
| 607 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq4.txt')) |
---|
| 608 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff4.txt')) |
---|
| 609 | pars.update(volfraction=0.3, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) |
---|
| 610 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 611 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 612 | compare("sasfit ellipsoid P(Q) 10% Rp 30% Vf", target, actual); plt.show() |
---|
| 613 | #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.3 |
---|
| 614 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq5.txt')) |
---|
| 615 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff5.txt')) |
---|
| 616 | pars.update(volfraction=0.3, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) |
---|
| 617 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 618 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 619 | compare("sasfit ellipsoid P(Q) 30% Rp 30% Vf", target, actual); plt.show() |
---|
| 620 | #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.3 |
---|
| 621 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq6.txt')) |
---|
| 622 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff6.txt')) |
---|
| 623 | pars.update(volfraction=0.3, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) |
---|
| 624 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 625 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 626 | compare("sasfit ellipsoid P(Q) 60% Rp 30% Vf", target, actual); plt.show() |
---|
| 627 | |
---|
| 628 | #N=1,s=2,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.1354236254,.6 |
---|
| 629 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq7.txt')) |
---|
| 630 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff7.txt')) |
---|
| 631 | pars.update(volfraction=0.6, radius_polar_pd=0.1, radius_equatorial_pd=0, radius_effective=13.1354236254) |
---|
| 632 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 633 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 634 | compare("sasfit ellipsoid P(Q) 10% Rp 60% Vf", target, actual); plt.show() |
---|
| 635 | #N=1,s=6,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.0901197149,.6 |
---|
| 636 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq8.txt')) |
---|
| 637 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff8.txt')) |
---|
| 638 | pars.update(volfraction=0.6, radius_polar_pd=0.3, radius_equatorial_pd=0, radius_effective=13.0901197149) |
---|
| 639 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 640 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 641 | compare("sasfit ellipsoid P(Q) 30% Rp 60% Vf", target, actual); plt.show() |
---|
| 642 | #N=1,s=12,X0=20,distr polar Rp=20,Re=10,eta_core=4,eta_solv=1, hardsphere ,13.336060917,.6 |
---|
| 643 | Q, SQ = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sq9.txt')) |
---|
| 644 | Q, SQ_EFF = load_sasfit(data_file('sasfit_polydisperse_ellipsoid_sqeff9.txt')) |
---|
| 645 | pars.update(volfraction=0.6, radius_polar_pd=0.6, radius_equatorial_pd=0, radius_effective=13.336060917) |
---|
| 646 | actual = ellipsoid_pe(Q, norm='sasfit', **pars) |
---|
| 647 | target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) |
---|
| 648 | compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show() |
---|
[757a3ff] | 649 | COMPARISON[('sasfit','ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian |
---|
[2cefd79] | 650 | |
---|
| 651 | def main(): |
---|
| 652 | key = tuple(sys.argv[1:]) |
---|
| 653 | if key not in COMPARISON: |
---|
| 654 | print("usage: sasfit_compare.py [sasview|sasfit|yun] [sphere|ellipsoid] [gaussian|schulz]") |
---|
| 655 | return |
---|
| 656 | comparison = COMPARISON[key] |
---|
| 657 | comparison() |
---|
| 658 | |
---|
| 659 | if __name__ == "__main__": |
---|
| 660 | main() |
---|