source: sasmodels/sasmodels/models/cylinder_clone.c @ 994d77f

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

Convert double to float rather than using real

  • Property mode set to 100644
File size: 2.5 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, double radius, double length, double theta, double phi);
4
5
6// twovd = 2 * volume * delta_rho
7// besarg = q * R * sin(alpha)
8// siarg = q * L/2 * cos(alpha)
9double _cyl(double twovd, double besarg, double siarg, double alpha);
10double _cyl(double twovd, double besarg, double siarg, double alpha)
11{
12    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
13    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
14    return twovd*si*bj;
15}
16
17double form_volume(double radius, double length)
18{
19    return M_PI*radius*radius*length;
20}
21double Iq(double q,
22    double sldCyl,
23    double sldSolv,
24    double radius,
25    double length)
26{
27    const double qr = q*radius;
28    const double qh = q*0.5*length;
29    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length);
30    double total = 0.0;
31    // double lower=0, upper=M_PI_2;
32    for (int i=0; i<76 ;i++) {
33        // translate a point in [-1,1] to a point in [lower,upper]
34        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
35        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
36        double sn, cn;
37        SINCOS(alpha, sn, cn);
38        const double fq = _cyl(twovd, qr*sn, qh*cn, alpha);
39        total += Gauss76Wt[i] * fq * fq * sn;
40    }
41    // translate dx in [-1,1] to dx in [lower,upper]
42    //const double form = (upper-lower)/2.0*total;
43    return 1.0e8 * total * M_PI_4;
44}
45
46double Iqxy(double qx, double qy,
47    double sldCyl,
48    double sldSolv,
49    double radius,
50    double length,
51    double cyl_theta,
52    double cyl_phi)
53{
54    double sn, cn; // slots to hold sincos function output
55
56    // Compute angle alpha between q and the cylinder axis
57    SINCOS(cyl_theta*M_PI_180, sn, cn);
58    // # The following correction factor exists in sasview, but it can't be
59    // # right, so we are leaving it out for now.
60    const double spherical_integration = fabs(cn)*M_PI_2;
61    const double q = sqrt(qx*qx+qy*qy);
62    const double cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q);
63    const double alpha = acos(cos_val);
64
65    const double qr = q*radius;
66    const double qh = q*0.5*length;
67    const double twovd = 2.0*(sldCyl-sldSolv)*form_volume(radius, length);
68    SINCOS(alpha, sn, cn);
69    const double fq = _cyl(twovd, qr*sn, qh*cn, alpha);
70    return 1.0e8 * fq * fq * spherical_integration;
71}
Note: See TracBrowser for help on using the repository browser.