Changeset 688d315 in sasmodels


Ignore:
Timestamp:
Nov 20, 2017 11:45:30 AM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
c11d09f, 0e55afe
Parents:
4073633 (diff), 8c7d5d5 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'ticket-776-orientation' into ticket-786

Files:
8 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_parallelepiped.py

    r4073633 r688d315  
    215215 
    216216# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
    217 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    218 tests = [[{}, 0.2, 0.533149288477], 
    219          [{}, [0.2], [0.533149288477]], 
    220          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
    221          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
    222         ] 
    223 del tests  # TODO: fix the tests 
    224 del qx, qy  # not necessary to delete, but cleaner 
     217if 0:  # pak: model rewrite; need to update tests 
     218    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     219    tests = [[{}, 0.2, 0.533149288477], 
     220            [{}, [0.2], [0.533149288477]], 
     221            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
     222            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
     223            ] 
     224    del qx, qy  # not necessary to delete, but cleaner 

  • explore/asymint.py

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

  • sasmodels/compare.py

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

  • sasmodels/models/core_shell_parallelepiped.c

    r904cd9c rfc0b7aa  
     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) 
    2845    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    … …  
    3047    const double mu = 0.5 * q * length_b; 
    3148 
    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; 
     49    // Scale sides by B 
     50    const double a_over_b = length_a / length_b; 
     51    const double c_over_b = length_c / length_b; 
    3852 
    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 
     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; 
    4656 
    4757    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 
     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; 
    5573 
    5674    // 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; 
     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); 
    6979 
    7080    // outer integral (with gauss points), integration limits = 0, 1 
    71     double outer_total = 0; //initialize integral 
    72  
     81    double outer_sum = 0; //initialize integral 
    7382    for( int i=0; i<76; i++) { 
    7483        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    7584        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    7685 
    77         // inner integral (with gauss points), integration limits = 0, 1 
    78         double inner_total = 0.0; 
    79         double inner_total_crim = 0.0; 
     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; 
    8090        for(int j=0; j<76; j++) { 
    8191            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    8292            double sin_uu, cos_uu; 
    8393            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); 
     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); 
    8898 
    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; 
     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 
    92110 
    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; 
     111            inner_sum += Gauss76Wt[j] * f * f; 
    96112        } 
    97         inner_total *= 0.5; 
    98         inner_total_crim *= 0.5; 
     113        inner_sum *= 0.5; 
    99114        // 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); 
     115        outer_sum += Gauss76Wt[i] * inner_sum; 
    103116    } 
    104     outer_total *= 0.5; 
     117    outer_sum *= 0.5; 
    105118 
    106119    //convert from [1e-12 A-1] to [cm-1] 
    107     return 1.0e-4 * outer_total; 
     120    return 1.0e-4 * outer_sum; 
    108121} 
    109122 
    … …  
    129142 
    130143    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); 
     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; 
    139159 
    140160    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    141161    // 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); 
     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); 
    152171 
    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); 
     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 
    160183 
    161184    return 1.0e-4 * f * f; 

  • 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

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

  • 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

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