source: sasmodels/sasmodels/models/core_shell_bicelle.c @ 6d5601c

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 6d5601c was 6d5601c, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

use 4/3 pi constant when computing R_eff

  • 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    case 1: // equivalent sphere
58        return radius_from_volume(radius, thick_rim, thick_face, length);
59    case 2: // outer rim radius
60        return radius + thick_rim;
61    case 3: // half outer thickness
62        return 0.5*length + thick_face;
63    case 4: // half diagonal
64        return radius_from_diagonal(radius,thick_rim,thick_face,length);
65    }
66}
67
68static void
69Fq(double q,
70    double *F1,
71    double *F2,
72    double radius,
73    double thick_radius,
74    double thick_face,
75    double length,
76    double sld_core,
77    double sld_face,
78    double sld_rim,
79    double sld_solvent)
80{
81    // set up the integration end points
82    const double uplim = M_PI_4;
83    const double halflength = 0.5*length;
84
85    double total_F1 = 0.0;
86    double total_F2 = 0.0;
87    for(int i=0;i<GAUSS_N;i++) {
88        double theta = (GAUSS_Z[i] + 1.0)*uplim;
89        double sin_theta, cos_theta; // slots to hold sincos function output
90        SINCOS(theta, sin_theta, cos_theta);
91        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
92                                   halflength, sld_core, sld_face, sld_rim, sld_solvent);
93        total_F1 += GAUSS_W[i]*form*sin_theta;
94        total_F2 += GAUSS_W[i]*form*form*sin_theta;
95    }
96    // Correct for integration range
97    total_F1 *= uplim;
98    total_F2 *= uplim;
99
100    *F1 = 1.0e-2*total_F1;
101    *F2 = 1.0e-4*total_F2;
102}
103
104static double
105Iqac(double qab, double qc,
106    double radius,
107    double thick_rim,
108    double thick_face,
109    double length,
110    double core_sld,
111    double face_sld,
112    double rim_sld,
113    double solvent_sld)
114{
115    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
116                           0.5*length, core_sld, face_sld, rim_sld,
117                           solvent_sld);
118    return 1.0e-4*fq*fq;
119}
Note: See TracBrowser for help on using the repository browser.