1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
static double
form_volume(double radius, double thick_rim, double thick_face, double length)
{
return M_PI*square(radius+thick_rim)*(length+2.0*thick_face);
}
static double
bicelle_kernel(double qab,
double qc,
double radius,
double thick_radius,
double thick_face,
double halflength,
double sld_core,
double sld_face,
double sld_rim,
double sld_solvent)
{
const double dr1 = sld_core-sld_face;
const double dr2 = sld_rim-sld_solvent;
const double dr3 = sld_face-sld_rim;
const double vol1 = M_PI*square(radius)*2.0*(halflength);
const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face);
const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face);
const double be1 = sas_2J1x_x((radius)*qab);
const double be2 = sas_2J1x_x((radius+thick_radius)*qab);
const double si1 = sas_sinx_x((halflength)*qc);
const double si2 = sas_sinx_x((halflength+thick_face)*qc);
const double t = vol1*dr1*si1*be1 +
vol2*dr2*si2*be2 +
vol3*dr3*si2*be1;
return t;
}
static double
radius_from_excluded_volume(double radius, double thick_rim, double thick_face, double length)
{
const double radius_tot = radius + thick_rim;
const double length_tot = length + 2.0*thick_face;
return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot)));
}
static double
radius_from_volume(double radius, double thick_rim, double thick_face, double length)
{
const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length);
return cbrt(volume_bicelle/M_4PI_3);
}
static double
radius_from_diagonal(double radius, double thick_rim, double thick_face, double length)
{
const double radius_tot = radius + thick_rim;
const double length_tot = length + 2.0*thick_face;
return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot);
}
static double
radius_effective(int mode, double radius, double thick_rim, double thick_face, double length)
{
switch (mode) {
default:
case 1: // equivalent cylinder excluded volume
return radius_from_excluded_volume(radius, thick_rim, thick_face, length);
case 2: // equivalent sphere
return radius_from_volume(radius, thick_rim, thick_face, length);
case 3: // outer rim radius
return radius + thick_rim;
case 4: // half outer thickness
return 0.5*length + thick_face;
case 5: // half diagonal
return radius_from_diagonal(radius,thick_rim,thick_face,length);
}
}
static void
Fq(double q,
double *F1,
double *F2,
double radius,
double thick_radius,
double thick_face,
double length,
double sld_core,
double sld_face,
double sld_rim,
double sld_solvent)
{
// set up the integration end points
const double uplim = M_PI_4;
const double halflength = 0.5*length;
double total_F1 = 0.0;
double total_F2 = 0.0;
for(int i=0;i<GAUSS_N;i++) {
double theta = (GAUSS_Z[i] + 1.0)*uplim;
double sin_theta, cos_theta; // slots to hold sincos function output
SINCOS(theta, sin_theta, cos_theta);
double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face,
halflength, sld_core, sld_face, sld_rim, sld_solvent);
total_F1 += GAUSS_W[i]*form*sin_theta;
total_F2 += GAUSS_W[i]*form*form*sin_theta;
}
// Correct for integration range
total_F1 *= uplim;
total_F2 *= uplim;
*F1 = 1.0e-2*total_F1;
*F2 = 1.0e-4*total_F2;
}
static double
Iqac(double qab, double qc,
double radius,
double thick_rim,
double thick_face,
double length,
double core_sld,
double face_sld,
double rim_sld,
double solvent_sld)
{
double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face,
0.5*length, core_sld, face_sld, rim_sld,
solvent_sld);
return 1.0e-4*fq*fq;
}
|
