Changeset 49eb251 in sasmodels for explore/asymint.py


Ignore:
Timestamp:
Nov 6, 2017 10:49:46 AM (5 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
74768cb, d8ac2ad
Parents:
393facf
Message:

add core_shell_parallelepiped to asymint (and fix parallelepiped)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/asymint.py

    r1820208 r49eb251  
    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(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) 
     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) 
    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 
     97def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv): 
     98    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
     99    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc) 
     100    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc) 
     101    drV0 = CONTRAST*a*b*c 
     102    dra, drb, drc = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT 
     103    Aa, Ab, Ac = b*c, a*c, a*b 
     104    Ta, Tb, Tc = a + 2*da, b + 2*db, c + 2*dc 
     105    drVa, drVb, drVc = dra*a*Aa, drb*b*Ab, drc*c*Ac 
     106    drVTa, drVTb, drVTc = dra*Ta*Aa, drb*Tb*Ab, drc*Tc*Ac 
     107    def Fq(qa, qb, qc): 
     108        siA = env.sas_sinx_x(a*qa/2) 
     109        siB = env.sas_sinx_x(b*qb/2) 
     110        siC = env.sas_sinx_x(c*qc/2) 
     111        siAt = env.sas_sinx_x(Ta*qa/2) 
     112        siBt = env.sas_sinx_x(Tb*qb/2) 
     113        siCt = env.sas_sinx_x(Tc*qc/2) 
     114        return (drV0*siA*siB*siC 
     115            + (drVTa*siAt-drVa*siA)*siB*siC 
     116            + siA*(drVTb*siBt-drVb*siB)*siC 
     117            + siA*siB*(drVTc*siCt-drVc*siC)) 
     118    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
     119    volume = a*b*c + 2*da*Aa + 2*db*Ab + 2*dc*Ac 
     120    norm = 1/(volume*10000) 
    95121    return norm, Fq 
    96122 
     
    184210    NORM, KERNEL = make_parallelepiped(A, B, C) 
    185211    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     212elif shape == 'core_shell_parallelepiped': 
     213    #A, B, C = 4450, 14000, 47 
     214    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
     215    A, B, C = 6800, 114, 1380 
     216    DA, DB, DC = 2300, 21, 58 
     217    SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
     218    if 1: # C shortest 
     219        B, C = C, B 
     220        DB, DC = DC, DB 
     221        SLDB, SLDC = SLDC, SLDB 
     222    elif 0: # C longest 
     223        A, C = C, A 
     224        DA, DC = DC, DA 
     225        SLDA, SLDC = SLDC, SLDA 
     226    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     227    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    186228elif shape == 'paracrystal': 
    187229    LATTICE = 'bcc' 
     
    342384    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    343385    print("gauss-500", *gauss_quad_2d(Q, n=500)) 
     386    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
     387    print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
    344388    #gridded_2d(Q, n=2**8+1) 
    345389    gridded_2d(Q, n=2**10+1) 
    346     #gridded_2d(Q, n=2**13+1) 
     390    #gridded_2d(Q, n=2**12+1) 
    347391    #gridded_2d(Q, n=2**15+1) 
    348     if shape != 'paracrystal':  # adaptive forms are too slow! 
     392    if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
     393        # adaptive forms on models for which the calculations are fast enough 
    349394        print("dblquad", *scipy_dblquad_2d(Q)) 
    350395        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
Note: See TracChangeset for help on using the changeset viewer.