Changeset 0f113fb in sasmodels


Ignore:
Timestamp:
Dec 4, 2017 6:18:21 AM (6 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:
c11d09f, a189283
Parents:
10ee838 (diff), 791281c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-786

Files:
24 edited

Legend:

Unmodified
Added
Removed

  • extra/pylint.rc

    rb669b49 r6e68289  
    2121# List of plugins (as comma separated values of python modules names) to load, 
    2222# usually to register additional checkers. 
    23 load-plugins=pylint_numpy,pylint_pyopencl,pylint_sas 
     23#load-plugins=pylint_numpy,pylint_scipy,pylint_pyopencl,pylint_sas 
    2424 
    2525# Use multiple processes to speed up Pylint. 
    … …  
    280280# (useful for modules/projects where namespaces are manipulated during runtime 
    281281# and thus existing member attributes cannot be deduced by static analysis 
    282 ignored-modules=numpy,np,numpy.random, 
    283     bumps,sas, 
     282ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special 
    284283 
    285284# List of classes names for which member attributes should not be checked 

  • sasmodels/models/_spherepy.py

    r2d81cfe ref07e95  
    4040John Wiley and Sons, New York, (1955) 
    4141 
    42 *2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.* 
     42* **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 
    4343""" 
    4444 

  • sasmodels/models/be_polyelectrolyte.py

    r2d81cfe ref07e95  
    6767* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    6868* **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 
    69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** 
    70   October 07, 2016 
     69* **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** October 07, 2016 
    7170""" 
    7271 

  • sasmodels/models/fractal_core_shell.py

    r2d81cfe ref07e95  
    5555 
    5656* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    57 * **Last Modified by:** Paul Butler and Paul Kienzle **on:** November 27, 2016 
    58 * **Last Reviewed by:** Paul Butler and Paul Kienzle **on:** November 27, 2016 
     57* **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 
     58* **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 
    5959""" 
    6060 

  • sasmodels/models/multilayer_vesicle.py

    r2d81cfe ref07e95  
    107107* **Last Modified by:** Paul Kienzle **Date:** Feb 7, 2017 
    108108* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 
    109  
    110109""" 
    111110 

  • sasmodels/models/parallelepiped.py

    r2d81cfe ref07e95  
    167167---------------------------- 
    168168 
    169 * **Author:** This model is based on form factor calculations implemented 
    170   in a c-library provided by the NIST Center for Neutron Research (Kline, 2006). 
     169* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    171170* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017 
    172171* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017 

  • sasmodels/models/polymer_micelle.py

    r2d81cfe r791281c  
    2828.. math:: 
    2929    P(q) &= N^2\beta^2_s\Phi(qr)^2 + N\beta^2_cP_c(q) 
    30             + 2N^2\beta_s\beta_cS_{sc}s_c(q) + N(N-1)\beta_c^2S_{cc}(q) \\ 
     30            + 2N^2\beta_s\beta_cS_{sc}(q) + N(N-1)\beta_c^2S_{cc}(q) \\ 
    3131    \beta_s &= V_\text{core}(\rho_\text{core} - \rho_\text{solvent}) \\ 
    3232    \beta_c &= V_\text{corona}(\rho_\text{corona} - \rho_\text{solvent}) 
    … …  
    6969 
    7070J Pedersen, *J. Appl. Cryst.*, 33 (2000) 637-640 
     71 
     72* **Modified by:** Richard Heenan **Date:** March 20, 2016 
     73* **Verified by:** Paul Kienzle **Date:** November 29, 2017 
     74* **Description modified by:** Paul Kienzle **Date:** November 29, 2017 
     75* **Description reviewed by:** Steve King **Date:** November 30, 2017 
    7176""" 
    7277 

  • sasmodels/models/pringle.py

    r2d81cfe ref07e95  
    4343Derivation by Stefan Alexandru Rautu. 
    4444 
    45 **Author:** Andrew Jackson **on:** 2008 
    46  
    47 **Last Modified by:** Wojciech Wpotrzebowski **on:** March 20, 2016 
    48  
    49 **Last Reviewed by:** Andrew Jackson **on:** September 26, 2016 
     45* **Author:** Andrew Jackson **Date:** 2008 
     46* **Last Modified by:** Wojciech Wpotrzebowski **Date:** March 20, 2016 
     47* **Last Reviewed by:** Andrew Jackson **Date:** September 26, 2016 
    5048""" 
    5149 

  • sasmodels/models/raspberry.py

    r2d81cfe ref07e95  
    102102Science*, 343(1) (2010) 36-41 
    103103 
    104 **Author:** Andrew Jackson **on:** 2008 
    105  
    106 **Modified by:** Andrew Jackson **on:** March 20, 2016 
    107  
    108 **Reviewed by:** Andrew Jackson **on:** March 20, 2016 
     104* **Author:** Andrew Jackson **Date:** 2008 
     105* **Modified by:** Andrew Jackson **Date:** March 20, 2016 
     106* **Reviewed by:** Andrew Jackson **Date:** March 20, 2016 
    109107""" 
    110108 

  • sasmodels/models/sphere.py

    r2d81cfe ref07e95  
    4040John Wiley and Sons, New York, (1955) 
    4141 
    42 *2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.* 
     42* **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 
    4343""" 
    4444 

  • sasmodels/models/spinodal.py

    r2d81cfe ref07e95  
    2626 
    2727* **Author:** Dirk Honecker **Date:** Oct 7, 2016 
    28 * **Last Modified by:** 
    29 * **Last Reviewed by:** 
    3028""" 
    3129 

  • sasmodels/models/stacked_disks.py

    r2d81cfe ref07e95  
    106106 
    107107* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    108 * **Last Modified by:** Paul Butler and Paul Kienzle **on:** November 26, 2016 
    109 * **Last Reviewed by:** Paul Butler and Paul Kienzle **on:** November 26, 2016 
     108* **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016 
     109* **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016 
    110110""" 
    111111 

  • sasmodels/models/triaxial_ellipsoid.py

    r2d81cfe r37f08d2  
    162162    Returns the effective radius used in the S*P calculation 
    163163    """ 
    164     import numpy as np 
    165164    from .ellipsoid import ER as ellipsoid_ER 
    166165 

  • sasmodels/models/two_lorentzian.py

    r2d81cfe ref07e95  
    2727None. 
    2828 
    29 **Author:** NIST IGOR/DANSE **on:** pre 2010 
    30  
    31 **Last Modified by:** Piotr rozyczko **on:** January 29, 2016 
    32  
    33 **Last Reviewed by:** Paul Butler **on:** March 21, 2016 
     29* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     30* **Last Modified by:** Piotr rozyczko **Date:** January 29, 2016 
     31* **Last Reviewed by:** Paul Butler **Date:** March 21, 2016 
    3432""" 
    3533 

  • sasmodels/models/two_power_law.py

    r2d81cfe ref07e95  
    3737None. 
    3838 
    39 **Author:** NIST IGOR/DANSE **on:** pre 2010 
    40  
    41 **Last Modified by:** Wojciech Wpotrzebowski **on:** February 18, 2016 
    42  
    43 **Last Reviewed by:** Paul Butler **on:** March 21, 2016 
     39* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     40* **Last Modified by:** Wojciech Wpotrzebowski **Date:** February 18, 2016 
     41* **Last Reviewed by:** Paul Butler **Date:** March 21, 2016 
    4442""" 
    4543 

  • sasmodels/models/vesicle.py

    r2d81cfe ref07e95  
    5656Sons, New York, (1955) 
    5757 
    58 **Author:** NIST IGOR/DANSE **on:** pre 2010 
    59  
    60 **Last Modified by:** Paul Butler **on:** March 20, 2016 
    61  
    62 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 
     58* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
     59* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
     60* **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
    6361""" 
    6462 

  • 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)) 

  • sasmodels/compare.py

    r2d81cfe r0e55afe  
    693693        data = empty_data2D(q, resolution=res) 
    694694        data.accuracy = opts['accuracy'] 
    695         set_beam_stop(data, 0.0004) 
     695        set_beam_stop(data, qmin) 
    696696        index = ~data.mask 
    697697    else: 
    … …  
    12741274        if model_info != model_info2: 
    12751275            pars2 = randomize_pars(model_info2, pars2) 
    1276             limit_dimensions(model_info, pars2, maxdim) 
     1276            limit_dimensions(model_info2, pars2, maxdim) 
    12771277            # Share values for parameters with the same name 
    12781278            for k, v in pars.items(): 

  • sasmodels/models/core_shell_parallelepiped.c

    r904cd9c r4493288  
     1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with 
     2// c endcaps and b overlapping a.  With the proper choice of parameters, 
     3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness 
     4// and subtracting 2*thickness from length, this should match the hollow 
     5// rectangular prism.)  Set it to 0 for the documented behaviour. 
     6#define OVERLAPPING 0 
    17static double 
    28form_volume(double length_a, double length_b, double length_c, 
    39    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    410{ 
    5     //return length_a * length_b * length_c; 
    6     return length_a * length_b * length_c + 
    7            2.0 * thick_rim_a * length_b * length_c + 
    8            2.0 * thick_rim_b * length_a * length_c + 
    9            2.0 * thick_rim_c * length_a * length_b; 
     11    return 
     12#if OVERLAPPING 
     13        // Hollow rectangular prism only includes the volume of the shell 
     14        // so uncomment the next line when comparing.  Solid rectangular 
     15        // prism, or parallelepiped want filled cores, so comment when 
     16        // comparing. 
     17        //-length_a * length_b * length_c + 
     18        (length_a + 2.0*thick_rim_a) * 
     19        (length_b + 2.0*thick_rim_b) * 
     20        (length_c + 2.0*thick_rim_c); 
     21#else 
     22        length_a * length_b * length_c + 
     23        2.0 * thick_rim_a * length_b * length_c + 
     24        2.0 * length_a * thick_rim_b * length_c + 
     25        2.0 * length_a * length_b * thick_rim_c; 
     26#endif 
    1027} 
    1128 
    … …  
    2441    double thick_rim_c) 
    2542{ 
    26     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
     43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    2744    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    28     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     45    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     46    // Code rewritten (PAK) 
    2947 
    30     const double mu = 0.5 * q * length_b; 
     48    const double half_q = 0.5*q; 
    3149 
    32     //calculate volume before rescaling (in original code, but not used) 
    33     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
    34     //double vol = length_a * length_b * length_c + 
    35     //       2.0 * thick_rim_a * length_b * length_c + 
    36     //       2.0 * thick_rim_b * length_a * length_c + 
    37     //       2.0 * thick_rim_c * length_a * length_b; 
     50    const double tA = length_a + 2.0*thick_rim_a; 
     51    const double tB = length_b + 2.0*thick_rim_b; 
     52    const double tC = length_c + 2.0*thick_rim_c; 
    3853 
    39     // Scale sides by B 
    40     const double a_scaled = length_a / length_b; 
    41     const double c_scaled = length_c / length_b; 
    42  
    43     double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    44     double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
    45     double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    46  
    47     double Vin = length_a * length_b * length_c; 
    48     //double Vot = (length_a * length_b * length_c + 
    49     //            2.0 * thick_rim_a * length_b * length_c + 
    50     //            2.0 * length_a * thick_rim_b * length_c + 
    51     //            2.0 * length_a * length_b * thick_rim_c); 
    52     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    53     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    54     double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
    55  
    56     // Scale factors (note that drC is not used later) 
    57     const double drho0 = (core_sld-solvent_sld); 
    58     const double drhoA = (arim_sld-solvent_sld); 
    59     const double drhoB = (brim_sld-solvent_sld); 
    60     const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
    61  
    62  
    63     // Precompute scale factors for combining cross terms from the shape 
    64     const double scale23 = drhoA*V1; 
    65     const double scale14 = drhoB*V2; 
    66     const double scale24 = drhoC*V3; 
    67     const double scale11 = drho0*Vin; 
    68     const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
     54    // Scale factors 
     55    const double dr0 = (core_sld-solvent_sld); 
     56    const double drA = (arim_sld-solvent_sld); 
     57    const double drB = (brim_sld-solvent_sld); 
     58    const double drC = (crim_sld-solvent_sld); 
    6959 
    7060    // outer integral (with gauss points), integration limits = 0, 1 
    71     double outer_total = 0; //initialize integral 
     61    double outer_sum = 0; //initialize integral 
     62    for( int i=0; i<76; i++) { 
     63        const double cos_alpha = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    7265 
    73     for( int i=0; i<76; i++) { 
    74         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    75         double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
     67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
     68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     69        double inner_sum = 0.0; 
     70        for(int j=0; j<76; j++) { 
     71            const double beta = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     72            double sin_beta, cos_beta; 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
     75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
     77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
    7678 
    77         // inner integral (with gauss points), integration limits = 0, 1 
    78         double inner_total = 0.0; 
    79         double inner_total_crim = 0.0; 
    80         for(int j=0; j<76; j++) { 
    81             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    82             double sin_uu, cos_uu; 
    83             SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    84             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    85             const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    86             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    87             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     79#if OVERLAPPING 
     80            const double f = dr0*siA*siB*siC 
     81                + drA*(siAt-siA)*siB*siC 
     82                + drB*siAt*(siBt-siB)*siC 
     83                + drC*siAt*siBt*(siCt-siC); 
     84#else 
     85            const double f = dr0*siA*siB*siC 
     86                + drA*(siAt-siA)*siB*siC 
     87                + drB*siA*(siBt-siB)*siC 
     88                + drC*siA*siB*(siCt-siC); 
     89#endif 
    8890 
    89             // Expression in libCylinder.c (neither drC nor Vot are used) 
    90             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
    91             const double form_crim = scale11*si1*si2; 
    92  
    93             //  correct FF : sum of square of phase factors 
    94             inner_total += Gauss76Wt[j] * form * form; 
    95             inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
     91            inner_sum += Gauss76Wt[j] * f * f; 
    9692        } 
    97         inner_total *= 0.5; 
    98         inner_total_crim *= 0.5; 
     93        inner_sum *= 0.5; 
    9994        // now sum up the outer integral 
    100         const double si = sas_sinx_x(mu * c_scaled * sigma); 
    101         const double si_crim = sas_sinx_x(mu * tc * sigma); 
    102         outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
     95        outer_sum += Gauss76Wt[i] * inner_sum; 
    10396    } 
    104     outer_total *= 0.5; 
     97    outer_sum *= 0.5; 
    10598 
    10699    //convert from [1e-12 A-1] to [cm-1] 
    107     return 1.0e-4 * outer_total; 
     100    return 1.0e-4 * outer_sum; 
    108101} 
    109102 
    … …  
    128121    const double drC = crim_sld-solvent_sld; 
    129122 
    130     double Vin = length_a * length_b * length_c; 
    131     double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc) 
    132     double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    133     double V3 = 2.0 * length_a * length_b * thick_rim_c; 
    134     // As for the 1D case, Vot is not used 
    135     //double Vot = (length_a * length_b * length_c + 
    136     //              2.0 * thick_rim_a * length_b * length_c + 
    137     //              2.0 * length_a * thick_rim_b * length_c + 
    138     //              2.0 * length_a * length_b * thick_rim_c); 
    139  
    140123    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    141124    // the scaling by B. 
    142     double ta = length_a + 2.0*thick_rim_a; 
    143     double tb = length_b + 2.0*thick_rim_b; 
    144     double tc = length_c + 2.0*thick_rim_c; 
    145     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    146     double siA = sas_sinx_x(0.5*length_a*qa); 
    147     double siB = sas_sinx_x(0.5*length_b*qb); 
    148     double siC = sas_sinx_x(0.5*length_c*qc); 
    149     double siAt = sas_sinx_x(0.5*ta*qa); 
    150     double siBt = sas_sinx_x(0.5*tb*qb); 
    151     double siCt = sas_sinx_x(0.5*tc*qc); 
     125    const double tA = length_a + 2.0*thick_rim_a; 
     126    const double tB = length_b + 2.0*thick_rim_b; 
     127    const double tC = length_c + 2.0*thick_rim_c; 
     128    const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
     129    const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
     130    const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
     131    const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
     132    const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
     133    const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
    152134 
    153  
    154     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    155     // in the 1D code, but should be checked! 
    156     double f = ( dr0*siA*siB*siC*Vin 
    157                + drA*(siAt-siA)*siB*siC*V1 
    158                + drB*siA*(siBt-siB)*siC*V2 
    159                + drC*siA*siB*(siCt-siC)*V3); 
     135#if OVERLAPPING 
     136    const double f = dr0*siA*siB*siC 
     137        + drA*(siAt-siA)*siB*siC 
     138        + drB*siAt*(siBt-siB)*siC 
     139        + drC*siAt*siBt*(siCt-siC); 
     140#else 
     141    const double f = dr0*siA*siB*siC 
     142        + drA*(siAt-siA)*siB*siC 
     143        + drB*siA*(siBt-siB)*siC 
     144        + drC*siA*siB*(siCt-siC); 
     145#endif 
    160146 
    161147    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2d81cfe r10ee838  
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    66The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. However at this time the 1D 
    8 calculation does **NOT** actually calculate a c face rim despite the presence 
    9 of the parameter. Some other aspects of the 1D calculation may be wrong. 
    10  
    11 .. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    2115where $\langle \ldots \rangle$ is an average over all possible orientations 
    2216of the rectangular solid. 
    23  
    2417 
    2518The function calculated is the form factor of the rectangular solid below. 
    … …  
    4134    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4235 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently. 
     36**meaning that there are "gaps" at the corners of the solid.** 
    4537 
    4638The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4739core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and shell, in the same manner as a core-shell model. 
    49  
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta) 
    54                }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
    55     - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta) 
    56            }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
    57     \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta) 
    58                }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
    59     \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta) 
    60                }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    61  
    62 .. note:: 
    63  
    64     Why does t_B not appear in the above equation? 
    65     For the calculation of the form factor to be valid, the sides of the solid 
    66     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    67     If this inequality is not satisfied, the model will not report an error, 
    68     but the calculation will not be correct and thus the result wrong. 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
     52    &= (\rho_\text{core}-\rho_\text{solvent}) 
     53       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
     54    &+ (\rho_\text{A}-\rho_\text{solvent}) 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     56    &+ (\rho_\text{B}-\rho_\text{solvent}) 
     57        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     58    &+ (\rho_\text{C}-\rho_\text{solvent}) 
     59        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 
     60 
     61with 
     62 
     63.. math:: 
     64 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     66 
     67and 
     68 
     69.. math:: 
     70 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
     74 
     75 
     76where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
     77are the scattering length of the parallelepiped core, and the rectangular 
     78slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
     79is the scattering length of the solvent. 
    6980 
    7081FITTING NOTES 
     82~~~~~~~~~~~~~ 
     83 
    7184If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    72 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    73 However, **no interparticle interference effects are included in this 
    74 calculation.** 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
    7587 
    7688There are many parameters in this model. Hold as many fixed as possible with 
    7789known values, or you will certainly end up at a solution that is unphysical. 
    7890 
    79 Constraints must be applied during fitting to ensure that the inequality 
    80 $A < B < C$ is not violated. The calculation will not report an error, 
    81 but the results will not be correct. 
    82  
    8391The returned value is in units of |cm^-1|, on absolute scale. 
    8492 
    8593NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8694based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    87 and length $(C+2t_C)$ values, after appropriately 
    88 sorting the three dimensions to give an oblate or prolate particle, to give an 
    89 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9098 
    9199For 2d data the orientation of the particle is required, described using 
    92 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    93 of the calculation and angular dispersions see :ref:`orientation` . 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
    94102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    95103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
     104 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    96109 
    97110.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    114127    Equations (1), (13-14). (in German) 
    115128.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    116    lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     129   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
    117130   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    118131   =26379>`_ 
    … …  
    175188        Return equivalent radius (ER) 
    176189    """ 
    177  
    178     # surface average radius (rough approximation) 
    179     surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 
    180  
    181     height = length_c + 2.0*thick_rim_c 
    182  
    183     ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 
    184     return 0.5 * (ddd) ** (1. / 3.) 
     190    from .parallelepiped import ER as ER_p 
     191 
     192    a = length_a + 2*thick_rim_a 
     193    b = length_b + 2*thick_rim_b 
     194    c = length_c + 2*thick_rim_c 
     195    return ER_p(a, b, c) 
    185196 
    186197# VR defaults to 1.0 
    … …  
    216227            psi_pd=10, psi_pd_n=1) 
    217228 
    218 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     229# rkh 7/4/17 add random unit test for 2d, note make all params different, 
     230# 2d values not tested against other codes or models 
    219231if 0:  # pak: model rewrite; need to update tests 
    220232    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/hollow_rectangular_prism.c

    r1e7b0db0 r8de1477  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a,  
     2double Iq(double q, double sld, double solvent_sld, double length_a, 
    33          double b2a_ratio, double c2a_ratio, double thickness); 
    44 
    … …  
    3737    const double v2a = 0.0; 
    3838    const double v2b = M_PI_2;  //phi integration limits 
    39      
     39 
    4040    double outer_sum = 0.0; 
    4141    for(int i=0; i<76; i++) { 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
     86 
     87double Iqxy(double qa, double qb, double qc, 
     88    double sld, 
     89    double solvent_sld, 
     90    double length_a, 
     91    double b2a_ratio, 
     92    double c2a_ratio, 
     93    double thickness) 
     94{ 
     95    const double length_b = length_a * b2a_ratio; 
     96    const double length_c = length_a * c2a_ratio; 
     97    const double a_half = 0.5 * length_a; 
     98    const double b_half = 0.5 * length_b; 
     99    const double c_half = 0.5 * length_c; 
     100    const double vol_total = length_a * length_b * length_c; 
     101    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 
     102 
     103    // Amplitude AP from eqn. (13) 
     104 
     105    const double termA1 = sas_sinx_x(qa * a_half); 
     106    const double termA2 = sas_sinx_x(qa * (a_half-thickness)); 
     107 
     108    const double termB1 = sas_sinx_x(qb * b_half); 
     109    const double termB2 = sas_sinx_x(qb * (b_half-thickness)); 
     110 
     111    const double termC1 = sas_sinx_x(qc * c_half); 
     112    const double termC2 = sas_sinx_x(qc * (c_half-thickness)); 
     113 
     114    const double AP1 = vol_total * termA1 * termB1 * termC1; 
     115    const double AP2 = vol_core * termA2 * termB2 * termC2; 
     116 
     117    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     118    const double delrho = sld - solvent_sld; 
     119 
     120    // Convert from [1e-12 A-1] to [cm-1] 
     121    return 1.0e-4 * square(delrho * (AP1-AP2)); 
     122} 

  • sasmodels/models/hollow_rectangular_prism.py

    r2d81cfe r0e55afe  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7 It computes only the 1D scattering, not the 2D. 
     7 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 **The 2D scattering intensity is not computed by this model.** 
     68For 2d data the orientation of the particle is required, described using 
     69angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     70of the calculation and angular dispersions see :ref:`orientation` . 
     71The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 
     72$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 
     73 
     74For 2d, constraints must be applied during fitting to ensure that the inequality 
     75$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     76but the results may be not correct. 
     77 
     78.. figure:: img/parallelepiped_angle_definition.png 
     79 
     80    Definition of the angles for oriented hollow rectangular prism. 
     81    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     82    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism. 
     83    The neutron or X-ray beam is along the $z$ axis. 
     84 
     85.. figure:: img/parallelepiped_angle_projection.png 
     86 
     87    Examples of the angles for oriented hollow rectangular prisms against the 
     88    detector plane. 
    6989 
    7090 
    … …  
    113133              ["thickness", "Ang", 1, [0, inf], "volume", 
    114134               "Thickness of parallelepiped"], 
     135              ["theta", "degrees", 0, [-360, 360], "orientation", 
     136               "c axis to beam angle"], 
     137              ["phi", "degrees", 0, [-360, 360], "orientation", 
     138               "rotation about beam"], 
     139              ["psi", "degrees", 0, [-360, 360], "orientation", 
     140               "rotation about c axis"], 
    115141             ] 
    116142 

  • sasmodels/models/rectangular_prism.c

    r1e7b0db0 r8de1477  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a,  
     2double Iq(double q, double sld, double solvent_sld, double length_a, 
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2626    const double v2a = 0.0; 
    2727    const double v2b = M_PI_2;  //phi integration limits 
    28      
     28 
    2929    double outer_sum = 0.0; 
    3030    for(int i=0; i<76; i++) { 
    … …  
    5353    double answer = 0.5*(v1b-v1a)*outer_sum; 
    5454 
    55     // Normalize by Pi (Eqn. 16).  
    56     // The term (ABC)^2 does not appear because it was introduced before on  
     55    // Normalize by Pi (Eqn. 16). 
     56    // The term (ABC)^2 does not appear because it was introduced before on 
    5757    // the definitions of termA, termB, termC. 
    58     // The factor 2 appears because the theta integral has been defined between  
     58    // The factor 2 appears because the theta integral has been defined between 
    5959    // 0 and pi/2, instead of 0 to pi. 
    6060    answer /= M_PI_2; //Form factor P(q) 
    … …  
    6464    answer *= square((sld-solvent_sld)*volume); 
    6565 
    66     // Convert from [1e-12 A-1] to [cm-1]  
     66    // Convert from [1e-12 A-1] to [cm-1] 
    6767    answer *= 1.0e-4; 
    6868 
    6969    return answer; 
    7070} 
     71 
     72 
     73double Iqxy(double qa, double qb, double qc, 
     74    double sld, 
     75    double solvent_sld, 
     76    double length_a, 
     77    double b2a_ratio, 
     78    double c2a_ratio) 
     79{ 
     80    const double length_b = length_a * b2a_ratio; 
     81    const double length_c = length_a * c2a_ratio; 
     82    const double a_half = 0.5 * length_a; 
     83    const double b_half = 0.5 * length_b; 
     84    const double c_half = 0.5 * length_c; 
     85    const double volume = length_a * length_b * length_c; 
     86 
     87    // Amplitude AP from eqn. (13) 
     88 
     89    const double termA = sas_sinx_x(qa * a_half); 
     90    const double termB = sas_sinx_x(qb * b_half); 
     91    const double termC = sas_sinx_x(qc * c_half); 
     92 
     93    const double AP = termA * termB * termC; 
     94 
     95    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     96    const double delrho = sld - solvent_sld; 
     97 
     98    // Convert from [1e-12 A-1] to [cm-1] 
     99    return 1.0e-4 * square(volume * delrho * AP); 
     100} 

  • sasmodels/models/rectangular_prism.py

    r2d81cfe r0e55afe  
    1212the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 
    1313to *a* will generate a distribution of cubes of different sizes). 
    14 Note also that, contrary to :ref:`parallelepiped`, it does not compute 
    15 the 2D scattering. 
    1614 
    1715 
    … …  
    2624that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2725and not to the usual convention used for example in the 
    28 :ref:`parallelepiped` model. As the present model does not compute 
    29 the 2D scattering, this has no further consequences. 
     26:ref:`parallelepiped` model. 
    3027 
    3128In this model the scattering from a massive parallelepiped with an 
    … …  
    6562units) *scale* represents the volume fraction (which is unitless). 
    6663 
    67 **The 2D scattering intensity is not computed by this model.** 
     64For 2d data the orientation of the particle is required, described using 
     65angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     66of the calculation and angular dispersions see :ref:`orientation` . 
     67The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 
     68$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 
     69 
     70For 2d, constraints must be applied during fitting to ensure that the inequality 
     71$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     72but the results may be not correct. 
     73 
     74.. figure:: img/parallelepiped_angle_definition.png 
     75 
     76    Definition of the angles for oriented core-shell parallelepipeds. 
     77    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     78    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 
     79    The neutron or X-ray beam is along the $z$ axis. 
     80 
     81.. figure:: img/parallelepiped_angle_projection.png 
     82 
     83    Examples of the angles for oriented rectangular prisms against the 
     84    detector plane. 
     85 
    6886 
    6987 
    … …  
    108126              ["c2a_ratio", "", 1, [0, inf], "volume", 
    109127               "Ratio sides c/a"], 
     128              ["theta", "degrees", 0, [-360, 360], "orientation", 
     129               "c axis to beam angle"], 
     130              ["phi", "degrees", 0, [-360, 360], "orientation", 
     131               "rotation about beam"], 
     132              ["psi", "degrees", 0, [-360, 360], "orientation", 
     133               "rotation about c axis"], 
    110134             ] 
    111135 
Note: See TracChangeset for help on using the changeset viewer.