Changeset d8ac2ad in sasmodels


Ignore:
Timestamp:
Nov 6, 2017 2:51:30 PM (7 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, 4073633
Parents:
49eb251 (diff), eb87965 (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:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/orientation/orientation.rst

    r3d40839 r82592da  
    6262care for large ranges of angle. 
    6363 
    64 Note that the form factors for asymmetric particles are also performing 
    65 numerical integrations over one or more variables, so care should be taken, 
    66 especially with very large particles or more extreme aspect ratios. Users can 
    67 experiment with the values of *Npts* and *Nsigs*, the number of steps used in the 
    68 integration and the range spanned in number of standard deviations. The 
    69 standard deviation is entered in units of degrees. For a rectangular 
    70 (uniform) distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard 
    71 deviations (this may be changed soon). 
     64.. note:: 
     65    Note that the form factors for oriented particles are also performing 
     66    numerical integrations over one or more variables, so care should be taken, 
     67    especially with very large particles or more extreme aspect ratios. In such  
     68    cases results may not be accurate, particularly at very high Q, unless the model 
     69    has been specifically coded to use limiting forms of the scattering equations. 
     70     
     71    For best numerical results keep the $\theta$ distribution narrower than the $\phi$  
     72    distribution. Thus for asymmetric particles, such as elliptical_cylinder, you may  
     73    need to reorder the sizes of the three axes to acheive the desired result.  
     74    This is due to the issues of mapping a rectangular distribution onto the  
     75    surface of a sphere. 
    7276 
    73 Where appropriate, for best numerical results, keep $a < b < c$ and the 
    74 $\theta$ distribution narrower than the $\phi$ distribution. 
     77Users can experiment with the values of *Npts* and *Nsigs*, the number of steps  
     78used in the integration and the range spanned in number of standard deviations.  
     79The standard deviation is entered in units of degrees. For a "rectangular"  
     80distribution the full width should be $\pm \sqrt(3)$ ~ 1.73 standard deviations.  
     81The new "uniform" distribution avoids this by letting you directly specify the  
     82half width. 
     83 
     84The angular distributions will be truncated outside of the range -180 to +180  
     85degrees, so beware of using saying a broad Gaussian distribution with large value 
     86of *Nsigs*, as the array of *Npts* may be truncated to many fewer points than would  
     87give a good integration,as well as becoming rather meaningless. (At some point  
     88in the future the actual polydispersity arrays may be made available to the user  
     89for inspection.) 
    7590 
    7691Some more detailed technical notes are provided in the developer section of 
    … …  
    7994*Document History* 
    8095 
    81 | 2017-10-27 Richard Heenan 
     96| 2017-11-06 Richard Heenan 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rbecded3 r82592da  
    9292 
    9393    // Compute effective radius in rotated coordinates 
    94     const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
    95     const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
    96                                    + square((r_minor+thick_rim)*qb)); 
     94    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 
     95    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 
     96                                   + square((r_minor+thick_rim)*qa)); 
    9797    const double be1 = sas_2J1x_x( qr_hat ); 
    9898    const double be2 = sas_2J1x_x( qrshell_hat ); 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r129bdc4 r82592da  
    9191          double sigma) 
    9292{ 
    93     // THIS NEEDS TESTING 
     93    // integrated 2d seems to match 1d reasonably well, except perhaps at very high Q 
    9494    // Vol1,2,3 and dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 
    9595    const double dr1 = -rhor - rhoh + rhoc + rhosolv; 
    … …  
    103103 
    104104    // Compute effective radius in rotated coordinates 
    105     const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     105    const double qr_hat = sqrt(square(r_major*qb) + square(r_minor*qa)); 
    106106    // does this need to be changed for the "missing corners" where there there is no "belt" ? 
    107     const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
    108                                    + square((r_minor+thick_rim)*qb)); 
     107    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qb) 
     108                                   + square((r_minor+thick_rim)*qa)); 
    109109    const double be1 = sas_2J1x_x( qr_hat ); 
    110110    const double be2 = sas_2J1x_x( qrshell_hat ); 

  • sasmodels/models/elliptical_cylinder.c

    rbecded3 r82592da  
    2828        //const double arg = radius_minor*sin_val; 
    2929        double inner_sum=0; 
    30         for(int j=0;j<20;j++) { 
    31             //20 gauss points for the inner integral 
    32             const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<76;j++) { 
     31            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
     32            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3333            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3434            const double be = sas_2J1x_x(q*r); 
    35             inner_sum += Gauss20Wt[j] * be * be; 
     35            inner_sum += Gauss76Wt[j] * be * be; 
    3636        } 
    3737        //now calculate the value of the inner integral 
    … …  
    6161    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    6262    // Given:    radius_major = r_ratio * radius_minor 
    63     const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     63    const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa)); 
    6464    const double be = sas_2J1x_x(qr); 
    6565    const double si = sas_sinx_x(qc*0.5*length); 

  • 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/core_shell_parallelepiped.py

    r904cd9c r393facf  
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    66The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. However at this time 
    8 the 1D calculation does **NOT** actually calculate a c face rim despite the presence of 
    9 the parameter. Some other aspects of the 1D calculation may be wrong. 
    10  
    11 .. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     7"rim" can be different on each (pair) of faces. 
     8 
    149 
    1510The form factor is normalized by the particle volume $V$ such that 
    … …  
    4136    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4237 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently. 
     38**meaning that there are "gaps" at the corners of the solid.** 
    4539 
    4640The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4741core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and shell, in the same manner as a core-shell model. 
    49  
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\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. 
     42amplitudes of the core and the slabs on the edges. 
     43 
     44the scattering amplitude is computed for a particular orientation of the core-shell 
     45parallelepiped with respect to the scattering vector and then averaged over all 
     46possible orientations, where $\alpha$ is the angle between the $z$ axis and the longest axis $C$ 
     47of the parallelepiped, $\beta$ is the angle between projection of the particle in the $xy$ detector plane and the $y$ axis. 
     48 
     49.. 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 
     57with 
     58 
     59.. math:: 
     60 
     61    S(x) = \frac{\sin \tfrac{1}{2}Q x}{\tfrac{1}{2}Q x} 
     62 
     63where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ are 
     64the scattering length of the parallelepiped core, and the rectangular slabs of 
     65thickness $t_A$, $t_B$ and $t_C$, respectively. 
     66$\rho_\text{solvent}$ is the scattering length of the solvent. 
    6567 
    6668FITTING NOTES 
    … …  
    7375known values, or you will certainly end up at a solution that is unphysical. 
    7476 
    75 Constraints must be applied during fitting to ensure that the inequality 
    76 $A < B < C$ is not violated. The calculation will not report an error, 
    77 but the results will not be correct. 
    78  
    7977The returned value is in units of |cm^-1|, on absolute scale. 
    8078 
    … …  
    9189$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    9290 
     91For 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, 
     93but the results may be not correct. 
     94 
    9395.. figure:: img/parallelepiped_angle_definition.png 
    9496 
    9597    Definition of the angles for oriented core-shell parallelepipeds. 
    9698    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
    97     rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 
     99    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the parallelepiped. 
    98100    The neutron or X-ray beam is along the $z$ axis. 
    99101 
    … …  
    109111    Equations (1), (13-14). (in German) 
    110112.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    111    lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     113   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
    112114   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    113115   =26379>`_ 

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