Changes in explore/asymint.py [a1c32c2:1820208] in sasmodels


Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/asymint.py

    ra1c32c2 r1820208  
    8686    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
    8787    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) 
    9191        return siA * siB * siC 
    9292    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
    9393    volume = a*b*c 
    9494    norm = CONTRAST**2*volume/10000 
    95     return norm, Fq 
    96  
    97 def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv): 
    98     overlapping = False 
    99     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 = CONTRAST 
    103     drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT 
    104     tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 
    105     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*siC 
    114                     + drA*(siAt-siA)*siB*siC 
    115                     + drB*siAt*(siBt-siB)*siC 
    116                     + drC*siAt*siBt*(siCt-siC)) 
    117         else: 
    118             return (dr0*siA*siB*siC 
    119                     + drA*(siAt-siA)*siB*siC 
    120                     + drB*siA*(siBt-siB)*siC 
    121                     + 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*dc 
    125     else: 
    126         volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc 
    127     norm = 1/(volume*10000) 
    12895    return norm, Fq 
    12996 
     
    217184    NORM, KERNEL = make_parallelepiped(A, B, C) 
    218185    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    219 elif shape == 'core_shell_parallelepiped': 
    220     #A, B, C = 4450, 14000, 47 
    221     #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    222     A, B, C = 6800, 114, 1380 
    223     DA, DB, DC = 2300, 21, 58 
    224     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,3 
    226     #SLD_SOLVENT,CONTRAST = 0, 4 
    227     if 1: # C shortest 
    228         B, C = C, B 
    229         DB, DC = DC, DB 
    230         SLDB, SLDC = SLDC, SLDB 
    231     elif 0: # C longest 
    232         A, C = C, A 
    233         DA, DC = DC, DA 
    234         SLDA, SLDC = SLDC, SLDA 
    235     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) 
    237186elif shape == 'paracrystal': 
    238187    LATTICE = 'bcc' 
     
    393342    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    394343    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)) 
    397344    #gridded_2d(Q, n=2**8+1) 
    398345    gridded_2d(Q, n=2**10+1) 
    399     #gridded_2d(Q, n=2**12+1) 
     346    #gridded_2d(Q, n=2**13+1) 
    400347    #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! 
    403349        print("dblquad", *scipy_dblquad_2d(Q)) 
    404350        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
Note: See TracChangeset for help on using the changeset viewer.