Changeset 2a0b2b1 in sasmodels


Ignore:
Timestamp:
Apr 14, 2017 8:30:29 AM (8 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:
fdd56a1
Parents:
9901384
Message:

restructure all 2D models to work with (qa,qb,qc) = rotate(qx,qy) rather than working with angles directly in preparation for revised jitter algorithm

Location:
sasmodels/models
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/barbell.c

    r592343f r2a0b2b1  
    1010//barbell kernel - same as dumbell 
    1111static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     12_bell_kernel(double qab, double qc, double h, double radius_bell, 
     13             double half_length) 
    1414{ 
    1515    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626    //    m = q R cos(alpha) 
    2727    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     28    const double m = radius_bell*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    3232    for (int i = 0; i < 76; i++){ 
    3333        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
    … …  
    4444 
    4545static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, 
     47    double radius_bell, double radius, double half_length) 
    4948{ 
    50     const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    51     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    52     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     49    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     50    const double bj = sas_2J1x_x(radius*qab); 
     51    const double si = sas_sinx_x(half_length*qc); 
    5352    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5453    const double Aq = bell_fq + cyl_fq; 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     85        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8786        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8887    } 
    … …  
    103102    double q, sin_alpha, cos_alpha; 
    104103    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     104    const double qab = q*sin_alpha; 
     105    const double qc = q*cos_alpha; 
    105106 
    106107    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     108    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    108109 
    109110    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r50beefe r2a0b2b1  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2_sq_bcc(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Rewriting equations for efficiency, accuracy and readability, and so 
     5    // code is reusable between 1D and 2D models. 
     6    const double a1 = +qa - qc + qb; 
     7    const double a2 = +qa + qc - qb; 
     8    const double a3 = -qa + qc + qb; 
    69 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     10    const double half_dnn = 0.5*dnn; 
     11    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     12 
     13#if 1 
     14    // Numerator: (1 - exp(a)^2)^3 
     15    //         => (-(exp(2a) - 1))^3 
     16    //         => -expm1(2a)^3 
     17    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) 
     18    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 
     19    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 
     20    const double exp_arg = exp(-arg); 
     21    const double Sq = -cube(expm1(-2.0*arg)) 
     22        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0)); 
     25#else 
     26    // Alternate form, which perhaps is more approachable 
     27    const double sinh_qd = sinh(arg); 
     28    const double cosh_qd = cosh(arg); 
     29    const double Sq = sinh_qd/(cosh_qd - cos(half_dnn*a1)) 
     30                    * sinh_qd/(cosh_qd - cos(half_dnn*a2)) 
     31                    * sinh_qd/(cosh_qd - cos(half_dnn*a3)); 
     32#endif 
     33 
     34    return Sq; 
     35} 
    1036 
    1137 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     38// occupied volume fraction calculated from lattice symmetry and sphere radius 
     39static double 
     40_bcc_volume_fraction(double radius, double dnn) 
     41{ 
     42    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2043} 
    2144 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     45static double 
     46form_volume(double radius) 
     47{ 
    4448    return sphere_volume(radius); 
    4549} 
    4650 
    4751 
    48 double Iq(double q, double dnn, 
     52static double Iq(double q, double dnn, 
    4953  double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     54  double sld, double solvent_sld) 
     55{ 
     56    // translate a point in [-1,1] to a point in [0, 2 pi] 
     57    const double phi_m = M_PI; 
     58    const double phi_b = M_PI; 
     59    // translate a point in [-1,1] to a point in [0, pi] 
     60    const double theta_m = M_PI_2; 
     61    const double theta_b = M_PI_2; 
    5162 
    52         //Volume fraction calculated from lattice symmetry and sphere radius 
    53         const double s1 = dnn/sqrt(0.75); 
    54         const double latticescale = 2.0*sphere_volume(radius/s1); 
     63    double outer_sum = 0.0; 
     64    for(int i=0; i<150; i++) { 
     65        double inner_sum = 0.0; 
     66        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     67        double sin_theta, cos_theta; 
     68        SINCOS(theta, sin_theta, cos_theta); 
     69        const double qc = q*cos_theta; 
     70        const double qab = q*sin_theta; 
     71        for(int j=0;j<150;j++) { 
     72            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     73            double sin_phi, cos_phi; 
     74            SINCOS(phi, sin_phi, cos_phi); 
     75            const double qa = qab*cos_phi; 
     76            const double qb = qab*sin_phi; 
     77            const double fq = _sq_bcc(qa, qb, qc, dnn, d_factor); 
     78            inner_sum += Gauss150Wt[j] * fq; 
     79        } 
     80        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     81        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     82    } 
     83    outer_sum *= theta_m; 
     84    const double Sq = outer_sum/(4.0*M_PI); 
     85    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5586 
    56     const double va = 0.0; 
    57     const double vb = 2.0*M_PI; 
    58     const double vaj = 0.0; 
    59     const double vbj = M_PI; 
    60  
    61     double summ = 0.0; 
    62     double answer = 0.0; 
    63         for(int i=0; i<150; i++) { 
    64                 //setup inner integral over the ellipsoidal cross-section 
    65                 double summj=0.0; 
    66                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    67                 for(int j=0;j<150;j++) { 
    68                         //20 gauss points for the inner integral 
    69                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    70                         double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    71                         summj += yyy; 
    72                 } 
    73                 //now calculate the value of the inner integral 
    74                 double answer = (vbj-vaj)/2.0*summj; 
    75  
    76                 //now calculate outer integral 
    77                 summ = summ+(Gauss150Wt[i] * answer); 
    78         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    79  
    80         answer = (vb-va)/2.0*summ; 
    81         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    82  
    83     return answer; 
     87    return _bcc_volume_fraction(radius, dnn) * Pq * Sq; 
    8488} 
    8589 
    8690 
    87 double Iqxy(double qx, double qy, 
     91static double Iqxy(double qx, double qy, 
    8892    double dnn, double d_factor, double radius, 
    8993    double sld, double solvent_sld, 
    … …  
    9296    double q, zhat, yhat, xhat; 
    9397    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     98    const double qa = q*xhat; 
     99    const double qb = q*yhat; 
     100    const double qc = q*zhat; 
    94101 
    95     const double a1 = +xhat - zhat + yhat; 
    96     const double a2 = +xhat + zhat - yhat; 
    97     const double a3 = -xhat + zhat + yhat; 
    98  
    99     const double qd = 0.5*q*dnn; 
    100     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    101     const double tanh_qd = tanh(arg); 
    102     const double cosh_qd = cosh(arg); 
    103     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    105                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    106  
    107     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    108     //the occupied volume of the lattice 
    109     const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    110     return lattice_scale * Fq; 
     102    q = sqrt(qa*qa + qb*qb + qc*qc); 
     103    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     104    const double Sq = _sq_bcc(qa, qb, qc, dnn, d_factor); 
     105    return _bcc_volume_fraction(radius, dnn) * Pq * Sq; 
    111106} 

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r2a0b2b1  
    8181.. figure:: img/parallelepiped_angle_definition.png 
    8282 
    83     Orientation of the crystal with respect to the scattering plane, when  
     83    Orientation of the crystal with respect to the scattering plane, when 
    8484    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8585 

  • sasmodels/models/capped_cylinder.c

    r592343f r2a0b2b1  
    1414//   radius_cap is the radius of the lens 
    1515//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     16//   theta is the angle of the cylinder wrt q. 
    1717static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     18_cap_kernel(double qab, double qc, double h, double radius_cap, 
     19            double half_length) 
    2020{ 
    2121    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626 
    2727    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     28    //    cos (q (R t - h + L/2) cos(theta)) 
    2929    // so turn it into: 
    3030    //    cos (m t + b) 
    3131    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     32    //    m = q R cos(theta) 
     33    //    b = q(L/2-h) cos(theta) 
     34    const double m = radius_cap*qc; // cos argument slope 
     35    const double b = (half_length-h)*qc; // cos argument intercept 
     36    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3737    double total = 0.0; 
    3838    for (int i=0; i<76 ;i++) { 
    3939        const double t = Gauss76Z[i]*zm + zb; 
    4040        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
    … …  
    5050 
    5151static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     52_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5453{ 
    55     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    56     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    57     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     54    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     55    const double bj = sas_2J1x_x(radius*qab); 
     56    const double si = sas_sinx_x(half_length*qc); 
    5857    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5958    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    101100    double total = 0.0; 
    102101    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    107         const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
    108         // sin_alpha for spherical coord integration 
    109         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     102        const double theta = Gauss76Z[i]*zm + zb; 
     103        double sin_theta, cos_theta; // slots to hold sincos function output 
     104        SINCOS(theta, sin_theta, cos_theta); 
     105        const double qab = q*sin_theta; 
     106        const double qc = q*cos_theta; 
     107        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     108        // scale by sin_theta for spherical coord integration 
     109        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110110    } 
    111111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    125125    double q, sin_alpha, cos_alpha; 
    126126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     127    const double qab = q*sin_alpha; 
     128    const double qc = q*cos_alpha; 
    127129 
    128130    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     131    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130132 
    131133    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    rb260926 r2a0b2b1  
    1 double form_volume(double radius, double thick_rim, double thick_face, double length); 
    2 double Iq(double q, 
    3           double radius, 
    4           double thick_rim, 
    5           double thick_face, 
    6           double length, 
    7           double core_sld, 
    8           double face_sld, 
    9           double rim_sld, 
    10           double solvent_sld); 
    11  
    12  
    13 double Iqxy(double qx, double qy, 
    14           double radius, 
    15           double thick_rim, 
    16           double thick_face, 
    17           double length, 
    18           double core_sld, 
    19           double face_sld, 
    20           double rim_sld, 
    21           double solvent_sld, 
    22           double theta, 
    23           double phi); 
    24  
    25  
    26 double form_volume(double radius, double thick_rim, double thick_face, double length) 
     1static double 
     2form_volume(double radius, double thick_rim, double thick_face, double length) 
    273{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
     4    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
    295} 
    306 
    317static double 
    32 bicelle_kernel(double q, 
    33               double rad, 
    34               double radthick, 
    35               double facthick, 
    36               double halflength, 
    37               double rhoc, 
    38               double rhoh, 
    39               double rhor, 
    40               double rhosolv, 
    41               double sin_alpha, 
    42               double cos_alpha) 
     8bicelle_kernel(double qab, 
     9    double qc, 
     10    double radius, 
     11    double thick_radius, 
     12    double thick_face, 
     13    double halflength, 
     14    double sld_core, 
     15    double sld_face, 
     16    double sld_rim, 
     17    double sld_solvent) 
    4318{ 
    44     const double dr1 = rhoc-rhoh; 
    45     const double dr2 = rhor-rhosolv; 
    46     const double dr3 = rhoh-rhor; 
    47     const double vol1 = M_PI*square(rad)*2.0*(halflength); 
    48     const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
    49     const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
     19    const double dr1 = sld_core-sld_face; 
     20    const double dr2 = sld_rim-sld_solvent; 
     21    const double dr3 = sld_face-sld_rim; 
     22    const double vol1 = M_PI*square(radius)*2.0*(halflength); 
     23    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 
     24    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 
    5025 
    51     const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
    52     const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
    53     const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
    54     const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
     26    const double be1 = sas_2J1x_x((radius)*qab); 
     27    const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 
     28    const double si1 = sas_sinx_x((halflength)*qc); 
     29    const double si2 = sas_sinx_x((halflength+thick_face)*qc); 
    5530 
    5631    const double t = vol1*dr1*si1*be1 + 
    … …  
    5833                     vol3*dr3*si2*be1; 
    5934 
    60     const double retval = t*t; 
    61  
    62     return retval; 
    63  
     35    return t; 
    6436} 
    6537 
    6638static double 
    67 bicelle_integration(double q, 
    68                    double rad, 
    69                    double radthick, 
    70                    double facthick, 
    71                    double length, 
    72                    double rhoc, 
    73                    double rhoh, 
    74                    double rhor, 
    75                    double rhosolv) 
     39Iq(double q, 
     40    double radius, 
     41    double thick_radius, 
     42    double thick_face, 
     43    double length, 
     44    double sld_core, 
     45    double sld_face, 
     46    double sld_rim, 
     47    double sld_solvent) 
    7648{ 
    7749    // set up the integration end points 
    … …  
    7951    const double halflength = 0.5*length; 
    8052 
    81     double summ = 0.0; 
     53    double total = 0.0; 
    8254    for(int i=0;i<N_POINTS_76;i++) { 
    83         double alpha = (Gauss76Z[i] + 1.0)*uplim; 
    84         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    85         SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
    88                              sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
     56        double sin_theta, cos_theta; // slots to hold sincos function output 
     57        SINCOS(theta, sin_theta, cos_theta); 
     58        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     59                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    9061    } 
    9162 
    9263    // calculate value of integral to return 
    93     double answer = uplim*summ; 
    94     return answer; 
     64    double answer = total*uplim; 
     65    return 1.0e-4*answer; 
    9566} 
    9667 
    9768static double 
    98 bicelle_kernel_2d(double qx, double qy, 
    99           double radius, 
    100           double thick_rim, 
    101           double thick_face, 
    102           double length, 
    103           double core_sld, 
    104           double face_sld, 
    105           double rim_sld, 
    106           double solvent_sld, 
    107           double theta, 
    108           double phi) 
     69Iqxy(double qx, double qy, 
     70    double radius, 
     71    double thick_rim, 
     72    double thick_face, 
     73    double length, 
     74    double core_sld, 
     75    double face_sld, 
     76    double rim_sld, 
     77    double solvent_sld, 
     78    double theta, 
     79    double phi) 
    10980{ 
    11081    double q, sin_alpha, cos_alpha; 
    11182    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     83    const double qab = q*sin_alpha; 
     84    const double qc = q*cos_alpha; 
    11285 
    113     double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
     86    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    11487                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     88                           solvent_sld); 
     89    return 1.0e-4*fq*fq; 
    11790} 
    118  
    119 double Iq(double q, 
    120           double radius, 
    121           double thick_rim, 
    122           double thick_face, 
    123           double length, 
    124           double core_sld, 
    125           double face_sld, 
    126           double rim_sld, 
    127           double solvent_sld) 
    128 { 
    129     double intensity = bicelle_integration(q, radius, thick_rim, thick_face, 
    130                        length, core_sld, face_sld, rim_sld, solvent_sld); 
    131     return intensity*1.0e-4; 
    132 } 
    133  
    134  
    135 double Iqxy(double qx, double qy, 
    136           double radius, 
    137           double thick_rim, 
    138           double thick_face, 
    139           double length, 
    140           double core_sld, 
    141           double face_sld, 
    142           double rim_sld, 
    143           double solvent_sld, 
    144           double theta, 
    145           double phi) 
    146 { 
    147     double intensity = bicelle_kernel_2d(qx, qy, 
    148                       radius, 
    149                       thick_rim, 
    150                       thick_face, 
    151                       length, 
    152                       core_sld, 
    153                       face_sld, 
    154                       rim_sld, 
    155                       solvent_sld, 
    156                       theta, 
    157                       phi); 
    158  
    159     return intensity; 
    160 } 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rdedcf34 r2a0b2b1  
    1717        double thick_face, 
    1818        double length, 
    19         double rhoc, 
    20         double rhoh, 
    21         double rhor, 
    22         double rhosolv) 
     19        double sld_core, 
     20        double sld_face, 
     21        double sld_rim, 
     22        double sld_solvent) 
    2323{ 
    24     double si1,si2,be1,be2; 
    2524     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2625     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    27      //    const double uplim = M_PI_4; 
    2826    const double halfheight = 0.5*length; 
    29     //const double va = 0.0; 
    30     //const double vb = 1.0; 
    31     // inner integral limits 
    32     //const double vaj=0.0; 
    33     //const double vbj=M_PI; 
    34  
    3527    const double r_major = r_minor * x_core; 
    3628    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3729    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    38     const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
    39     const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    40     const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    41     //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
    42     //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    43     //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     30    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     31    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     32    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     33    const double dr1 = vol1*(sld_core-sld_face); 
     34    const double dr2 = vol2*(sld_rim-sld_solvent); 
     35    const double dr3 = vol3*(sld_face-sld_rim); 
    4436 
    4537    //initialize integral 
    … …  
    4739    for(int i=0;i<76;i++) { 
    4840        //setup inner integral over the ellipsoidal cross-section 
    49         // since we generate these lots of times, why not store them somewhere? 
    50         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    51         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    52         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    53         double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    56         si1 = sas_sinx_x(sinarg1); 
    57         si2 = sas_sinx_x(sinarg2); 
     41        //const double va = 0.0; 
     42        //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; 
     45        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     46        const double qab = q*sin_theta; 
     47        const double qc = q*cos_theta; 
     48        const double si1 = sas_sinx_x(halfheight*qc); 
     49        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     50        double inner_sum=0.0; 
    5851        for(int j=0;j<76;j++) { 
    5952            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    60             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    61             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    62             const double rr = sqrt(r2A - r2B*cos(beta)); 
    63             double besarg1 = q*rr*sin_alpha; 
    64             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    65             be1 = sas_2J1x_x(besarg1); 
    66             be2 = sas_2J1x_x(besarg2); 
    67             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    68                                               dr2*si2*be2 + 
    69                                               dr3*si2*be1); 
     53            // inner integral limits 
     54            //const double vaj=0.0; 
     55            //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; 
     58            const double rr = sqrt(r2A - r2B*cos(phi)); 
     59            const double be1 = sas_2J1x_x(rr*qab); 
     60            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     61            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     62 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    7064        } 
    7165        //now calculate outer integral 
    … …  
    8377          double thick_face, 
    8478          double length, 
    85           double rhoc, 
    86           double rhoh, 
    87           double rhor, 
    88           double rhosolv, 
     79          double sld_core, 
     80          double sld_face, 
     81          double sld_rim, 
     82          double sld_solvent, 
    8983          double theta, 
    9084          double phi, 
    9185          double psi) 
    9286{ 
    93        // THIS NEEDS TESTING 
    9487    double q, xhat, yhat, zhat; 
    9588    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    96     const double dr1 = rhoc-rhoh; 
    97     const double dr2 = rhor-rhosolv; 
    98     const double dr3 = rhoh-rhor; 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
     92 
     93    const double dr1 = sld_core-sld_face; 
     94    const double dr2 = sld_rim-sld_solvent; 
     95    const double dr3 = sld_face-sld_rim; 
    9996    const double r_major = r_minor*x_core; 
    10097    const double halfheight = 0.5*length; 
    … …  
    104101 
    105102    // Compute effective radius in rotated coordinates 
    106     const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
    107     const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
    108                                    + square((r_minor+thick_rim)*yhat)); 
    109     const double be1 = sas_2J1x_x( q*r_hat ); 
    110     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    111     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    112     const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq; 
     103    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     104    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
     105                                   + square((r_minor+thick_rim)*qb)); 
     106    const double be1 = sas_2J1x_x( qr_hat ); 
     107    const double be2 = sas_2J1x_x( qrshell_hat ); 
     108    const double si1 = sas_sinx_x( halfheight*qc ); 
     109    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
     110    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1; 
     111    return 1.0e-4 * fq*fq; 
    115112} 
    116  

  • sasmodels/models/core_shell_cylinder.c

    r592343f r2a0b2b1  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
    5     double radius, double thickness, double length, double theta, double phi); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
    … …  
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // 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; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6460    double q, sin_alpha, cos_alpha; 
    6561    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     62    const double qab = q*sin_alpha; 
     63    const double qc = q*cos_alpha; 
    6664 
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     65    const double core_r = radius; 
     66    const double core_h = 0.5*length; 
    6967    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     68    const double shell_r = (radius + thickness); 
     69    const double shell_h = (0.5*length + thickness); 
    7270    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7371 
    74     const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha) 
    75         + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
     72    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     73        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7674    return 1.0e-4 * fq * fq; 
    7775} 

  • sasmodels/models/core_shell_ellipsoid.c

    r0a3d9b2 r2a0b2b1  
    1 double form_volume(double radius_equat_core, 
    2                    double polar_core, 
    3                    double equat_shell, 
    4                    double polar_shell); 
    5 double Iq(double q, 
    6           double radius_equat_core, 
    7           double x_core, 
    8           double thick_shell, 
    9           double x_polar_shell, 
    10           double core_sld, 
    11           double shell_sld, 
    12           double solvent_sld); 
    131 
     2// Converted from Igor function gfn4, using the same pattern as ellipsoid 
     3// for evaluating the parts of the integral. 
     4//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN 
     5//                       BY (53) & (58-59) IN CHEN AND 
     6//                       KOTLARCHYK REFERENCE 
     7// 
     8//       <OBLATE ELLIPSOID> 
     9static double 
     10_cs_ellipsoid_kernel(double qab, double qc, 
     11    double equat_core, double polar_core, 
     12    double equat_shell, double polar_shell, 
     13    double sld_core_shell, double sld_shell_solvent) 
     14{ 
     15    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc)); 
     16    const double si_core = sas_3j1x_x(qr_core); 
     17    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core; 
     18    const double fq_core = si_core*volume_core*sld_core_shell; 
    1419 
    15 double Iqxy(double qx, double qy, 
    16           double radius_equat_core, 
    17           double x_core, 
    18           double thick_shell, 
    19           double x_polar_shell, 
    20           double core_sld, 
    21           double shell_sld, 
    22           double solvent_sld, 
    23           double theta, 
    24           double phi); 
     20    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc)); 
     21    const double si_shell = sas_3j1x_x(qr_shell); 
     22    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell; 
     23    const double fq_shell = si_shell*volume_shell*sld_shell_solvent; 
    2524 
     25    return fq_core + fq_shell; 
     26} 
    2627 
    27 double form_volume(double radius_equat_core, 
    28                    double x_core, 
    29                    double thick_shell, 
    30                    double x_polar_shell) 
     28static double 
     29form_volume(double radius_equat_core, 
     30    double x_core, 
     31    double thick_shell, 
     32    double x_polar_shell) 
    3133{ 
    3234    const double equat_shell = radius_equat_core + thick_shell; 
    … …  
    3739 
    3840static double 
    39 core_shell_ellipsoid_xt_kernel(double q, 
    40           double radius_equat_core, 
    41           double x_core, 
    42           double thick_shell, 
    43           double x_polar_shell, 
    44           double core_sld, 
    45           double shell_sld, 
    46           double solvent_sld) 
     41Iq(double q, 
     42    double radius_equat_core, 
     43    double x_core, 
     44    double thick_shell, 
     45    double x_polar_shell, 
     46    double core_sld, 
     47    double shell_sld, 
     48    double solvent_sld) 
    4749{ 
    48     const double lolim = 0.0; 
    49     const double uplim = 1.0; 
    50  
    51  
    52     const double delpc = core_sld - shell_sld; //core - shell 
    53     const double delps = shell_sld - solvent_sld; //shell - solvent 
    54  
     50    const double sld_core_shell = core_sld - shell_sld; 
     51    const double sld_shell_solvent = shell_sld - solvent_sld; 
    5552 
    5653    const double polar_core = radius_equat_core*x_core; 
    … …  
    5855    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5956 
    60     double summ = 0.0;   //initialize intergral 
     57    // translate from [-1, 1] => [0, 1] 
     58    const double m = 0.5; 
     59    const double b = 0.5; 
     60    double total = 0.0;     //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    62         double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    63         double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 
    64                           polar_shell, delpc, delps, q); 
    65         summ += Gauss76Wt[i] * yyy; 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
     63        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     64        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
     65            radius_equat_core, polar_core, 
     66            equat_shell, polar_shell, 
     67            sld_core_shell, sld_shell_solvent); 
     68        total += Gauss76Wt[i] * fq * fq; 
    6669    } 
    67     summ *= 0.5*(uplim-lolim); 
     70    total *= m; 
    6871 
    6972    // convert to [cm-1] 
    70     return 1.0e-4 * summ; 
     73    return 1.0e-4 * total; 
    7174} 
    7275 
    7376static double 
    74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    75           double radius_equat_core, 
    76           double x_core, 
    77           double thick_shell, 
    78           double x_polar_shell, 
    79           double core_sld, 
    80           double shell_sld, 
    81           double solvent_sld, 
    82           double theta, 
    83           double phi) 
     77Iqxy(double qx, double qy, 
     78    double radius_equat_core, 
     79    double x_core, 
     80    double thick_shell, 
     81    double x_polar_shell, 
     82    double core_sld, 
     83    double shell_sld, 
     84    double solvent_sld, 
     85    double theta, 
     86    double phi) 
    8487{ 
    8588    double q, sin_alpha, cos_alpha; 
    8689    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     90    const double qab = q*sin_alpha; 
     91    const double qc = q*cos_alpha; 
    8792 
    88     const double sldcs = core_sld - shell_sld; 
    89     const double sldss = shell_sld- solvent_sld; 
     93    const double sld_core_shell = core_sld - shell_sld; 
     94    const double sld_shell_solvent = shell_sld - solvent_sld; 
    9095 
    9196    const double polar_core = radius_equat_core*x_core; 
    … …  
    9398    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9499 
    95     // Call the IGOR library function to get the kernel: 
    96     // MUST use gfn4 not gf2 because of the def of params. 
    97     double answer = gfn4(cos_alpha, 
    98                   radius_equat_core, 
    99                   polar_core, 
    100                   equat_shell, 
    101                   polar_shell, 
    102                   sldcs, 
    103                   sldss, 
    104                   q); 
     100    double fq = _cs_ellipsoid_kernel(qab, qc, 
     101                  radius_equat_core, polar_core, 
     102                  equat_shell, polar_shell, 
     103                  sld_core_shell, sld_shell_solvent); 
    105104 
    106105    //convert to [cm-1] 
    107     answer *= 1.0e-4; 
    108  
    109     return answer; 
     106    return 1.0e-4 * fq * fq; 
    110107} 
    111  
    112 double Iq(double q, 
    113           double radius_equat_core, 
    114           double x_core, 
    115           double thick_shell, 
    116           double x_polar_shell, 
    117           double core_sld, 
    118           double shell_sld, 
    119           double solvent_sld) 
    120 { 
    121     double intensity = core_shell_ellipsoid_xt_kernel(q, 
    122            radius_equat_core, 
    123            x_core, 
    124            thick_shell, 
    125            x_polar_shell, 
    126            core_sld, 
    127            shell_sld, 
    128            solvent_sld); 
    129  
    130     return intensity; 
    131 } 
    132  
    133  
    134 double Iqxy(double qx, double qy, 
    135           double radius_equat_core, 
    136           double x_core, 
    137           double thick_shell, 
    138           double x_polar_shell, 
    139           double core_sld, 
    140           double shell_sld, 
    141           double solvent_sld, 
    142           double theta, 
    143           double phi) 
    144 { 
    145     double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    146                        radius_equat_core, 
    147                        x_core, 
    148                        thick_shell, 
    149                        x_polar_shell, 
    150                        core_sld, 
    151                        shell_sld, 
    152                        solvent_sld, 
    153                        theta, 
    154                        phi); 
    155  
    156     return intensity; 
    157 } 

  • sasmodels/models/core_shell_ellipsoid.py

    r9802ab3 r2a0b2b1  
    2525ellipsoid. This may have some undesirable effects if the aspect ratio of the 
    2626ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$ 
    27 - which assumes spheres - will not in any case be valid.  Generating a  
    28 custom product model will enable separate effective volume fraction and effective  
     27- which assumes spheres - will not in any case be valid.  Generating a 
     28custom product model will enable separate effective volume fraction and effective 
    2929radius in the $S(q)$. 
    3030 
    … …  
    4444 
    4545.. math:: 
    46     \begin{align}     
     46    \begin{align} 
    4747    F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ 
    4848    &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) 
    49     \end{align}  
     49    \end{align} 
    5050 
    5151where 
    52   
     52 
    5353.. math:: 
    5454 
    … …  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
    79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     79For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    8080see the :ref:`elliptical-cylinder` model for further information. 
    8181 
    … …  
    139139# pylint: enable=bad-whitespace, line-too-long 
    140140 
    141 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    142           "core_shell_ellipsoid.c"] 
     141source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    143142 
    144143def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c r2a0b2b1  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    … …  
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
    … …  
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36 
    3737    const double mu = 0.5 * q * length_b; 
    38      
     38 
    3939    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     41    //double vol = length_a * length_b * length_c + 
     42    //       2.0 * thick_rim_a * length_b * length_c + 
    4343    //       2.0 * thick_rim_b * length_a * length_c + 
    4444    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     45 
    4646    // Scale sides by B 
    4747    const double a_scaled = length_a / length_b; 
    … …  
    101101            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102102            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     104 
    105105            //  correct FF : sum of square of phase factors 
    106106            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    136136    double q, zhat, yhat, xhat; 
    137137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138    const double qa = q*xhat; 
     139    const double qb = q*yhat; 
     140    const double qc = q*zhat; 
    138141 
    139142    // cspkernel in csparallelepiped recoded here 
    … …  
    160163    double tc = length_a + 2.0*thick_rim_c; 
    161164    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     165    double siA = sas_sinx_x(0.5*length_a*qa); 
     166    double siB = sas_sinx_x(0.5*length_b*qb); 
     167    double siC = sas_sinx_x(0.5*length_c*qc); 
     168    double siAt = sas_sinx_x(0.5*ta*qa); 
     169    double siBt = sas_sinx_x(0.5*tb*qb); 
     170    double siCt = sas_sinx_x(0.5*tc*qc); 
     171 
    169172 
    170173    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    174177               + drB*siA*(siBt-siB)*siC*V2 
    175178               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     179 
    177180    return 1.0e-4 * f * f; 
    178181} 

  • sasmodels/models/cylinder.c

    r592343f r2a0b2b1  
    1 double form_volume(double radius, double length); 
    2 double fq(double q, double sn, double cn,double radius, double length); 
    3 double orient_avg_1D(double q, double radius, double length); 
    4 double Iq(double q, double sld, double solvent_sld, double radius, double length); 
    5 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    6     double radius, double length, double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius<0 || v.length<0) 
    92 
    10 double form_volume(double radius, double length) 
     3static double 
     4form_volume(double radius, double length) 
    115{ 
    126    return M_PI*radius*radius*length; 
    137} 
    148 
    15 double fq(double q, double sn, double cn, double radius, double length) 
     9static double 
     10fq(double qab, double qc, double radius, double length) 
    1611{ 
    17     // precompute qr and qh to save time in the loop 
    18     const double qr = q*radius; 
    19     const double qh = q*0.5*length;  
    20     return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 
     12    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    2113} 
    2214 
    23 double orient_avg_1D(double q, double radius, double length) 
     15static double 
     16orient_avg_1D(double q, double radius, double length) 
    2417{ 
    2518    // translate a point in [-1,1] to a point in [0, pi/2] 
    2619    const double zm = M_PI_4; 
    27     const double zb = M_PI_4;  
     20    const double zb = M_PI_4; 
    2821 
    2922    double total = 0.0; 
    3023    for (int i=0; i<76 ;i++) { 
    31         const double alpha = Gauss76Z[i]*zm + zb; 
    32         double sn, cn; // slots to hold sincos function output 
    33         // alpha(theta,phi) the projection of the cylinder on the detector plane 
    34         SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn; 
     24        const double theta = Gauss76Z[i]*zm + zb; 
     25        double sin_theta, cos_theta; // slots to hold sincos function output 
     26        // theta (theta,phi) the projection of the cylinder on the detector plane 
     27        SINCOS(theta , sin_theta, cos_theta); 
     28        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3630    } 
    3731    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3933} 
    4034 
    41 double Iq(double q, 
     35static double 
     36Iq(double q, 
    4237    double sld, 
    4338    double solvent_sld, 
    … …  
    4944} 
    5045 
    51  
    52 double Iqxy(double qx, double qy, 
     46static double 
     47Iqxy(double qx, double qy, 
    5348    double sld, 
    5449    double solvent_sld, 
    … …  
    6055    double q, sin_alpha, cos_alpha; 
    6156    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    62     //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 
     57    const double qab = q*sin_alpha; 
     58    const double qc = q*cos_alpha; 
    6359    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    64     const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
     60    const double form = fq(qab, qc, radius, length); 
    6561    return 1.0e-4 * square(s * form); 
    6662} 

  • sasmodels/models/ellipsoid.c

    r3b571ae r2a0b2b1  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
    … …  
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qx, double qy, 
    4442    double sld, 
    4543    double sld_solvent, 
    … …  
    5149    double q, sin_alpha, cos_alpha; 
    5250    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     51    const double qab = q*sin_alpha; 
     52    const double qc = q*cos_alpha; 
     53 
     54    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     55    const double f = sas_3j1x_x(qr); 
    5656    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5757 
    5858    return 1.0e-4 * square(f * s); 
    5959} 
    60  

  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r2a0b2b1  
    1 double form_volume(double radius_minor, double r_ratio, double length); 
    2 double Iq(double q, double radius_minor, double r_ratio, double length, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    5             double sld, double solvent_sld, double theta, double phi, double psi); 
    6  
    7  
    8 double 
     1static double 
    92form_volume(double radius_minor, double r_ratio, double length) 
    103{ 
    … …  
    125} 
    136 
    14 double 
     7static double 
    158Iq(double q, double radius_minor, double r_ratio, double length, 
    169   double sld, double solvent_sld) 
    … …  
    6154 
    6255 
    63 double 
     56static double 
    6457Iqxy(double qx, double qy, 
    6558     double radius_minor, double r_ratio, double length, 
    … …  
    6962    double q, xhat, yhat, zhat; 
    7063    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     64    const double qa = q*xhat; 
     65    const double qb = q*yhat; 
     66    const double qc = q*zhat; 
    7167 
    7268    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7369    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    75     const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     70    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     71    const double be = sas_2J1x_x(qr); 
     72    const double si = sas_sinx_x(qc*0.5*length); 
    7773    const double Aq = be * si; 
    7874    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r2a0b2b1  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2_sq_fcc(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Rewriting equations for efficiency, accuracy and readability, and so 
     5    // code is reusable between 1D and 2D models. 
     6    const double a1 = qb + qa; 
     7    const double a2 = qa + qc; 
     8    const double a3 = qb + qc; 
    69 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     10    const double half_dnn = 0.5*dnn; 
     11    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     12 
     13    // Numerator: (1 - exp(a)^2)^3 
     14    //         => (-(exp(2a) - 1))^3 
     15    //         => -expm1(2a)^3 
     16    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) 
     17    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 
     18    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 
     19    const double exp_arg = exp(-arg); 
     20    const double Sq = -cube(expm1(-2.0*arg)) 
     21        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0)); 
     24 
     25    return Sq; 
     26} 
    927 
    1028 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     29// occupied volume fraction calculated from lattice symmetry and sphere radius 
     30static double 
     31_fcc_volume_fraction(double radius, double dnn) 
     32{ 
     33    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1934} 
    2035 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     36static double 
     37form_volume(double radius) 
     38{ 
    4339    return sphere_volume(radius); 
    4440} 
    4541 
    4642 
    47 double Iq(double q, double dnn, 
     43static double Iq(double q, double dnn, 
    4844  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     45  double sld, double solvent_sld) 
     46{ 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI; 
     49    const double phi_b = M_PI; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_2; 
     52    const double theta_b = M_PI_2; 
    5053 
    51         //Volume fraction calculated from lattice symmetry and sphere radius 
    52         const double s1 = dnn*sqrt(2.0); 
    53         const double latticescale = 4.0*sphere_volume(radius/s1); 
     54    double outer_sum = 0.0; 
     55    for(int i=0; i<150; i++) { 
     56        double inner_sum = 0.0; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        double sin_theta, cos_theta; 
     59        SINCOS(theta, sin_theta, cos_theta); 
     60        const double qc = q*cos_theta; 
     61        const double qab = q*sin_theta; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     64            double sin_phi, cos_phi; 
     65            SINCOS(phi, sin_phi, cos_phi); 
     66            const double qa = qab*cos_phi; 
     67            const double qb = qab*sin_phi; 
     68            const double fq = _sq_fcc(qa, qb, qc, dnn, d_factor); 
     69            inner_sum += Gauss150Wt[j] * fq; 
     70        } 
     71        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73    } 
     74    outer_sum *= theta_m; 
     75    const double Sq = outer_sum/(4.0*M_PI); 
     76    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5477 
    55     const double va = 0.0; 
    56     const double vb = 2.0*M_PI; 
    57     const double vaj = 0.0; 
    58     const double vbj = M_PI; 
    59  
    60     double summ = 0.0; 
    61     double answer = 0.0; 
    62         for(int i=0; i<150; i++) { 
    63                 //setup inner integral over the ellipsoidal cross-section 
    64                 double summj=0.0; 
    65                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    66                 for(int j=0;j<150;j++) { 
    67                         //20 gauss points for the inner integral 
    68                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    69                         double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    70                         summj += yyy; 
    71                 } 
    72                 //now calculate the value of the inner integral 
    73                 double answer = (vbj-vaj)/2.0*summj; 
    74  
    75                 //now calculate outer integral 
    76                 summ = summ+(Gauss150Wt[i] * answer); 
    77         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    78  
    79         answer = (vb-va)/2.0*summ; 
    80         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    81  
    82     return answer; 
     78    return _fcc_volume_fraction(radius, dnn) * Pq * Sq; 
     79} 
    8380 
    8481 
    85 } 
    86  
    87 double Iqxy(double qx, double qy, 
     82static double Iqxy(double qx, double qy, 
    8883    double dnn, double d_factor, double radius, 
    8984    double sld, double solvent_sld, 
    … …  
    9287    double q, zhat, yhat, xhat; 
    9388    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
    9492 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
    98     const double qd = 0.5*q*dnn; 
    99     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    100     const double tanh_qd = tanh(arg); 
    101     const double cosh_qd = cosh(arg); 
    102     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    103                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    105  
    106     //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 
    107  
    108     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    109     //the occupied volume of the lattice 
    110     const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    111     return lattice_scale * Fq; 
     93    q = sqrt(qa*qa + qb*qb + qc*qc); 
     94    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     95    const double Sq = _sq_fcc(qa, qb, qc, dnn, d_factor); 
     96    return _fcc_volume_fraction(radius, dnn) * Pq * Sq; 
    11297} 

  • sasmodels/models/hollow_cylinder.c

    r592343f r2a0b2b1  
    11double form_volume(double radius, double thickness, double length); 
    22double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
     3    double solvent_sld); 
    44double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
     5    double solvent_sld, double theta, double phi); 
    66 
    77//#define INVALID(v) (v.radius_core >= v.radius) 
    … …  
    1414} 
    1515 
    16  
    1716static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     17_fq(double qab, double qc, 
     18    double radius, double thickness, double length) 
    2019{ 
    21     const double qs = q*sin_val; 
    22     const double lam1 = sas_2J1x_x((radius+thickness)*qs); 
    23     const double lam2 = sas_2J1x_x(radius*qs); 
     20    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     21    const double lam2 = sas_2J1x_x(radius*qab); 
    2422    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     23    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     24    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     25    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     26    const double t2 = sas_sinx_x(0.5*length*qc); 
    2927    return psi*t2; 
    3028} 
    … …  
    4341{ 
    4442    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     43    const double upper = 1.0;        //limits of numerical integral 
    4644 
    47     double summ = 0.0;                  //initialize intergral 
     45    double summ = 0.0;            //initialize intergral 
    4846    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     47        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     48        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     49        const double form = _fq(q*sin_theta, q*cos_theta, 
     50                                radius, thickness, length); 
     51        summ += Gauss76Wt[i] * form * form; 
    5452    } 
    5553 
    … …  
    6664    double q, sin_alpha, cos_alpha; 
    6765    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    68     const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
    69         sin_alpha, cos_alpha); 
     66    const double qab = q*sin_alpha; 
     67    const double qc = q*cos_alpha; 
     68 
     69    const double form = _fq(qab, qc, radius, thickness, length); 
    7070 
    7171    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7373} 
    74  

  • sasmodels/models/parallelepiped.c

    rd605080 r2a0b2b1  
    2020{ 
    2121    const double mu = 0.5 * q * length_b; 
    22      
     22 
    2323    // Scale sides by B 
    2424    const double a_scaled = length_a / length_b; 
    2525    const double c_scaled = length_c / length_b; 
    26          
     26 
    2727    // outer integral (with gauss points), integration limits = 0, 1 
    2828    double outer_total = 0; //initialize integral 
    … …  
    6969    double q, xhat, yhat, zhat; 
    7070    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     71    const double qa = q*xhat; 
     72    const double qb = q*yhat; 
     73    const double qc = q*zhat; 
    7174 
    72     const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
    73     const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
    74     const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
     75    const double siA = sas_sinx_x(0.5*length_a*qa); 
     76    const double siB = sas_sinx_x(0.5*length_b*qb); 
     77    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7578    const double V = form_volume(length_a, length_b, length_c); 
    7679    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe r2a0b2b1  
    1 double form_volume(double radius); 
     1static double 
     2_sq_sc(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Rewriting equations for efficiency, accuracy and readability, and so 
     5    // code is reusable between 1D and 2D models. 
     6    const double a1 = qa; 
     7    const double a2 = qb; 
     8    const double a3 = qc; 
    29 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     10    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    911 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     12    // Numerator: (1 - exp(a)^2)^3 
     13    //         => (-(exp(2a) - 1))^3 
     14    //         => -expm1(2a)^3 
     15    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) 
     16    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 
     17    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 
     18    const double exp_arg = exp(-arg); 
     19    const double Sq = -cube(expm1(-2.0*arg)) 
     20        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     21          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
    1923 
    20 double form_volume(double radius) 
     24    return Sq; 
     25} 
     26 
     27// occupied volume fraction calculated from lattice symmetry and sphere radius 
     28static double 
     29_sc_volume_fraction(double radius, double dnn) 
     30{ 
     31    return sphere_volume(radius/dnn); 
     32} 
     33 
     34static double 
     35form_volume(double radius) 
    2136{ 
    2237    return sphere_volume(radius); 
    2338} 
    2439 
    25 static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     40 
     41static double Iq(double q, double dnn, 
     42  double d_factor, double radius, 
     43  double sld, double solvent_sld) 
    2744{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     45    // translate a point in [-1,1] to a point in [0, 2 pi] 
     46    const double phi_m = M_PI_4; 
     47    const double phi_b = M_PI_4; 
     48    // translate a point in [-1,1] to a point in [0, pi] 
     49    const double theta_m = M_PI_4; 
     50    const double theta_b = M_PI_4; 
    3051 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3352 
    34         double temp6 = snt; 
    35         double temp7 = -1.0*temp3*snt*cnp; 
    36         double temp8 = temp3*snt*snp; 
    37         double temp9 = temp3*cnt; 
    38         double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 
    39                                (1.0-temp4*cos((temp8))+temp5)* 
    40                                (1.0-temp4*cos((temp9))+temp5)); 
    41         return (result); 
     53    double outer_sum = 0.0; 
     54    for(int i=0; i<150; i++) { 
     55        double inner_sum = 0.0; 
     56        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     57        double sin_theta, cos_theta; 
     58        SINCOS(theta, sin_theta, cos_theta); 
     59        const double qc = q*cos_theta; 
     60        const double qab = q*sin_theta; 
     61        for(int j=0;j<150;j++) { 
     62            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     63            double sin_phi, cos_phi; 
     64            SINCOS(phi, sin_phi, cos_phi); 
     65            const double qa = qab*cos_phi; 
     66            const double qb = qab*sin_phi; 
     67            const double fq = _sq_sc(qa, qb, qc, dnn, d_factor); 
     68            inner_sum += Gauss150Wt[j] * fq; 
     69        } 
     70        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     71        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     72    } 
     73    outer_sum *= theta_m; 
     74    const double Sq = outer_sum/M_PI_2; 
     75    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     76 
     77    return _sc_volume_fraction(radius, dnn) * Pq * Sq; 
    4278} 
    4379 
    44 static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
    46 { 
    47     //Function to calculate integrand values for simple cubic structure 
    4880 
    49         double da = d_factor*dnn; 
    50         double temp1 = qq*qq*da*da; 
    51         double temp2 = cube(-expm1(-temp1)); 
    52         double temp3 = qq*dnn; 
    53         double temp4 = 2.0*exp(-0.5*temp1); 
    54         double temp5 = exp(-1.0*temp1); 
    55  
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    57  
    58         return(integrand); 
    59 } 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
     81static double Iqxy(double qx, double qy, 
     82    double dnn, double d_factor, double radius, 
     83    double sld, double solvent_sld, 
     84    double theta, double phi, double psi) 
    11285{ 
    11386    double q, zhat, yhat, xhat; 
    11487    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     88    const double qa = q*xhat; 
     89    const double qb = q*yhat; 
     90    const double qc = q*zhat; 
    11591 
    116     const double qd = q*dnn; 
    117     const double arg = 0.5*square(qd*d_factor); 
    118     const double tanh_qd = tanh(arg); 
    119     const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
    123  
    124     const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    125     //the occupied volume of the lattice 
    126     const double lattice_scale = sphere_volume(radius/dnn); 
    127     return lattice_scale * Fq; 
     92    q = sqrt(qa*qa + qb*qb + qc*qc); 
     93    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     94    const double Sq = _sq_sc(qa, qb, qc, dnn, d_factor); 
     95    return _sc_volume_fraction(radius, dnn) * Pq * Sq; 
    12896} 

  • sasmodels/models/sc_paracrystal.py

    r69e1afc r2a0b2b1  
    152152    [{}, 0.001, 10.3048], 
    153153    [{}, 0.215268, 0.00814889], 
    154     [{}, (0.414467), 0.001313289], 
     154    [{}, 0.414467, 0.001313289], 
    155155    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 
    156156    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 
    157157    ] 
    158  
    159  

  • sasmodels/models/stacked_disks.c

    r19f996b r2a0b2b1  
    11static double stacked_disks_kernel( 
    2     double q, 
     2    double qab, 
     3    double qc, 
    34    double halfheight, 
    45    double thick_layer, 
    … …  
    910    double layer_sld, 
    1011    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1312    double d) 
    1413 
    … …  
    2019    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2120 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     21    const double besarg1 = radius*qab; 
     22    //const double besarg2 = radius*qab; 
    2423 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     24    const double sinarg1 = halfheight*qc; 
     25    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2726 
    2827    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4342 
    4443    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     44    double qd_cos_alpha = d*qc; 
    4645    //d*cos_alpha is the projection of d onto q (in other words the component 
    4746    //of d that is parallel to q. 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     85        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8786                           halfheight, 
    8887                           thick_layer, 
    … …  
    9392                           layer_sld, 
    9493                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9794                           d); 
    9895        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    155152    double q, sin_alpha, cos_alpha; 
    156153    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     154    const double qab = q*sin_alpha; 
     155    const double qc = q*cos_alpha; 
    157156 
    158157    double d = 2.0 * thick_layer + thick_core; 
    159158    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     159    double answer = stacked_disks_kernel(qab, qc, 
    161160                     halfheight, 
    162161                     thick_layer, 
    … …  
    167166                     layer_sld, 
    168167                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171168                     d); 
    172169 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r2a0b2b1  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
    … …  
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qx, double qy, 
    5249    double sld, 
    5350    double sld_solvent, 
    … …  
    6158    double q, xhat, yhat, zhat; 
    6259    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    const double qa = q*xhat; 
     61    const double qb = q*yhat; 
     62    const double qc = q*zhat; 
    6363 
    64     const double r = sqrt(square(radius_equat_minor*xhat) 
    65                           + square(radius_equat_major*yhat) 
    66                           + square(radius_polar*zhat)); 
    67     const double fq = sas_3j1x_x(q*r); 
    68     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     64    const double qr = sqrt(square(radius_equat_minor*qa) 
     65                           + square(radius_equat_major*qb) 
     66                           + square(radius_polar*qc)); 
     67    const double fq = sas_3j1x_x(qr); 
     68    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     69    const double drho = (sld - sld_solvent); 
    6970 
    70     return 1.0e-4 * square(s * fq); 
     71    return 1.0e-4 * square(vol * drho * fq); 
    7172} 
    72  
Note: See TracChangeset for help on using the changeset viewer.