source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 5024a56

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 5024a56 was d42dd4a, checked in by pkienzle, 5 years ago

fix compiler warnings for CUDA

  • Property mode set to 100644
File size: 3.5 KB
Line 
1static double
2form_volume(double radius, double thick_rim, double thick_face, double length)
3{
4    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face);
5}
6
7static double
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)
18{
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);
25
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);
30
31    const double t = vol1*dr1*si1*be1 +
32                     vol2*dr2*si2*be2 +
33                     vol3*dr3*si2*be1;
34
35    return t;
36}
37
38static double
39radius_from_volume(double radius, double thick_rim, double thick_face, double length)
40{
41    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length);
42    return cbrt(volume_bicelle/M_4PI_3);
43}
44
45static double
46radius_from_diagonal(double radius, double thick_rim, double thick_face, double length)
47{
48    const double radius_tot = radius + thick_rim;
49    const double length_tot = length + 2.0*thick_face;
50    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot);
51}
52
53static double
54effective_radius(int mode, double radius, double thick_rim, double thick_face, double length)
55{
56    switch (mode) {
57    default:
58    case 1: // equivalent sphere
59        return radius_from_volume(radius, thick_rim, thick_face, length);
60    case 2: // outer rim radius
61        return radius + thick_rim;
62    case 3: // half outer thickness
63        return 0.5*length + thick_face;
64    case 4: // half diagonal
65        return radius_from_diagonal(radius,thick_rim,thick_face,length);
66    }
67}
68
69static void
70Fq(double q,
71    double *F1,
72    double *F2,
73    double radius,
74    double thick_radius,
75    double thick_face,
76    double length,
77    double sld_core,
78    double sld_face,
79    double sld_rim,
80    double sld_solvent)
81{
82    // set up the integration end points
83    const double uplim = M_PI_4;
84    const double halflength = 0.5*length;
85
86    double total_F1 = 0.0;
87    double total_F2 = 0.0;
88    for(int i=0;i<GAUSS_N;i++) {
89        double theta = (GAUSS_Z[i] + 1.0)*uplim;
90        double sin_theta, cos_theta; // slots to hold sincos function output
91        SINCOS(theta, sin_theta, cos_theta);
92        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
93                                   halflength, sld_core, sld_face, sld_rim, sld_solvent);
94        total_F1 += GAUSS_W[i]*form*sin_theta;
95        total_F2 += GAUSS_W[i]*form*form*sin_theta;
96    }
97    // Correct for integration range
98    total_F1 *= uplim;
99    total_F2 *= uplim;
100
101    *F1 = 1.0e-2*total_F1;
102    *F2 = 1.0e-4*total_F2;
103}
104
105static double
106Iqac(double qab, double qc,
107    double radius,
108    double thick_rim,
109    double thick_face,
110    double length,
111    double core_sld,
112    double face_sld,
113    double rim_sld,
114    double solvent_sld)
115{
116    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
117                           0.5*length, core_sld, face_sld, rim_sld,
118                           solvent_sld);
119    return 1.0e-4*fq*fq;
120}
Note: See TracBrowser for help on using the repository browser.