source: sasmodels/sasmodels/models/core_shell_cylinder.c @ 0d6e865

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0d6e865 was 0d6e865, checked in by dirk, 7 years ago

Rewriting selected models in spherical coordinates. This fixes ticket #492 for these models.

  • Property mode set to 100644
File size: 3.4 KB
Line 
1double form_volume(double radius, double thickness, double length);
2double Iq(double q, double core_sld, double shell_sld, double solvent_sld,
3    double radius, double thickness, double length);
4double 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
7// twovd = 2 * volume * delta_rho
8// besarg = q * R * sin(alpha)
9// siarg = q * L/2 * cos(alpha)
10double _cyl(double twovd, double besarg, double siarg);
11double _cyl(double twovd, double besarg, double siarg)
12{
13    const double bj = (besarg == 0.0 ? 0.5 : 0.5*sas_J1c(besarg));
14    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
15    return twovd*si*bj;
16}
17
18double form_volume(double radius, double thickness, double length)
19{
20    return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness);
21}
22
23double Iq(double q,
24    double core_sld,
25    double shell_sld,
26    double solvent_sld,
27    double radius,
28    double thickness,
29    double length)
30{
31    // precalculate constants
32    const double core_qr = q*radius;
33    const double core_qh = q*0.5*length;
34    const double core_twovd = 2.0 * form_volume(radius,0,length)
35                            * (core_sld-shell_sld);
36    const double shell_qr = q*(radius + thickness);
37    const double shell_qh = q*(0.5*length + thickness);
38    const double shell_twovd = 2.0 * form_volume(radius,thickness,length)
39                             * (shell_sld-solvent_sld);
40    double total = 0.0;
41    // double lower=0, upper=M_PI_2;
42    for (int i=0; i<76 ;i++) {
43        // translate a point in [-1,1] to a point in [lower,upper]
44        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
45        double sn, cn;
46        const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2);
47        SINCOS(alpha, sn, cn);
48        const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn)
49            + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn);
50        total += Gauss76Wt[i] * fq * fq * sn;
51    }
52    // translate dx in [-1,1] to dx in [lower,upper]
53    //const double form = (upper-lower)/2.0*total;
54    return 1.0e-4 * total * M_PI_4;
55}
56
57
58double Iqxy(double qx, double qy,
59    double core_sld,
60    double shell_sld,
61    double solvent_sld,
62    double radius,
63    double thickness,
64    double length,
65    double theta,
66    double phi)
67{
68    double sn, cn; // slots to hold sincos function output
69
70    // Compute angle alpha between q and the cylinder axis
71    SINCOS(phi*M_PI_180, sn, cn);
72    // # The following correction factor exists in sasview, but it can't be
73    // # right, so we are leaving it out for now.
74    // const double correction = fabs(cn)*M_PI_2;
75    const double q = sqrt(qx*qx+qy*qy);
76    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
77    const double alpha = acos(cos_val);
78
79    const double core_qr = q*radius;
80    const double core_qh = q*0.5*length;
81    const double core_twovd = 2.0 * form_volume(radius,0,length)
82                            * (core_sld-shell_sld);
83    const double shell_qr = q*(radius + thickness);
84    const double shell_qh = q*(0.5*length + thickness);
85    const double shell_twovd = 2.0 * form_volume(radius,thickness,length)
86                             * (shell_sld-solvent_sld);
87
88    SINCOS(alpha, sn, cn);
89    const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn)
90        + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn);
91    return 1.0e-4 * fq * fq;
92}
Note: See TracBrowser for help on using the repository browser.