Changeset 925b3b5 in sasmodels


Ignore:
Timestamp:
Jan 16, 2018 6:26:28 PM (7 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
c1bccff
Parents:
ec8d4ac (diff), a1c32c2 (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.
git-author:
Paul Butler <butlerpd@…> (01/16/18 18:26:28)
git-committer:
GitHub <noreply@…> (01/16/18 18:26:28)
Message:

Merge pull request #55 from SasView?/ticket-786

core_shell_parallelepiped: rewrite so that integrated 2D matches 1D. Fixes #786

Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • explore/asymint.py

    r1820208 ra1c32c2  
    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    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) 
    95128    return norm, Fq 
    96129 
    … …  
    184217    NORM, KERNEL = make_parallelepiped(A, B, C) 
    185218    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     219elif 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) 
    186237elif shape == 'paracrystal': 
    187238    LATTICE = 'bcc' 
    … …  
    342393    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    343394    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)) 
    344397    #gridded_2d(Q, n=2**8+1) 
    345398    gridded_2d(Q, n=2**10+1) 
    346     #gridded_2d(Q, n=2**13+1) 
     399    #gridded_2d(Q, n=2**12+1) 
    347400    #gridded_2d(Q, n=2**15+1) 
    348     if shape != 'paracrystal':  # adaptive forms are too slow! 
     401    if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
     402        # adaptive forms on models for which the calculations are fast enough 
    349403        print("dblquad", *scipy_dblquad_2d(Q)) 
    350404        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 

  • .travis.yml

    rce8c388 r2d09df1  
    66    env: 
    77    - PY=2.7 
     8    - DEPLOY=True 
    89  - os: linux 
    910    env: 
    … …  
    5152  on: 
    5253    branch: master 
     54    condition: $DEPLOY = True 
    5355notifications: 
    5456  slack: 

  • sasmodels/kernel_iq.c

    r6aee3ab rec8d4ac  
    173173static double 
    174174qac_apply( 
    175     QACRotation rotation, 
     175    QACRotation *rotation, 
    176176    double qx, double qy, 
    177177    double *qa_out, double *qc_out) 
    178178{ 
    179     const double dqc = rotation.R31*qx + rotation.R32*qy; 
     179    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    180180    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    181181    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    … …  
    246246static double 
    247247qabc_apply( 
    248     QABCRotation rotation, 
     248    QABCRotation *rotation, 
    249249    double qx, double qy, 
    250250    double *qa_out, double *qb_out, double *qc_out) 
    251251{ 
    252     *qa_out = rotation.R11*qx + rotation.R12*qy; 
    253     *qb_out = rotation.R21*qx + rotation.R22*qy; 
    254     *qc_out = rotation.R31*qx + rotation.R32*qy; 
     252    *qa_out = rotation->R11*qx + rotation->R12*qy; 
     253    *qb_out = rotation->R21*qx + rotation->R22*qy; 
     254    *qc_out = rotation->R31*qx + rotation->R32*qy; 
    255255} 
    256256 
    … …  
    453453  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
    454454  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 
    455   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     455  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc) 
    456456  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
    457457 
    … …  
    467467  local_values.table.psi = 0.; 
    468468  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
    469   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     469  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    470470  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    471471#endif 
Note: See TracChangeset for help on using the changeset viewer.