source: sasmodels/sasmodels/models/cylinder.c @ 50e1e40

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 50e1e40 was 50e1e40, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

use lib functions to simplify barbell, capped_cylinder, cylinder, ellipsoid, triaxial_ellipsoid; fix barbell calculation

  • Property mode set to 100644
File size: 1.9 KB
Line 
1double form_volume(double radius, double length);
2double Iq(double q, double sld, double solvent_sld, double radius, double length);
3double Iqxy(double qx, double qy, double sld, double solvent_sld,
4    double radius, double length, double theta, double phi);
5
6double form_volume(double radius, double length)
7{
8    return M_PI*radius*radius*length;
9}
10
11double Iq(double q,
12    double sld,
13    double solvent_sld,
14    double radius,
15    double length)
16{
17    // TODO: return NaN if radius<0 or length<0?
18    // precompute qr and qh to save time in the loop
19    const double qr = q*radius;
20    const double qh = q*0.5*length;
21
22    // translate a point in [-1,1] to a point in [0, pi/2]
23    const double zm = M_PI_4;
24    const double zb = M_PI_4;
25
26    double total = 0.0;
27    for (int i=0; i<76 ;i++) {
28        const double alpha = Gauss76Z[i]*zm + zb;
29        double sn, cn;
30        SINCOS(alpha, sn, cn);
31        const double fq = sinc(qh*cn) * J1c(qr*sn);
32        total += Gauss76Wt[i] * fq*fq * sn;
33    }
34    // translate dx in [-1,1] to dx in [lower,upper]
35    const double form = total*zm;
36    const double s = (sld - solvent_sld) * form_volume(radius, length);
37    return 1.0e-4 * s * s * form;
38}
39
40
41double Iqxy(double qx, double qy,
42    double sld,
43    double solvent_sld,
44    double radius,
45    double length,
46    double theta,
47    double phi)
48{
49    // TODO: return NaN if radius<0 or length<0?
50    double sn, cn; // slots to hold sincos function output
51
52    // Compute angle alpha between q and the cylinder axis
53    SINCOS(theta*M_PI_180, sn, cn);
54    const double q = sqrt(qx*qx + qy*qy);
55    const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q);
56    const double alpha = acos(cos_val);
57
58    SINCOS(alpha, sn, cn);
59    const double fq = sinc(q*0.5*length*cn) * J1c(q*radius*sn);
60    const double s = (sld-solvent_sld) * form_volume(radius, length);
61    return 1.0e-4 * square(s * fq);
62}
Note: See TracBrowser for help on using the repository browser.