Changeset 7b0abf8 in sasmodels
 Timestamp:
 Aug 7, 2018 9:07:11 PM (4 months ago)
 Branches:
 F1F2models_grethe, beta_approx, py3, ticket1015gpumemerror, ticket1157, ticket608userdefinedweights, ticket_1156
 Children:
 c036ddb
 Parents:
 502c7b8
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/beta/sasfit_compare.py
r502c7b8 r7b0abf8 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)) 6 SASMODELS_DIR = r"C:\Source\sasmodels" 5 SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 7 6 sys.path.insert(0, SASMODELS_DIR) 8 7 … … 41 40 y = x[index]**2 42 41 c1 = 1.0/10.0 43 c2 = 42 c2 = +1.0/280.0 44 43 c3 = 1.0/15120.0 45 c4 = 44 c4 = +1.0/1330560.0 46 45 c5 = 1.0/172972800.0 47 46 retvalue[index] = 1.0 + y*(c1 + y*(c2 + y*(c3 + y*(c4 + y*c5)))) … … 67 66 #gives the hardsphere structure factor that sasview uses 68 67 def _hardsphere_simple(q, radius_effective, volfraction): 69 CUTOFFHS =0.0568 CUTOFFHS = 0.05 70 69 if fabs(radius_effective) < 1.E12: 71 HARDSPH =1.070 HARDSPH = 1.0 72 71 return HARDSPH 73 X = 1.0/( 74 D = X*X75 A = (1.+2.*volfraction)*D76 A *= A77 X =fabs(q*radius_effective*2.0)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) 78 77 if X < 5.E06: 79 HARDSPH =1./A78 HARDSPH = 1./A 80 79 return HARDSPH 81 X2 = X*X80 X2 = X*X 82 81 B = (1.0 +0.5*volfraction)*D 83 82 B *= B 84 83 B *= 6.*volfraction 85 G =0.5*volfraction*A84 G = 0.5*volfraction*A 86 85 if X < CUTOFFHS: 87 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 88 HARDSPH = 1./(1. + volfraction*FF )87 HARDSPH = 1./(1. + volfraction*FF ) 89 88 return HARDSPH 90 X4 =X2*X289 X4 = X2*X2 91 90 S, C = sin(X), cos(X) 92 FF = ((G*( (4.*X2 24.)*X*S (X4 12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S (X22.)*C 2.) )/X + A*(SX*C))/X93 HARDSPH = 1./(1. + 24.*volfraction*FF/X2)91 FF = ((G*( (4.*X2 24.)*X*S (X4 12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S (X22.)*C 2.) )/X + A*(SX*C))/X 92 HARDSPH = 1./(1. + 24.*volfraction*FF/X2) 94 93 return HARDSPH 95 94 … … 107 106 nsigmas = NSIGMA_GAUSS 108 107 x = np.linspace(centersigma*nsigmas, center+sigma*nsigmas, num=N_GAUSS) 109 x = x[(x >= lb) & (x <= ub)]108 x = x[(x >= lb) & (x <= ub)] 110 109 px = np.exp((xcenter)**2 / (2.0 * sigma * sigma)) 111 110 return x, px … … 119 118 nsigmas = NSIGMA_SCHULZ 120 119 x = np.linspace(centersigma*nsigmas, center+sigma*nsigmas, num=N_SCHULZ) 121 x = x[(x >= lb) & (x <= ub)]120 x = x[(x >= lb) & (x <= ub)] 122 121 R = x/center 123 122 z = (center/sigma)**2 … … 137 136 else: 138 137 ee = 2*radius_polar 139 if (radius_polar * radius_equatorial != 0):138 if radius_polar * radius_equatorial != 0: 140 139 bd = 1.0  ee 141 140 e1 = np.sqrt(ee) … … 148 147 return 0.5*ddd**(1.0 / 3.0) 149 148 150 def ellipsoid_volume(radius_polar, radius_equatorial):149 def ellipsoid_volume(radius_polar, radius_equatorial): 151 150 volume = (4./3.)*pi*radius_polar*radius_equatorial**2 152 151 return volume … … 210 209 for k, Rpk in enumerate(Rp_val): 211 210 for i, Rei in enumerate(Re_val): 212 theory = ellipsoid_theta(q, Rpk,Rei,sld,sld_solvent)211 theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent) 213 212 volume = ellipsoid_volume(Rpk, Rei) 214 213 weight = Rp_prob[k]*Re_prob[i] … … 217 216 F1 += theory.F1*weight 218 217 F2 += theory.F2*weight 219 radius_eff += weight*ER_ellipsoid(Rpk, Rei)218 radius_eff += weight*ER_ellipsoid(Rpk, Rei) 220 219 F1 /= total_weight 221 220 F2 /= total_weight … … 344 343 #Sq = Iq/Pq 345 344 #Iq = None#= Sq = None 346 r =I._kernel.results345 r = I._kernel.results 347 346 return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 348 347 … … 383 382 model = 'sphere' 384 383 pars = dict( 385 radius=20, sld=4,sld_solvent=1,384 radius=20, sld=4, sld_solvent=1, 386 385 background=0, 387 386 radius_pd=.1, radius_pd_type=pd_type, … … 393 392 title = " ".join(("sasmodels", model, pd_type)) 394 393 compare(title, target, actual) 395 COMPARISON[('sasview', 'sphere','gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian')396 COMPARISON[('sasview', 'sphere','schulz')] = lambda: compare_sasview_sphere(pd_type='schulz')394 COMPARISON[('sasview', 'sphere', 'gaussian')] = lambda: compare_sasview_sphere(pd_type='gaussian') 395 COMPARISON[('sasview', 'sphere', 'schulz')] = lambda: compare_sasview_sphere(pd_type='schulz') 397 396 398 397 def compare_sasview_ellipsoid(pd_type='gaussian'): … … 400 399 model = 'ellipsoid' 401 400 pars = dict( 402 radius_polar=20, radius_equatorial=400,sld=4,sld_solvent=1,401 radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, 403 402 background=0, 404 403 radius_polar_pd=.1, radius_polar_pd_type=pd_type, … … 412 411 title = " ".join(("sasmodels", model, pd_type)) 413 412 compare(title, target, actual) 414 COMPARISON[('sasview', 'ellipsoid','gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian')415 COMPARISON[('sasview', 'ellipsoid','schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz')413 COMPARISON[('sasview', 'ellipsoid', 'gaussian')] = lambda: compare_sasview_ellipsoid(pd_type='gaussian') 414 COMPARISON[('sasview', 'ellipsoid', 'schulz')] = lambda: compare_sasview_ellipsoid(pd_type='schulz') 416 415 417 416 def compare_yun_ellipsoid_mono(): … … 437 436 #compare(title, target, actual, fields="P S I Seff Ibeta") 438 437 compare(title, target, actual) 439 COMPARISON[('yun', 'ellipsoid','gaussian')] = compare_yun_ellipsoid_mono440 COMPARISON[('yun', 'ellipsoid','schulz')] = compare_yun_ellipsoid_mono438 COMPARISON[('yun', 'ellipsoid', 'gaussian')] = compare_yun_ellipsoid_mono 439 COMPARISON[('yun', 'ellipsoid', 'schulz')] = compare_yun_ellipsoid_mono 441 440 442 441 def compare_yun_sphere_gauss(): … … 449 448 } 450 449 451 data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T450 data = np.loadtxt(data_file('testPolydisperseGaussianSphere.dat'), skiprows=2).T 452 451 Q = data[0] 453 452 F1 = data[1] … … 460 459 title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) 461 460 compare(title, target, actual) 462 data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T461 data = np.loadtxt(data_file('testPolydisperseGaussianSphere2.dat'), skiprows=2).T 463 462 pars.update(radius_pd=0.15) 464 463 Q = data[0] … … 472 471 title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) 473 472 compare(title, target, actual) 474 COMPARISON[('yun', 'sphere','gaussian')] = compare_yun_sphere_gauss473 COMPARISON[('yun', 'sphere', 'gaussian')] = compare_yun_sphere_gauss 475 474 476 475 … … 493 492 compare(title, target, actual) 494 493 #compare(title, target, actual, fields="P") 495 COMPARISON[('sasfit', 'sphere','gaussian')] = compare_sasfit_sphere_gauss494 COMPARISON[('sasfit', 'sphere', 'gaussian')] = compare_sasfit_sphere_gauss 496 495 497 496 def compare_sasfit_sphere_schulz(): … … 512 511 title = " ".join(("sasfit", "sphere", "pd=10% schulz")) 513 512 compare(title, target, actual) 514 COMPARISON[('sasfit', 'sphere','schulz')] = compare_sasfit_sphere_schulz513 COMPARISON[('sasfit', 'sphere', 'schulz')] = compare_sasfit_sphere_schulz 515 514 516 515 def compare_sasfit_ellipsoid_schulz(): … … 533 532 title = " ".join(("sasfit", "ellipsoid", "pd=10% schulz")) 534 533 compare(title, target, actual) 535 COMPARISON[('sasfit', 'ellipsoid','schulz')] = compare_sasfit_ellipsoid_schulz534 COMPARISON[('sasfit', 'ellipsoid', 'schulz')] = compare_sasfit_ellipsoid_schulz 536 535 537 536 … … 653 652 target = Theory(Q=Q, S=SQ, Seff=SQ_EFF) 654 653 compare("sasfit ellipsoid P(Q) 60% Rp 60% Vf", target, actual); plt.show() 655 COMPARISON[('sasfit', 'ellipsoid','gaussian')] = compare_sasfit_ellipsoid_gaussian654 COMPARISON[('sasfit', 'ellipsoid', 'gaussian')] = compare_sasfit_ellipsoid_gaussian 656 655 657 656 def main():
Note: See TracChangeset
for help on using the changeset viewer.