source: sasmodels/sasmodels/models/core_shell_cylinder.c @ 99658f6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 99658f6 was 99658f6, checked in by grethevj, 10 months ago

updated ER functions including cylinder excluded volume, to match 4.x

  • Property mode set to 100644
File size: 4.0 KB
Line 
1// vd = volume * delta_rho
2// besarg = q * R * sin(theta)
3// siarg = q * L/2 * cos(theta)
4static double _cyl(double vd, double besarg, double siarg)
5{
6    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
7}
8
9static double
10form_volume(double radius, double thickness, double length)
11{
12    return M_PI*square(radius+thickness)*(length+2.0*thickness);
13}
14
15static double
16radius_from_excluded_volume(double radius, double thickness, double length)
17{
18    const double radius_tot = radius + thickness;
19    const double length_tot = length + 2.0*thickness;
20    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot)));
21}
22
23static double
24radius_from_volume(double radius, double thickness, double length)
25{
26    const double volume_outer_cyl = form_volume(radius,thickness,length);
27    return cbrt(volume_outer_cyl/M_4PI_3);
28}
29
30static double
31radius_from_diagonal(double radius, double thickness, double length)
32{
33    const double radius_outer = radius + thickness;
34    const double length_outer = length + 2.0*thickness;
35    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer);
36}
37
38static double
39effective_radius(int mode, double radius, double thickness, double length)
40{
41    switch (mode) {
42    default:
43    case 1: //cylinder excluded volume
44        return radius_from_excluded_volume(radius, thickness, length);
45    case 2: // equivalent volume sphere
46        return radius_from_volume(radius, thickness, length);
47    case 3: // outer radius
48        return radius + thickness;
49    case 4: // half outer length
50        return 0.5*length + thickness;
51    case 5: // half min outer length
52        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);
53    case 6: // half max outer length
54        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);
55    case 7: // half outer diagonal
56        return radius_from_diagonal(radius,thickness,length);
57    }
58}
59
60static void
61Fq(double q,
62    double *F1,
63    double *F2,
64    double core_sld,
65    double shell_sld,
66    double solvent_sld,
67    double radius,
68    double thickness,
69    double length)
70{
71    // precalculate constants
72    const double core_r = radius;
73    const double core_h = 0.5*length;
74    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
75    const double shell_r = (radius + thickness);
76    const double shell_h = (0.5*length + thickness);
77    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
78    double total_F1 = 0.0;
79    double total_F2 = 0.0;
80    for (int i=0; i<GAUSS_N ;i++) {
81        // translate a point in [-1,1] to a point in [0, pi/2]
82        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0;
83        double sin_theta, cos_theta;
84        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4;
85        SINCOS(theta, sin_theta,  cos_theta);
86        const double qab = q*sin_theta;
87        const double qc = q*cos_theta;
88        const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
89            + _cyl(shell_vd, shell_r*qab, shell_h*qc);
90        total_F1 += GAUSS_W[i] * fq * sin_theta;
91        total_F2 += GAUSS_W[i] * fq * fq * sin_theta;
92    }
93    // translate dx in [-1,1] to dx in [lower,upper]
94    //const double form = (upper-lower)/2.0*total;
95    *F1 = 1.0e-2 * total_F1 * M_PI_4;
96    *F2 = 1.0e-4 * total_F2 * M_PI_4;
97}
98
99static double
100Iqac(double qab, double qc,
101    double core_sld,
102    double shell_sld,
103    double solvent_sld,
104    double radius,
105    double thickness,
106    double length)
107{
108    const double core_r = radius;
109    const double core_h = 0.5*length;
110    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
111    const double shell_r = (radius + thickness);
112    const double shell_h = (0.5*length + thickness);
113    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
114
115    const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
116        + _cyl(shell_vd, shell_r*qab, shell_h*qc);
117    return 1.0e-4 * fq * fq;
118}
Note: See TracBrowser for help on using the repository browser.