Changes in / [925b3b5:ec8d4ac] in sasmodels


Ignore:
Files:
1 deleted
8 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)) 

  • sasmodels/compare.py

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

  • sasmodels/models/core_shell_parallelepiped.c

    r4493288 rc69d6d6  
    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 
    71static double 
    82form_volume(double length_a, double length_b, double length_c, 
    93    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    104{ 
    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 
     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; 
    2710} 
    2811 
    … …  
    4124    double thick_rim_c) 
    4225{ 
    43     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
     26    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    4427    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    45     // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    46     // Code rewritten (PAK) 
     28    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    4729 
    48     const double half_q = 0.5*q; 
     30    const double mu = 0.5 * q * length_b; 
    4931 
    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; 
     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; 
    5338 
    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); 
     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; 
    5969 
    6070    // outer integral (with gauss points), integration limits = 0, 1 
    61     double outer_sum = 0; //initialize integral 
     71    double outer_total = 0; //initialize integral 
     72 
    6273    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); 
     74        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     75        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    6576 
    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; 
     77        // inner integral (with gauss points), integration limits = 0, 1 
     78        double inner_total = 0.0; 
     79        double inner_total_crim = 0.0; 
    7080        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); 
     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); 
    7888 
    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 
     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; 
    9092 
    91             inner_sum += Gauss76Wt[j] * f * f; 
     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; 
    9296        } 
    93         inner_sum *= 0.5; 
     97        inner_total *= 0.5; 
     98        inner_total_crim *= 0.5; 
    9499        // now sum up the outer integral 
    95         outer_sum += Gauss76Wt[i] * inner_sum; 
     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); 
    96103    } 
    97     outer_sum *= 0.5; 
     104    outer_total *= 0.5; 
    98105 
    99106    //convert from [1e-12 A-1] to [cm-1] 
    100     return 1.0e-4 * outer_sum; 
     107    return 1.0e-4 * outer_total; 
    101108} 
    102109 
    … …  
    121128    const double drC = crim_sld-solvent_sld; 
    122129 
     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 
    123140    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    124141    // the scaling by B. 
    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); 
     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); 
    134152 
    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 
     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); 
    146160 
    147161    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r10ee838 r2d81cfe  
    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. 
     7"rim" can be different on each (pair) of faces. However at this time the 1D 
     8calculation does **NOT** actually calculate a c face rim despite the presence 
     9of 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. 
    814 
    915The form factor is normalized by the particle volume $V$ such that 
    … …  
    1521where $\langle \ldots \rangle$ is an average over all possible orientations 
    1622of the rectangular solid. 
     23 
    1724 
    1825The function calculated is the form factor of the rectangular solid below. 
    … …  
    3441    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    3542 
    36 **meaning that there are "gaps" at the corners of the solid.** 
     43**meaning that there are "gaps" at the corners of the solid.**  Again note that 
     44$t_C = 0$ currently. 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
     48amplitudes of the core and shell, in the same manner as a core-shell model. 
    4849 
    4950.. math:: 
    5051 
    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  
    61 with 
    62  
    63 .. math:: 
    64  
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
    66  
    67 and 
    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  
    76 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
    77 are the scattering length of the parallelepiped core, and the rectangular 
    78 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
    79 is the scattering length of the solvent. 
     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. 
    8069 
    8170FITTING NOTES 
    82 ~~~~~~~~~~~~~ 
    83  
    8471If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    85 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
    86 **no interparticle interference effects are included in this calculation.** 
     72value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
     73However, **no interparticle interference effects are included in this 
     74calculation.** 
    8775 
    8876There are many parameters in this model. Hold as many fixed as possible with 
    8977known values, or you will certainly end up at a solution that is unphysical. 
    9078 
     79Constraints must be applied during fitting to ensure that the inequality 
     80$A < B < C$ is not violated. The calculation will not report an error, 
     81but the results will not be correct. 
     82 
    9183The returned value is in units of |cm^-1|, on absolute scale. 
    9284 
    9385NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    9486based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    95 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96 to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     87and length $(C+2t_C)$ values, after appropriately 
     88sorting the three dimensions to give an oblate or prolate particle, to give an 
     89effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9890 
    9991For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
    101 details of the calculation and angular dispersions see :ref:`orientation`. 
     92angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     93of the calculation and angular dispersions see :ref:`orientation` . 
    10294The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    10395$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104  
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
    10996 
    11097.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    127114    Equations (1), (13-14). (in German) 
    128115.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    129    lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
     116   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
    130117   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    131118   =26379>`_ 
    … …  
    188175        Return equivalent radius (ER) 
    189176    """ 
    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) 
     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.) 
    196185 
    197186# VR defaults to 1.0 
    … …  
    227216            psi_pd=10, psi_pd_n=1) 
    228217 
    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 
     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 
    231219if 0:  # pak: model rewrite; need to update tests 
    232220    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/hollow_rectangular_prism.c

    r8de1477 r1e7b0db0  
    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  
    87 double 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 r2d81cfe  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7  
     7It computes only the 1D scattering, not the 2D. 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 For 2d data the orientation of the particle is required, described using 
    69 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    70 of the calculation and angular dispersions see :ref:`orientation` . 
    71 The 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  
    74 For 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, 
    76 but 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. 
     68**The 2D scattering intensity is not computed by this model.** 
    8969 
    9070 
    … …  
    133113              ["thickness", "Ang", 1, [0, inf], "volume", 
    134114               "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"], 
    141115             ] 
    142116 

  • sasmodels/models/rectangular_prism.c

    r8de1477 r1e7b0db0  
    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  
    73 double 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 r2d81cfe  
    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). 
     14Note also that, contrary to :ref:`parallelepiped`, it does not compute 
     15the 2D scattering. 
    1416 
    1517 
    … …  
    2426that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2527and not to the usual convention used for example in the 
    26 :ref:`parallelepiped` model. 
     28:ref:`parallelepiped` model. As the present model does not compute 
     29the 2D scattering, this has no further consequences. 
    2730 
    2831In this model the scattering from a massive parallelepiped with an 
    … …  
    6265units) *scale* represents the volume fraction (which is unitless). 
    6366 
    64 For 2d data the orientation of the particle is required, described using 
    65 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    66 of the calculation and angular dispersions see :ref:`orientation` . 
    67 The 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  
    70 For 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, 
    72 but 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  
     67**The 2D scattering intensity is not computed by this model.** 
    8668 
    8769 
    … …  
    126108              ["c2a_ratio", "", 1, [0, inf], "volume", 
    127109               "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"], 
    134110             ] 
    135111 
Note: See TracChangeset for help on using the changeset viewer.