Changes in explore/asymint.py [a1c32c2:1820208] in sasmodels
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/asymint.py
ra1c32c2 r1820208 86 86 a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 87 87 def Fq(qa, qb, qc): 88 siA = env.sas_sinx_x( a*qa/2)89 siB = env.sas_sinx_x( b*qb/2)90 siC = env.sas_sinx_x( c*qc/2)88 siA = env.sas_sinx_x(0.5*a*qa/2) 89 siB = env.sas_sinx_x(0.5*b*qb/2) 90 siC = env.sas_sinx_x(0.5*c*qc/2) 91 91 return siA * siB * siC 92 92 Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 93 93 volume = a*b*c 94 94 norm = CONTRAST**2*volume/10000 95 return norm, Fq96 97 def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):98 overlapping = False99 a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)100 da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)101 slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)102 dr0 = CONTRAST103 drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT104 tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc105 def Fq(qa, qb, qc):106 siA = a*env.sas_sinx_x(a*qa/2)107 siB = b*env.sas_sinx_x(b*qb/2)108 siC = c*env.sas_sinx_x(c*qc/2)109 siAt = tA*env.sas_sinx_x(tA*qa/2)110 siBt = tB*env.sas_sinx_x(tB*qb/2)111 siCt = tC*env.sas_sinx_x(tC*qc/2)112 if overlapping:113 return (dr0*siA*siB*siC114 + drA*(siAt-siA)*siB*siC115 + drB*siAt*(siBt-siB)*siC116 + drC*siAt*siBt*(siCt-siC))117 else:118 return (dr0*siA*siB*siC119 + drA*(siAt-siA)*siB*siC120 + drB*siA*(siBt-siB)*siC121 + drC*siA*siB*(siCt-siC))122 Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)123 if overlapping:124 volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc125 else:126 volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc127 norm = 1/(volume*10000)128 95 return norm, Fq 129 96 … … 217 184 NORM, KERNEL = make_parallelepiped(A, B, C) 218 185 NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 219 elif shape == 'core_shell_parallelepiped':220 #A, B, C = 4450, 14000, 47221 #A, B, C = 445, 140, 47 # integer for the sake of mpf222 A, B, C = 6800, 114, 1380223 DA, DB, DC = 2300, 21, 58224 SLDA, SLDB, SLDC = "5", "-0.3", "11.5"225 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3226 #SLD_SOLVENT,CONTRAST = 0, 4227 if 1: # C shortest228 B, C = C, B229 DB, DC = DC, DB230 SLDB, SLDC = SLDC, SLDB231 elif 0: # C longest232 A, C = C, A233 DA, DC = DC, DA234 SLDA, SLDC = SLDC, SLDA235 NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)236 NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)237 186 elif shape == 'paracrystal': 238 187 LATTICE = 'bcc' … … 393 342 print("gauss-150", *gauss_quad_2d(Q, n=150)) 394 343 print("gauss-500", *gauss_quad_2d(Q, n=500)) 395 print("gauss-1025", *gauss_quad_2d(Q, n=1025))396 print("gauss-2049", *gauss_quad_2d(Q, n=2049))397 344 #gridded_2d(Q, n=2**8+1) 398 345 gridded_2d(Q, n=2**10+1) 399 #gridded_2d(Q, n=2**1 2+1)346 #gridded_2d(Q, n=2**13+1) 400 347 #gridded_2d(Q, n=2**15+1) 401 if shape not in ('paracrystal', 'core_shell_parallelepiped'): 402 # adaptive forms on models for which the calculations are fast enough 348 if shape != 'paracrystal': # adaptive forms are too slow! 403 349 print("dblquad", *scipy_dblquad_2d(Q)) 404 350 print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
Note: See TracChangeset
for help on using the changeset viewer.