Changes in / [d8ac2ad:eb87965] in sasmodels


Ignore:
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • explore/asymint.py

    r49eb251 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     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) 
    12195    return norm, Fq 
    12296 
    … …  
    210184    NORM, KERNEL = make_parallelepiped(A, B, C) 
    211185    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    212 elif 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) 
    228186elif shape == 'paracrystal': 
    229187    LATTICE = 'bcc' 
    … …  
    384342    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    385343    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)) 
    388344    #gridded_2d(Q, n=2**8+1) 
    389345    gridded_2d(Q, n=2**10+1) 
    390     #gridded_2d(Q, n=2**12+1) 
     346    #gridded_2d(Q, n=2**13+1) 
    391347    #gridded_2d(Q, n=2**15+1) 
    392     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    393         # adaptive forms on models for which the calculations are fast enough 
     348    if shape != 'paracrystal':  # adaptive forms are too slow! 
    394349        print("dblquad", *scipy_dblquad_2d(Q)) 
    395350        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 

  • sasmodels/compare.py

    r376b0ee r3bfd924  
    788788        data = empty_data2D(q, resolution=res) 
    789789        data.accuracy = opts['accuracy'] 
    790         set_beam_stop(data, qmin) 
     790        set_beam_stop(data, 0.0004) 
    791791        index = ~data.mask 
    792792    else: 
    … …  
    13541354        if model_info != model_info2: 
    13551355            pars2 = randomize_pars(model_info2, pars2) 
    1356             limit_dimensions(model_info2, pars2, maxdim) 
     1356            limit_dimensions(model_info, pars2, maxdim) 
    13571357            # Share values for parameters with the same name 
    13581358            for k, v in pars.items(): 

  • sasmodels/models/core_shell_parallelepiped.c

    rfc0b7aa 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) 
    4528    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    … …  
    4730    const double mu = 0.5 * q * length_b; 
    4831 
     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; 
     38 
    4939    // Scale sides by B 
    50     const double a_over_b = length_a / length_b; 
    51     const double c_over_b = length_c / length_b; 
     40    const double a_scaled = length_a / length_b; 
     41    const double c_scaled = length_c / length_b; 
    5242 
    53     double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b; 
    54     double tB_over_b = 1+ 2.0*thick_rim_b/length_b; 
    55     double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b; 
     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 
    5646 
    5747    double Vin = length_a * length_b * length_c; 
    58 #if OVERLAPPING 
    59     const double capA_area = length_b*length_c; 
    60     const double capB_area = (length_a+2.*thick_rim_a)*length_c; 
    61     const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 
    62 #else 
    63     const double capA_area = length_b*length_c; 
    64     const double capB_area = length_a*length_c; 
    65     const double capC_area = length_a*length_b; 
    66 #endif 
    67     const double Va = length_a * capA_area; 
    68     const double Vb = length_b * capB_area; 
    69     const double Vc = length_c * capC_area; 
    70     const double Vat = Va + 2.0 * thick_rim_a * capA_area; 
    71     const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 
    72     const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 
     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 
    7355 
    7456    // Scale factors (note that drC is not used later) 
    75     const double dr0 = (core_sld-solvent_sld); 
    76     const double drA = (arim_sld-solvent_sld); 
    77     const double drB = (brim_sld-solvent_sld); 
    78     const double drC = (crim_sld-solvent_sld); 
     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; 
    7969 
    8070    // outer integral (with gauss points), integration limits = 0, 1 
    81     double outer_sum = 0; //initialize integral 
     71    double outer_total = 0; //initialize integral 
     72 
    8273    for( int i=0; i<76; i++) { 
    8374        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    8475        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    8576 
    86         // inner integral (with gauss points), integration limits = 0, pi/2 
    87         const double siC = sas_sinx_x(mu * sigma * c_over_b); 
    88         const double siCt = sas_sinx_x(mu * sigma * tC_over_b); 
    89         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; 
    9080        for(int j=0; j<76; j++) { 
    9181            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    9282            double sin_uu, cos_uu; 
    9383            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    94             const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b); 
    95             const double siB = sas_sinx_x(mu_proj * cos_uu ); 
    96             const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b); 
    97             const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); 
     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); 
    9888 
    99 #if OVERLAPPING 
    100             const double f = dr0*Vin*siA*siB*siC 
    101                 + drA*(Vat*siAt-Va*siA)*siB*siC 
    102                 + drB*siAt*(Vbt*siBt-Vb*siB)*siC 
    103                 + drC*siAt*siBt*(Vct*siCt-Vc*siC); 
    104 #else 
    105             const double f = dr0*Vin*siA*siB*siC 
    106                 + drA*(Vat*siAt-Va*siA)*siB*siC 
    107                 + drB*siA*(Vbt*siBt-Vb*siB)*siC 
    108                 + drC*siA*siB*(Vct*siCt-Vc*siC); 
    109 #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; 
    11092 
    111             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; 
    11296        } 
    113         inner_sum *= 0.5; 
     97        inner_total *= 0.5; 
     98        inner_total_crim *= 0.5; 
    11499        // now sum up the outer integral 
    115         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); 
    116103    } 
    117     outer_sum *= 0.5; 
     104    outer_total *= 0.5; 
    118105 
    119106    //convert from [1e-12 A-1] to [cm-1] 
    120     return 1.0e-4 * outer_sum; 
     107    return 1.0e-4 * outer_total; 
    121108} 
    122109 
    … …  
    142129 
    143130    double Vin = length_a * length_b * length_c; 
    144 #if OVERLAPPING 
    145     const double capA_area = length_b*length_c; 
    146     const double capB_area = (length_a+2.*thick_rim_a)*length_c; 
    147     const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 
    148 #else 
    149     const double capA_area = length_b*length_c; 
    150     const double capB_area = length_a*length_c; 
    151     const double capC_area = length_a*length_b; 
    152 #endif 
    153     const double Va = length_a * capA_area; 
    154     const double Vb = length_b * capB_area; 
    155     const double Vc = length_c * capC_area; 
    156     const double Vat = Va + 2.0 * thick_rim_a * capA_area; 
    157     const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 
    158     const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 
     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); 
    159139 
    160140    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    161141    // the scaling by B. 
    162     const double tA = length_a + 2.0*thick_rim_a; 
    163     const double tB = length_b + 2.0*thick_rim_b; 
    164     const double tC = length_c + 2.0*thick_rim_c; 
    165     const double siA = sas_sinx_x(0.5*length_a*qa); 
    166     const double siB = sas_sinx_x(0.5*length_b*qb); 
    167     const double siC = sas_sinx_x(0.5*length_c*qc); 
    168     const double siAt = sas_sinx_x(0.5*tA*qa); 
    169     const double siBt = sas_sinx_x(0.5*tB*qb); 
    170     const double siCt = 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); 
    171152 
    172 #if OVERLAPPING 
    173     const double f = dr0*Vin*siA*siB*siC 
    174         + drA*(Vat*siAt-Va*siA)*siB*siC 
    175         + drB*siAt*(Vbt*siBt-Vb*siB)*siC 
    176         + drC*siAt*siBt*(Vct*siCt-Vc*siC); 
    177 #else 
    178     const double f = dr0*Vin*siA*siB*siC 
    179         + drA*(Vat*siAt-Va*siA)*siB*siC 
    180         + drB*siA*(Vbt*siBt-Vb*siB)*siC 
    181         + drC*siA*siB*(Vct*siCt-Vc*siC); 
    182 #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); 
    183160 
    184161    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r393facf reda8b30  
    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. 
    8  
     7"rim" can be different on each (pair) of faces. However at this time 
     8the 1D calculation does **NOT** actually calculate a c face rim despite the presence of 
     9the 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. 
    914 
    1015The form factor is normalized by the particle volume $V$ such that 
    … …  
    3641    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    3742 
    38 **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. 
    3945 
    4046The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4147core-shell intensity being calculated as the square of the sum of the 
    42 amplitudes of the core and the slabs on the edges. 
    43  
    44 the scattering amplitude is computed for a particular orientation of the core-shell 
    45 parallelepiped with respect to the scattering vector and then averaged over all 
    46 possible orientations, where $\alpha$ is the angle between the $z$ axis and the longest axis $C$ 
    47 of the parallelepiped, $\beta$ is the angle between projection of the particle in the $xy$ detector plane and the $y$ axis. 
     48amplitudes of the core and shell, in the same manner as a core-shell model. 
    4849 
    4950.. math:: 
    50     \begin{align*} 
    51     F(Q)&=A B C (\rho_\text{core}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta)S(B \sin\alpha \cos\beta)S(C \cos\alpha) \\ 
    52     &+ 2t_A B C (\rho_\text{A}-\rho_\text{solvent})  \left[S((A+t_A) \sin\alpha \sin\beta)-S(A \sin\alpha \sin\beta)\right] S(B \sin\alpha \cos\beta) S(C \cos\alpha)\\ 
    53     &+ 2 A t_B C (\rho_\text{B}-\rho_\text{solvent})  S(A \sin\alpha \sin\beta) \left[S((B+t_B) \sin\alpha \cos\beta)-S(B \sin\alpha \cos\beta)\right] S(C \cos\alpha)\\ 
    54     &+ 2 A B t_C (\rho_\text{C}-\rho_\text{solvent}) S(A \sin\alpha \sin\beta) S(B \sin\alpha \cos\beta) \left[S((C+t_C) \cos\alpha)-S(C \cos\alpha)\right] 
    55     \end{align*} 
    56  
    57 with 
    58  
    59 .. math:: 
    60  
    61     S(x) = \frac{\sin \tfrac{1}{2}Q x}{\tfrac{1}{2}Q x} 
    62  
    63 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ are 
    64 the scattering length of the parallelepiped core, and the rectangular slabs of 
    65 thickness $t_A$, $t_B$ and $t_C$, respectively. 
    66 $\rho_\text{solvent}$ is the scattering length of the solvent. 
     51 
     52    F_{a}(Q,\alpha,\beta)= 
     53    \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
     54    - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
     55    \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
     56    \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
     57 
     58.. note:: 
     59 
     60    Why does t_B not appear in the above equation? 
     61    For the calculation of the form factor to be valid, the sides of the solid 
     62    MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
     63    If this inequality is not satisfied, the model will not report an error, 
     64    but the calculation will not be correct and thus the result wrong. 
    6765 
    6866FITTING NOTES 
    … …  
    7573known values, or you will certainly end up at a solution that is unphysical. 
    7674 
     75Constraints must be applied during fitting to ensure that the inequality 
     76$A < B < C$ is not violated. The calculation will not report an error, 
     77but the results will not be correct. 
     78 
    7779The returned value is in units of |cm^-1|, on absolute scale. 
    7880 
    … …  
    8991$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    9092 
    91 For 2d, constraints must be applied during fitting to ensure that the inequality 
    92 $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
    93 but the results may be not correct. 
    94  
    9593.. figure:: img/parallelepiped_angle_definition.png 
    9694 
    9795    Definition of the angles for oriented core-shell parallelepipeds. 
    9896    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
    99     rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the parallelepiped. 
     97    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 
    10098    The neutron or X-ray beam is along the $z$ axis. 
    10199 
    … …  
    111109    Equations (1), (13-14). (in German) 
    112110.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    113    lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
     111   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
    114112   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    115113   =26379>`_ 

  • 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

    r393facf r30b60d2  
    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 
    … …  
    201175         [{}, [0.2], [0.76687283098]], 
    202176        ] 
     177 

  • 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

    r393facf r31df0c9  
    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.