source: sasmodels/sasmodels/models/core_shell_bicelle.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 
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_excluded_volume(double radius, double thick_rim, double thick_face, double length)
40{
41    const double radius_tot = radius + thick_rim;
42    const double length_tot = length + 2.0*thick_face;
43    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot)));
44}
45
46static double
47radius_from_volume(double radius, double thick_rim, double thick_face, double length)
48{
49    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length);
50    return cbrt(volume_bicelle/M_4PI_3);
51}
52
53static double
54radius_from_diagonal(double radius, double thick_rim, double thick_face, double length)
55{
56    const double radius_tot = radius + thick_rim;
57    const double length_tot = length + 2.0*thick_face;
58    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot);
59}
60
61static double
62effective_radius(int mode, double radius, double thick_rim, double thick_face, double length)
63{
64    switch (mode) {
65    default:
66    case 1: // equivalent cylinder excluded volume
67        return radius_from_excluded_volume(radius, thick_rim, thick_face, length);
68    case 2: // equivalent sphere
69        return radius_from_volume(radius, thick_rim, thick_face, length);
70    case 3: // outer rim radius
71        return radius + thick_rim;
72    case 4: // half outer thickness
73        return 0.5*length + thick_face;
74    case 5: // half diagonal
75        return radius_from_diagonal(radius,thick_rim,thick_face,length);
76    }
77}
78
79static void
80Fq(double q,
81    double *F1,
82    double *F2,
83    double radius,
84    double thick_radius,
85    double thick_face,
86    double length,
87    double sld_core,
88    double sld_face,
89    double sld_rim,
90    double sld_solvent)
91{
92    // set up the integration end points
93    const double uplim = M_PI_4;
94    const double halflength = 0.5*length;
95
96    double total_F1 = 0.0;
97    double total_F2 = 0.0;
98    for(int i=0;i<GAUSS_N;i++) {
99        double theta = (GAUSS_Z[i] + 1.0)*uplim;
100        double sin_theta, cos_theta; // slots to hold sincos function output
101        SINCOS(theta, sin_theta, cos_theta);
102        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
103                                   halflength, sld_core, sld_face, sld_rim, sld_solvent);
104        total_F1 += GAUSS_W[i]*form*sin_theta;
105        total_F2 += GAUSS_W[i]*form*form*sin_theta;
106    }
107    // Correct for integration range
108    total_F1 *= uplim;
109    total_F2 *= uplim;
110
111    *F1 = 1.0e-2*total_F1;
112    *F2 = 1.0e-4*total_F2;
113}
114
115static double
116Iqac(double qab, double qc,
117    double radius,
118    double thick_rim,
119    double thick_face,
120    double length,
121    double core_sld,
122    double face_sld,
123    double rim_sld,
124    double solvent_sld)
125{
126    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
127                           0.5*length, core_sld, face_sld, rim_sld,
128                           solvent_sld);
129    return 1.0e-4*fq*fq;
130}
Note: See TracBrowser for help on using the repository browser.