Changeset 1941ec6 in sasmodels for sasmodels/models


Ignore:
Timestamp:
Dec 18, 2017 1:52:27 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Children:
67cc0ff
Parents:
8224d24 (diff), 93fbc34 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Paul Kienzle <pkienzle@…> (12/18/17 13:51:43)
git-committer:
Paul Kienzle <pkienzle@…> (12/18/17 13:52:27)
Message:

Merge remote-tracking branch 'upstream/ticket-1043'

Location:
sasmodels/models
Files:
1 added
33 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/_spherepy.py

    ref07e95 r108e70e  
    8888Iq.vectorized = True  # Iq accepts an array of q values 
    8989 
    90 def Iqxy(qx, qy, sld, sld_solvent, radius): 
    91     return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius) 
    92 Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values 
    93  
    9490def sesans(z, sld, sld_solvent, radius): 
    9591    """ 

  • sasmodels/models/barbell.c

    rbecded3 r108e70e  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < 76; i++){ 
    26         const double t = Gauss76Z[i]*zm + zb; 
     25    for (int i = 0; i < GAUSS_N; i++){ 
     26        const double t = GAUSS_Z[i]*zm + zb; 
    2727        const double radical = 1.0 - t*t; 
    2828        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    2929        const double Fq = cos(m*t + b) * radical * bj; 
    30         total += Gauss76Wt[i] * Fq; 
     30        total += GAUSS_W[i] * Fq; 
    3131    } 
    3232    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    7373    const double zb = M_PI_4; 
    7474    double total = 0.0; 
    75     for (int i = 0; i < 76; i++){ 
    76         const double alpha= Gauss76Z[i]*zm + zb; 
     75    for (int i = 0; i < GAUSS_N; i++){ 
     76        const double alpha= GAUSS_Z[i]*zm + zb; 
    7777        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7878        SINCOS(alpha, sin_alpha, cos_alpha); 
    7979        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     80        total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9090 
    9191static double 
    92 Iqxy(double qab, double qc, 
     92Iqac(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    rea60e08 r108e70e  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<150; i++) { 
     83    for(int i=0; i<GAUSS_N; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     85        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    8686        double sin_theta, cos_theta; 
    8787        SINCOS(theta, sin_theta, cos_theta); 
    8888        const double qc = q*cos_theta; 
    8989        const double qab = q*sin_theta; 
    90         for(int j=0;j<150;j++) { 
    91             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<GAUSS_N;j++) { 
     91            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    9292            double sin_phi, cos_phi; 
    9393            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    9595            const double qb = qab*sin_phi; 
    9696            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
    97             inner_sum += Gauss150Wt[j] * form; 
     97            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     100        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 
    … …  
    107107 
    108108 
    109 static double Iqxy(double qa, double qb, double qc, 
     109static double Iqabc(double qa, double qb, double qc, 
    110110    double dnn, double d_factor, double radius, 
    111111    double sld, double solvent_sld) 

  • sasmodels/models/capped_cylinder.c

    rbecded3 r108e70e  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
    33         const double t = Gauss76Z[i]*zm + zb; 
     32    for (int i=0; i<GAUSS_N; i++) { 
     33        const double t = GAUSS_Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    3535        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    37         total += Gauss76Wt[i] * Fq; 
     37        total += GAUSS_W[i] * Fq; 
    3838    } 
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9595    const double zb = M_PI_4; 
    9696    double total = 0.0; 
    97     for (int i=0; i<76 ;i++) { 
    98         const double theta = Gauss76Z[i]*zm + zb; 
     97    for (int i=0; i<GAUSS_N ;i++) { 
     98        const double theta = GAUSS_Z[i]*zm + zb; 
    9999        double sin_theta, cos_theta; // slots to hold sincos function output 
    100100        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    103103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104104        // scale by sin_theta for spherical coord integration 
    105         total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
     105        total += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    115115 
    116116static double 
    117 Iqxy(double qab, double qc, 
     117Iqac(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    rbecded3 r108e70e  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<N_POINTS_76;i++) { 
    55         double theta = (Gauss76Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<GAUSS_N;i++) { 
     55        double theta = (GAUSS_Z[i] + 1.0)*uplim; 
    5656        double sin_theta, cos_theta; // slots to hold sincos function output 
    5757        SINCOS(theta, sin_theta, cos_theta); 
    5858        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5959                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += Gauss76Wt[i]*fq*fq*sin_theta; 
     60        total += GAUSS_W[i]*fq*fq*sin_theta; 
    6161    } 
    6262 
    … …  
    6767 
    6868static double 
    69 Iqxy(double qab, double qc, 
     69Iqac(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r82592da r108e70e  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<76;i++) { 
     39    for(int i=0;i<GAUSS_N;i++) { 
    4040        //setup inner integral over the ellipsoidal cross-section 
    4141        //const double va = 0.0; 
    4242        //const double vb = 1.0; 
    43         //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     43        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
    4545        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4646        const double qab = q*sin_theta; 
    … …  
    4949        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    5050        double inner_sum=0.0; 
    51         for(int j=0;j<76;j++) { 
     51        for(int j=0;j<GAUSS_N;j++) { 
    5252            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    5353            // inner integral limits 
    5454            //const double vaj=0.0; 
    5555            //const double vbj=M_PI; 
    56             //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     56            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    5858            const double rr = sqrt(r2A - r2B*cos(phi)); 
    5959            const double be1 = sas_2J1x_x(rr*qab); 
    … …  
    6161            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6262 
    63             inner_sum += Gauss76Wt[j] * fq * fq; 
     63            inner_sum += GAUSS_W[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += Gauss76Wt[i] * inner_sum; 
     66        outer_sum += GAUSS_W[i] * inner_sum; 
    6767    } 
    6868 
    … …  
    7171 
    7272static double 
    73 Iqxy(double qa, double qb, double qc, 
     73Iqabc(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r82592da r108e70e  
    77        double length) 
    88{ 
    9     return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length +  
     9    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
    1111} 
    … …  
    4747    //initialize integral 
    4848    double outer_sum = 0.0; 
    49     for(int i=0;i<76;i++) { 
     49    for(int i=0;i<GAUSS_N;i++) { 
    5050        //setup inner integral over the ellipsoidal cross-section 
    5151        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
     52        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    5454        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    5555        double inner_sum=0; 
    … …  
    5858        si1 = sas_sinx_x(sinarg1); 
    5959        si2 = sas_sinx_x(sinarg2); 
    60         for(int j=0;j<76;j++) { 
     60        for(int j=0;j<GAUSS_N;j++) { 
    6161            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    62             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
     62            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    6464            const double rr = sqrt(r2A - r2B*cos(beta)); 
    6565            double besarg1 = q*rr*sin_alpha; 
    … …  
    6767            be1 = sas_2J1x_x(besarg1); 
    6868            be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
     69            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += Gauss76Wt[i] * inner_sum; 
     74        outer_sum += GAUSS_W[i] * inner_sum; 
    7575    } 
    7676 
    … …  
    7979 
    8080static double 
    81 Iqxy(double qa, double qb, double qc, 
     81Iqabc(double qa, double qb, double qc, 
    8282          double r_minor, 
    8383          double x_core, 
    … …  
    114114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
    115115} 
    116  

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r110f69c r108e70e  
    149149    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    150150    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
     151    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"], 
    151152    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    152153    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    153154    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"], 
    154     ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"] 
    155155    ] 
    156156 

  • sasmodels/models/core_shell_cylinder.c

    rbecded3 r108e70e  
    3030    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
     32    for (int i=0; i<GAUSS_N ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 
    3737        SINCOS(theta, sin_theta,  cos_theta); 
    3838        const double qab = q*sin_theta; 
    … …  
    4040        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4141            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += Gauss76Wt[i] * fq * fq * sin_theta; 
     42        total += GAUSS_W[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4848 
    4949 
    50 double Iqxy(double qab, double qc, 
     50double Iqac(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    rbecded3 r108e70e  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<76;i++) { 
    62         const double cos_theta = Gauss76Z[i]*m + b; 
     61    for(int i=0;i<GAUSS_N;i++) { 
     62        const double cos_theta = GAUSS_Z[i]*m + b; 
    6363        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    6464        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
    … …  
    6666            equat_shell, polar_shell, 
    6767            sld_core_shell, sld_shell_solvent); 
    68         total += Gauss76Wt[i] * fq * fq; 
     68        total += GAUSS_W[i] * fq * fq; 
    6969    } 
    7070    total *= m; 
    … …  
    7575 
    7676static double 
    77 Iqxy(double qab, double qc, 
     77Iqac(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    r904cd9c r108e70e  
     1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with 
     2// c endcaps and b overlapping a.  With the proper choice of parameters, 
     3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness 
     4// and subtracting 2*thickness from length, this should match the hollow 
     5// rectangular prism.)  Set it to 0 for the documented behaviour. 
     6#define OVERLAPPING 0 
    17static double 
    28form_volume(double length_a, double length_b, double length_c, 
    39    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    410{ 
    5     //return length_a * length_b * length_c; 
    6     return length_a * length_b * length_c + 
    7            2.0 * thick_rim_a * length_b * length_c + 
    8            2.0 * thick_rim_b * length_a * length_c + 
    9            2.0 * thick_rim_c * length_a * length_b; 
     11    return 
     12#if OVERLAPPING 
     13        // Hollow rectangular prism only includes the volume of the shell 
     14        // so uncomment the next line when comparing.  Solid rectangular 
     15        // prism, or parallelepiped want filled cores, so comment when 
     16        // comparing. 
     17        //-length_a * length_b * length_c + 
     18        (length_a + 2.0*thick_rim_a) * 
     19        (length_b + 2.0*thick_rim_b) * 
     20        (length_c + 2.0*thick_rim_c); 
     21#else 
     22        length_a * length_b * length_c + 
     23        2.0 * thick_rim_a * length_b * length_c + 
     24        2.0 * length_a * thick_rim_b * length_c + 
     25        2.0 * length_a * length_b * thick_rim_c; 
     26#endif 
    1027} 
    1128 
    … …  
    2441    double thick_rim_c) 
    2542{ 
    26     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
     43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    2744    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    28     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     45    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     46    // Code rewritten (PAK) 
    2947 
    30     const double mu = 0.5 * q * length_b; 
     48    const double half_q = 0.5*q; 
    3149 
    32     //calculate volume before rescaling (in original code, but not used) 
    33     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
    34     //double vol = length_a * length_b * length_c + 
    35     //       2.0 * thick_rim_a * length_b * length_c + 
    36     //       2.0 * thick_rim_b * length_a * length_c + 
    37     //       2.0 * thick_rim_c * length_a * length_b; 
     50    const double tA = length_a + 2.0*thick_rim_a; 
     51    const double tB = length_b + 2.0*thick_rim_b; 
     52    const double tC = length_c + 2.0*thick_rim_c; 
    3853 
    39     // Scale sides by B 
    40     const double a_scaled = length_a / length_b; 
    41     const double c_scaled = length_c / length_b; 
    42  
    43     double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    44     double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
    45     double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    46  
    47     double Vin = length_a * length_b * length_c; 
    48     //double Vot = (length_a * length_b * length_c + 
    49     //            2.0 * thick_rim_a * length_b * length_c + 
    50     //            2.0 * length_a * thick_rim_b * length_c + 
    51     //            2.0 * length_a * length_b * thick_rim_c); 
    52     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    53     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    54     double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
    55  
    56     // Scale factors (note that drC is not used later) 
    57     const double drho0 = (core_sld-solvent_sld); 
    58     const double drhoA = (arim_sld-solvent_sld); 
    59     const double drhoB = (brim_sld-solvent_sld); 
    60     const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
    61  
    62  
    63     // Precompute scale factors for combining cross terms from the shape 
    64     const double scale23 = drhoA*V1; 
    65     const double scale14 = drhoB*V2; 
    66     const double scale24 = drhoC*V3; 
    67     const double scale11 = drho0*Vin; 
    68     const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
     54    // Scale factors 
     55    const double dr0 = (core_sld-solvent_sld); 
     56    const double drA = (arim_sld-solvent_sld); 
     57    const double drB = (brim_sld-solvent_sld); 
     58    const double drC = (crim_sld-solvent_sld); 
    6959 
    7060    // outer integral (with gauss points), integration limits = 0, 1 
    71     double outer_total = 0; //initialize integral 
     61    double outer_sum = 0; //initialize integral 
     62    for( int i=0; i<GAUSS_N; i++) { 
     63        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    7265 
    73     for( int i=0; i<76; i++) { 
    74         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    75         double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
     67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
     68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     69        double inner_sum = 0.0; 
     70        for(int j=0; j<GAUSS_N; j++) { 
     71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     72            double sin_beta, cos_beta; 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
     75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
     77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
    7678 
    77         // inner integral (with gauss points), integration limits = 0, 1 
    78         double inner_total = 0.0; 
    79         double inner_total_crim = 0.0; 
    80         for(int j=0; j<76; j++) { 
    81             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    82             double sin_uu, cos_uu; 
    83             SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    84             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    85             const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    86             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    87             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     79#if OVERLAPPING 
     80            const double f = dr0*siA*siB*siC 
     81                + drA*(siAt-siA)*siB*siC 
     82                + drB*siAt*(siBt-siB)*siC 
     83                + drC*siAt*siBt*(siCt-siC); 
     84#else 
     85            const double f = dr0*siA*siB*siC 
     86                + drA*(siAt-siA)*siB*siC 
     87                + drB*siA*(siBt-siB)*siC 
     88                + drC*siA*siB*(siCt-siC); 
     89#endif 
    8890 
    89             // Expression in libCylinder.c (neither drC nor Vot are used) 
    90             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
    91             const double form_crim = scale11*si1*si2; 
    92  
    93             //  correct FF : sum of square of phase factors 
    94             inner_total += Gauss76Wt[j] * form * form; 
    95             inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
     91            inner_sum += GAUSS_W[j] * f * f; 
    9692        } 
    97         inner_total *= 0.5; 
    98         inner_total_crim *= 0.5; 
     93        inner_sum *= 0.5; 
    9994        // now sum up the outer integral 
    100         const double si = sas_sinx_x(mu * c_scaled * sigma); 
    101         const double si_crim = sas_sinx_x(mu * tc * sigma); 
    102         outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
     95        outer_sum += GAUSS_W[i] * inner_sum; 
    10396    } 
    104     outer_total *= 0.5; 
     97    outer_sum *= 0.5; 
    10598 
    10699    //convert from [1e-12 A-1] to [cm-1] 
    107     return 1.0e-4 * outer_total; 
     100    return 1.0e-4 * outer_sum; 
    108101} 
    109102 
    110103static double 
    111 Iqxy(double qa, double qb, double qc, 
     104Iqabc(double qa, double qb, double qc, 
    112105    double core_sld, 
    113106    double arim_sld, 
    … …  
    128121    const double drC = crim_sld-solvent_sld; 
    129122 
    130     double Vin = length_a * length_b * length_c; 
    131     double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc) 
    132     double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    133     double V3 = 2.0 * length_a * length_b * thick_rim_c; 
    134     // As for the 1D case, Vot is not used 
    135     //double Vot = (length_a * length_b * length_c + 
    136     //              2.0 * thick_rim_a * length_b * length_c + 
    137     //              2.0 * length_a * thick_rim_b * length_c + 
    138     //              2.0 * length_a * length_b * thick_rim_c); 
    139  
    140123    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    141124    // the scaling by B. 
    142     double ta = length_a + 2.0*thick_rim_a; 
    143     double tb = length_b + 2.0*thick_rim_b; 
    144     double tc = length_c + 2.0*thick_rim_c; 
    145     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    146     double siA = sas_sinx_x(0.5*length_a*qa); 
    147     double siB = sas_sinx_x(0.5*length_b*qb); 
    148     double siC = sas_sinx_x(0.5*length_c*qc); 
    149     double siAt = sas_sinx_x(0.5*ta*qa); 
    150     double siBt = sas_sinx_x(0.5*tb*qb); 
    151     double siCt = sas_sinx_x(0.5*tc*qc); 
     125    const double tA = length_a + 2.0*thick_rim_a; 
     126    const double tB = length_b + 2.0*thick_rim_b; 
     127    const double tC = length_c + 2.0*thick_rim_c; 
     128    const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
     129    const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
     130    const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
     131    const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
     132    const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
     133    const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
    152134 
    153  
    154     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    155     // in the 1D code, but should be checked! 
    156     double f = ( dr0*siA*siB*siC*Vin 
    157                + drA*(siAt-siA)*siB*siC*V1 
    158                + drB*siA*(siBt-siB)*siC*V2 
    159                + drC*siA*siB*(siCt-siC)*V3); 
     135#if OVERLAPPING 
     136    const double f = dr0*siA*siB*siC 
     137        + drA*(siAt-siA)*siB*siC 
     138        + drB*siAt*(siBt-siB)*siC 
     139        + drC*siAt*siBt*(siCt-siC); 
     140#else 
     141    const double f = dr0*siA*siB*siC 
     142        + drA*(siAt-siA)*siB*siC 
     143        + drB*siA*(siBt-siB)*siC 
     144        + drC*siA*siB*(siCt-siC); 
     145#endif 
    160146 
    161147    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2d81cfe r10ee838  
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    66The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. However at this time the 1D 
    8 calculation does **NOT** actually calculate a c face rim despite the presence 
    9 of the parameter. Some other aspects of the 1D calculation may be wrong. 
    10  
    11 .. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    2115where $\langle \ldots \rangle$ is an average over all possible orientations 
    2216of the rectangular solid. 
    23  
    2417 
    2518The function calculated is the form factor of the rectangular solid below. 
    … …  
    4134    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4235 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently. 
     36**meaning that there are "gaps" at the corners of the solid.** 
    4537 
    4638The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4739core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and shell, in the same manner as a core-shell model. 
    49  
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta) 
    54                }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
    55     - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta) 
    56            }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
    57     \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta) 
    58                }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
    59     \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta) 
    60                }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    61  
    62 .. note:: 
    63  
    64     Why does t_B not appear in the above equation? 
    65     For the calculation of the form factor to be valid, the sides of the solid 
    66     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    67     If this inequality is not satisfied, the model will not report an error, 
    68     but the calculation will not be correct and thus the result wrong. 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
     52    &= (\rho_\text{core}-\rho_\text{solvent}) 
     53       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
     54    &+ (\rho_\text{A}-\rho_\text{solvent}) 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     56    &+ (\rho_\text{B}-\rho_\text{solvent}) 
     57        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     58    &+ (\rho_\text{C}-\rho_\text{solvent}) 
     59        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 
     60 
     61with 
     62 
     63.. math:: 
     64 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     66 
     67and 
     68 
     69.. math:: 
     70 
     71    Q_A &= \sin\alpha \sin\beta \\ 
     72    Q_B &= \sin\alpha \cos\beta \\ 
     73    Q_C &= \cos\alpha 
     74 
     75 
     76where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
     77are the scattering length of the parallelepiped core, and the rectangular 
     78slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
     79is the scattering length of the solvent. 
    6980 
    7081FITTING NOTES 
     82~~~~~~~~~~~~~ 
     83 
    7184If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    72 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    73 However, **no interparticle interference effects are included in this 
    74 calculation.** 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
    7587 
    7688There are many parameters in this model. Hold as many fixed as possible with 
    7789known values, or you will certainly end up at a solution that is unphysical. 
    7890 
    79 Constraints must be applied during fitting to ensure that the inequality 
    80 $A < B < C$ is not violated. The calculation will not report an error, 
    81 but the results will not be correct. 
    82  
    8391The returned value is in units of |cm^-1|, on absolute scale. 
    8492 
    8593NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8694based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    87 and length $(C+2t_C)$ values, after appropriately 
    88 sorting the three dimensions to give an oblate or prolate particle, to give an 
    89 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9098 
    9199For 2d data the orientation of the particle is required, described using 
    92 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    93 of the calculation and angular dispersions see :ref:`orientation` . 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
    94102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    95103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
     104 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    96109 
    97110.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    114127    Equations (1), (13-14). (in German) 
    115128.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    116    lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     129   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
    117130   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    118131   =26379>`_ 
    … …  
    175188        Return equivalent radius (ER) 
    176189    """ 
    177  
    178     # surface average radius (rough approximation) 
    179     surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 
    180  
    181     height = length_c + 2.0*thick_rim_c 
    182  
    183     ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 
    184     return 0.5 * (ddd) ** (1. / 3.) 
     190    from .parallelepiped import ER as ER_p 
     191 
     192    a = length_a + 2*thick_rim_a 
     193    b = length_b + 2*thick_rim_b 
     194    c = length_c + 2*thick_rim_c 
     195    return ER_p(a, b, c) 
    185196 
    186197# VR defaults to 1.0 
    … …  
    216227            psi_pd=10, psi_pd_n=1) 
    217228 
    218 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     229# rkh 7/4/17 add random unit test for 2d, note make all params different, 
     230# 2d values not tested against other codes or models 
    219231if 0:  # pak: model rewrite; need to update tests 
    220232    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/cylinder.c

    rbecded3 r108e70e  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<76 ;i++) { 
    24         const double theta = Gauss76Z[i]*zm + zb; 
     23    for (int i=0; i<GAUSS_N ;i++) { 
     24        const double theta = GAUSS_Z[i]*zm + zb; 
    2525        double sin_theta, cos_theta; // slots to hold sincos function output 
    2626        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2727        SINCOS(theta , sin_theta, cos_theta); 
    2828        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += Gauss76Wt[i] * form * form * sin_theta; 
     29        total += GAUSS_W[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4545 
    4646static double 
    47 Iqxy(double qab, double qc, 
     47Iqac(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    rbecded3 r108e70e  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2525    const double zm = 0.5; 
    2626    const double zb = 0.5; 
    2727    double total = 0.0; 
    28     for (int i=0;i<76;i++) { 
    29         const double u = Gauss76Z[i]*zm + zb; 
     28    for (int i=0;i<GAUSS_N;i++) { 
     29        const double u = GAUSS_Z[i]*zm + zb; 
    3030        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3131        const double f = sas_3j1x_x(q*r); 
    32         total += Gauss76Wt[i] * f * f; 
     32        total += GAUSS_W[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3939 
    4040static double 
    41 Iqxy(double qab, double qc, 
     41Iqac(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r82592da r108e70e  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<76;i++) { 
     24    for(int i=0;i<GAUSS_N;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    2727        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2828        //const double arg = radius_minor*sin_val; 
    2929        double inner_sum=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; 
     30        for(int j=0;j<GAUSS_N;j++) { 
     31            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3332            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3433            const double be = sas_2J1x_x(q*r); 
    35             inner_sum += Gauss76Wt[j] * be * be; 
     34            inner_sum += GAUSS_W[j] * be * be; 
    3635        } 
    3736        //now calculate the value of the inner integral 
    … …  
    4039        //now calculate outer integral 
    4140        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    42         outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
     41        outer_sum += GAUSS_W[i] * inner_sum * si * si; 
    4342    } 
    4443    outer_sum *= 0.5*(vb-va); 
    … …  
    5554 
    5655static double 
    57 Iqxy(double qa, double qb, double qc, 
     56Iqabc(double qa, double qb, double qc, 
    5857     double radius_minor, double r_ratio, double length, 
    5958     double sld, double solvent_sld) 

  • sasmodels/models/elliptical_cylinder.py

    r2d81cfe ra261a83  
    121121# pylint: enable=bad-whitespace, line-too-long 
    122122 
    123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 
    124           "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    125124 
    126125demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    rf728001 r108e70e  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<150; i++) { 
     55    for(int i=0; i<GAUSS_N; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     57        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    5858        double sin_theta, cos_theta; 
    5959        SINCOS(theta, sin_theta, cos_theta); 
    6060        const double qc = q*cos_theta; 
    6161        const double qab = q*sin_theta; 
    62         for(int j=0;j<150;j++) { 
    63             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<GAUSS_N;j++) { 
     63            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    6464            double sin_phi, cos_phi; 
    6565            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6767            const double qb = qab*sin_phi; 
    6868            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    69             inner_sum += Gauss150Wt[j] * form; 
     69            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 
    … …  
    8080 
    8181 
    82 static double Iqxy(double qa, double qb, double qc, 
     82static double Iqabc(double qa, double qb, double qc, 
    8383    double dnn, double d_factor, double radius, 
    8484    double sld, double solvent_sld) 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r592343f r74768cb  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<N_POINTS_76;i++) { 
    20         const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<GAUSS_N;i++) { 
     20        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
    2121        double sn, cn; 
    2222        SINCOS(zi, sn, cn); 
    2323        const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 
    2424        const double yyy = sas_2J1x_x(arg); 
    25         sum += Gauss76Wt[i] * yyy * yyy; 
     25        sum += GAUSS_W[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

    rbecded3 r108e70e  
    3838 
    3939    double summ = 0.0;            //initialize intergral 
    40     for (int i=0;i<76;i++) { 
    41         const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     40    for (int i=0;i<GAUSS_N;i++) { 
     41        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
    4242        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4343        const double form = _fq(q*sin_theta, q*cos_theta, 
    4444                                radius, thickness, length); 
    45         summ += Gauss76Wt[i] * form * form; 
     45        summ += GAUSS_W[i] * form * form; 
    4646    } 
    4747 
    … …  
    5252 
    5353static double 
    54 Iqxy(double qab, double qc, 
     54Iqac(double qab, double qc, 
    5555    double radius, double thickness, double length, 
    5656    double sld, double solvent_sld) 

  • sasmodels/models/hollow_rectangular_prism.c

    r1e7b0db0 r108e70e  
    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; 
    41     for(int i=0; i<76; i++) { 
     41    for(int i=0; i<GAUSS_N; i++) { 
    4242 
    43         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    4444        double sin_theta, cos_theta; 
    4545        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    4949 
    5050        double inner_sum = 0.0; 
    51         for(int j=0; j<76; j++) { 
     51        for(int j=0; j<GAUSS_N; j++) { 
    5252 
    53             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    5454            double sin_phi, cos_phi; 
    5555            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6666            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    6767 
    68             inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
     68            inner_sum += GAUSS_W[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
     86 
     87double Iqabc(double qa, double qb, double qc, 
     88    double sld, 
     89    double solvent_sld, 
     90    double length_a, 
     91    double b2a_ratio, 
     92    double c2a_ratio, 
     93    double thickness) 
     94{ 
     95    const double length_b = length_a * b2a_ratio; 
     96    const double length_c = length_a * c2a_ratio; 
     97    const double a_half = 0.5 * length_a; 
     98    const double b_half = 0.5 * length_b; 
     99    const double c_half = 0.5 * length_c; 
     100    const double vol_total = length_a * length_b * length_c; 
     101    const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 
     102 
     103    // Amplitude AP from eqn. (13) 
     104 
     105    const double termA1 = sas_sinx_x(qa * a_half); 
     106    const double termA2 = sas_sinx_x(qa * (a_half-thickness)); 
     107 
     108    const double termB1 = sas_sinx_x(qb * b_half); 
     109    const double termB2 = sas_sinx_x(qb * (b_half-thickness)); 
     110 
     111    const double termC1 = sas_sinx_x(qc * c_half); 
     112    const double termC2 = sas_sinx_x(qc * (c_half-thickness)); 
     113 
     114    const double AP1 = vol_total * termA1 * termB1 * termC1; 
     115    const double AP2 = vol_core * termA2 * termB2 * termC2; 
     116 
     117    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     118    const double delrho = sld - solvent_sld; 
     119 
     120    // Convert from [1e-12 A-1] to [cm-1] 
     121    return 1.0e-4 * square(delrho * (AP1-AP2)); 
     122} 

  • sasmodels/models/hollow_rectangular_prism.py

    r2d81cfe r0e55afe  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7 It computes only the 1D scattering, not the 2D. 
     7 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 **The 2D scattering intensity is not computed by this model.** 
     68For 2d data the orientation of the particle is required, described using 
     69angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     70of the calculation and angular dispersions see :ref:`orientation` . 
     71The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 
     72$\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 
     73 
     74For 2d, constraints must be applied during fitting to ensure that the inequality 
     75$A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 
     76but the results may be not correct. 
     77 
     78.. figure:: img/parallelepiped_angle_definition.png 
     79 
     80    Definition of the angles for oriented hollow rectangular prism. 
     81    Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 
     82    rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism. 
     83    The neutron or X-ray beam is along the $z$ axis. 
     84 
     85.. figure:: img/parallelepiped_angle_projection.png 
     86 
     87    Examples of the angles for oriented hollow rectangular prisms against the 
     88    detector plane. 
    6989 
    7090 
    … …  
    113133              ["thickness", "Ang", 1, [0, inf], "volume", 
    114134               "Thickness of parallelepiped"], 
     135              ["theta", "degrees", 0, [-360, 360], "orientation", 
     136               "c axis to beam angle"], 
     137              ["phi", "degrees", 0, [-360, 360], "orientation", 
     138               "rotation about beam"], 
     139              ["psi", "degrees", 0, [-360, 360], "orientation", 
     140               "rotation about c axis"], 
    115141             ] 
    116142 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rab2aea8 r74768cb  
    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 
    … …  
    2929    const double v2a = 0.0; 
    3030    const double v2b = M_PI_2;  //phi integration limits 
    31      
     31 
    3232    double outer_sum = 0.0; 
    33     for(int i=0; i<76; i++) { 
    34         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<GAUSS_N; i++) { 
     34        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    3535 
    3636        double sin_theta, cos_theta; 
    … …  
    4444 
    4545        double inner_sum = 0.0; 
    46         for(int j=0; j<76; j++) { 
    47             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<GAUSS_N; j++) { 
     47            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    4848 
    4949            double sin_phi, cos_phi; 
    … …  
    6262                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6363 
    64             inner_sum += Gauss76Wt[j] * square(AL+AT); 
     64            inner_sum += GAUSS_W[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     68        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r994d77f r74768cb  
    77 * 
    88 */ 
     9 #ifdef GAUSS_N 
     10 # undef GAUSS_N 
     11 # undef GAUSS_Z 
     12 # undef GAUSS_W 
     13 #endif 
     14 #define GAUSS_N 150 
     15 #define GAUSS_Z Gauss150Z 
     16 #define GAUSS_W Gauss150Wt 
     17 
     18// Note: using array size 152 so that it is a multiple of 4 
    919 
    1020// Gaussians 
    11 constant double Gauss150Z[150]={ 
     21constant double Gauss150Z[152]={ 
    1222        -0.9998723404457334, 
    1323        -0.9993274305065947, 
    … …  
    159169        0.9983473449340834, 
    160170        0.9993274305065947, 
    161         0.9998723404457334 
     171        0.9998723404457334, 
     172        0., 
     173        0. 
    162174}; 
    163175 
    164 constant double Gauss150Wt[150]={ 
     176constant double Gauss150Wt[152]={ 
    165177        0.0003276086705538, 
    166178        0.0007624720924706, 
    … …  
    312324        0.0011976474864367, 
    313325        0.0007624720924706, 
    314         0.0003276086705538 
     326        0.0003276086705538, 
     327        0., 
     328        0. 
    315329}; 

  • sasmodels/models/lib/gauss20.c

    r994d77f r74768cb  
    77 * 
    88 */ 
     9 #ifdef GAUSS_N 
     10 # undef GAUSS_N 
     11 # undef GAUSS_Z 
     12 # undef GAUSS_W 
     13 #endif 
     14 #define GAUSS_N 20 
     15 #define GAUSS_Z Gauss20Z 
     16 #define GAUSS_W Gauss20Wt 
    917 
    1018// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r66d119f r74768cb  
    77 * 
    88 */ 
    9 #define N_POINTS_76 76 
     9 #ifdef GAUSS_N 
     10 # undef GAUSS_N 
     11 # undef GAUSS_Z 
     12 # undef GAUSS_W 
     13 #endif 
     14 #define GAUSS_N 76 
     15 #define GAUSS_Z Gauss76Z 
     16 #define GAUSS_W Gauss76Wt 
    1017 
    1118// Gaussians 
    12 constant double Gauss76Wt[N_POINTS_76]={ 
     19constant double Gauss76Wt[76]={ 
    1320        .00126779163408536,             //0 
    1421        .00294910295364247, 
    … …  
    8996}; 
    9097 
    91 constant double Gauss76Z[N_POINTS_76]={ 
     98constant double Gauss76Z[76]={ 
    9299        -.999505948362153,              //0 
    93100        -.997397786355355, 

  • sasmodels/models/line.py

    r2d81cfe r108e70e  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
     59 
    5960def Iqxy(qx, qy, *args): 
    6061    """ 
    … …  
    6970 
    7071Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values 
     72 
     73# uncomment the following to test Iqxy in C models 
     74#del Iq, Iqxy 
     75#c_code = """ 
     76#static double Iq(double q, double b, double m) { return m*q+b; } 
     77#static double Iqxy(double qx, double qy, double b, double m) 
     78#{ return (m*qx+b)*(m*qy+b); } 
     79#""" 
    7180 
    7281def random(): 

  • sasmodels/models/parallelepiped.c

    r9b7b23f r108e70e  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<76; i++) { 
    26         const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     25    for( int i=0; i<GAUSS_N; i++) { 
     26        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    2727        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    2828 
    … …  
    3030        // corresponding to angles from 0 to pi/2. 
    3131        double inner_total = 0.0; 
    32         for(int j=0; j<76; j++) { 
    33             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     32        for(int j=0; j<GAUSS_N; j++) { 
     33            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    3434            double sin_uu, cos_uu; 
    3535            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    3636            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    3737            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
    38             inner_total += Gauss76Wt[j] * square(si1 * si2); 
     38            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
    4040        inner_total *= 0.5; 
    4141 
    4242        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    43         outer_total += Gauss76Wt[i] * inner_total * si * si; 
     43        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 
    … …  
    5353 
    5454static double 
    55 Iqxy(double qa, double qb, double qc, 
     55Iqabc(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/pringle.c

    r1e7b0db0 r74768cb  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < 76; i++) { 
    32         r = Gauss76Z[i]*zm + zb; 
     31    for (int i=0; i < GAUSS_N; i++) { 
     32        r = GAUSS_Z[i]*zm + zb; 
    3333 
    3434        const double qrs = r*q_sin_psi; 
    3535        const double qrrc = r*r*q_cos_psi; 
    3636 
    37         double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
    3838        double S, C; 
    3939        SINCOS(alpha*qrrc, S, C); 
    … …  
    8686 
    8787    double sum = 0.0; 
    88     for (int i = 0; i < 76; i++) { 
    89         double psi = Gauss76Z[i]*zm + zb; 
     88    for (int i = 0; i < GAUSS_N; i++) { 
     89        double psi = GAUSS_Z[i]*zm + zb; 
    9090        double sin_psi, cos_psi; 
    9191        SINCOS(psi, sin_psi, cos_psi); 
    … …  
    9393        double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 
    9494        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 
    95         sum += Gauss76Wt[i] * pringle_kernel; 
     95        sum += GAUSS_W[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r1e7b0db0 r108e70e  
    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; 
    30     for(int i=0; i<76; i++) { 
    31         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<GAUSS_N; i++) { 
     31        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    3232        double sin_theta, cos_theta; 
    3333        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    3636 
    3737        double inner_sum = 0.0; 
    38         for(int j=0; j<76; j++) { 
    39             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<GAUSS_N; j++) { 
     39            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    4040            double sin_phi, cos_phi; 
    4141            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    4545            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
    4646            const double AP = termA * termB * termC; 
    47             inner_sum += Gauss76Wt[j] * AP * AP; 
     47            inner_sum += GAUSS_W[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     50        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    5151    } 
    5252 
    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 Iqabc(double qa, double qb, double qc, 
     74    double sld, 
     75    double solvent_sld, 
     76    double length_a, 
     77    double b2a_ratio, 
     78    double c2a_ratio) 
     79{ 
     80    const double length_b = length_a * b2a_ratio; 
     81    const double length_c = length_a * c2a_ratio; 
     82    const double a_half = 0.5 * length_a; 
     83    const double b_half = 0.5 * length_b; 
     84    const double c_half = 0.5 * length_c; 
     85    const double volume = length_a * length_b * length_c; 
     86 
     87    // Amplitude AP from eqn. (13) 
     88 
     89    const double termA = sas_sinx_x(qa * a_half); 
     90    const double termB = sas_sinx_x(qb * b_half); 
     91    const double termC = sas_sinx_x(qc * c_half); 
     92 
     93    const double AP = termA * termB * termC; 
     94 
     95    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
     96    const double delrho = sld - solvent_sld; 
     97 
     98    // Convert from [1e-12 A-1] to [cm-1] 
     99    return 1.0e-4 * square(volume * delrho * AP); 
     100} 

  • sasmodels/models/rectangular_prism.py

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

  • sasmodels/models/sc_paracrystal.c

    rf728001 r108e70e  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<150; i++) { 
     56    for(int i=0; i<GAUSS_N; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        const double theta = GAUSS_Z[i]*theta_m + theta_b; 
    5959        double sin_theta, cos_theta; 
    6060        SINCOS(theta, sin_theta, cos_theta); 
    6161        const double qc = q*cos_theta; 
    6262        const double qab = q*sin_theta; 
    63         for(int j=0;j<150;j++) { 
    64             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<GAUSS_N;j++) { 
     64            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    6565            double sin_phi, cos_phi; 
    6666            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6868            const double qb = qab*sin_phi; 
    6969            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
    70             inner_sum += Gauss150Wt[j] * form; 
     70            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 
    … …  
    8282 
    8383static double 
    84 Iqxy(double qa, double qb, double qc, 
     84Iqabc(double qa, double qb, double qc, 
    8585    double dnn, double d_factor, double radius, 
    8686    double sld, double solvent_sld) 

  • sasmodels/models/stacked_disks.c

    rbecded3 r108e70e  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<N_POINTS_76; i++) { 
    84         double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<GAUSS_N; i++) { 
     84        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
    8585        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8686        SINCOS(zi, sin_alpha, cos_alpha); 
    … …  
    9595                           solvent_sld, 
    9696                           d); 
    97         summ += Gauss76Wt[i] * yyy * sin_alpha; 
     97        summ += GAUSS_W[i] * yyy * sin_alpha; 
    9898    } 
    9999 
    … …  
    142142 
    143143static double 
    144 Iqxy(double qab, double qc, 
     144Iqac(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    rbecded3 r108e70e  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    23     for (int i=0;i<76;i++) { 
    24         //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
    25         const double phi = Gauss76Z[i]*zm + zb; 
     23    for (int i=0;i<GAUSS_N;i++) { 
     24        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 
     25        const double phi = GAUSS_Z[i]*zm + zb; 
    2626        const double pa_sinsq_phi = pa*square(sin(phi)); 
    2727 
    … …  
    2929        const double um = 0.5; 
    3030        const double ub = 0.5; 
    31         for (int j=0;j<76;j++) { 
     31        for (int j=0;j<GAUSS_N;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(Gauss76Z[j]*um + ub); 
     33            const double usq = square(GAUSS_Z[j]*um + ub); 
    3434            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    3535            const double fq = sas_3j1x_x(q*r); 
    36             inner += Gauss76Wt[j] * fq * fq; 
     36            inner += GAUSS_W[j] * fq * fq; 
    3737        } 
    38         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     38        outer += GAUSS_W[i] * inner;  // correcting for dx later 
    3939    } 
    4040    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    … …  
    4646 
    4747static double 
    48 Iqxy(double qa, double qb, double qc, 
     48Iqabc(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 
Note: See TracChangeset for help on using the changeset viewer.