source: sasmodels/sasmodels/models/core_shell_bicelle.c @ d277229

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d277229 was d277229, checked in by grethevj, 6 years ago

Models updated to include choices for effective interaction radii

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