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
131
132
133
134
135
|
static double
form_volume(double length_a, double length_b, double length_c)
{
return length_a * length_b * length_c;
}
static double
radius_from_excluded_volume(double length_a, double length_b, double length_c)
{
double r_equiv, length;
double lengths[3] = {length_a, length_b, length_c};
double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2]));
double length_1 = lengthmax;
double length_2 = lengthmax;
double length_3 = lengthmax;
for(int ilen=0; ilen<3; ilen++) {
if (lengths[ilen] < length_1) {
length_2 = length_1;
length_1 = lengths[ilen];
} else {
if (lengths[ilen] < length_2) {
length_2 = lengths[ilen];
}
}
}
if(length_2-length_1 > length_3-length_2) {
r_equiv = sqrt(length_2*length_3/M_PI);
length = length_1;
} else {
r_equiv = sqrt(length_1*length_2/M_PI);
length = length_3;
}
return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
}
static double
radius_effective(int mode, double length_a, double length_b, double length_c)
{
switch (mode) {
default:
case 1: // equivalent cylinder excluded volume
return radius_from_excluded_volume(length_a,length_b,length_c);
case 2: // equivalent volume sphere
return cbrt(length_a*length_b*length_c/M_4PI_3);
case 3: // half length_a
return 0.5 * length_a;
case 4: // half length_b
return 0.5 * length_b;
case 5: // half length_c
return 0.5 * length_c;
case 6: // equivalent circular cross-section
return sqrt(length_a*length_b/M_PI);
case 7: // half ab diagonal
return 0.5*sqrt(length_a*length_a + length_b*length_b);
case 8: // half diagonal
return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);
}
}
static void
Fq(double q,
double *F1,
double *F2,
double sld,
double solvent_sld,
double length_a,
double length_b,
double length_c)
{
const double mu = 0.5 * q * length_b;
// Scale sides by B
const double a_scaled = length_a / length_b;
const double c_scaled = length_c / length_b;
// outer integral (with gauss points), integration limits = 0, 1
double outer_total_F1 = 0.0; //initialize integral
double outer_total_F2 = 0.0; //initialize integral
for( int i=0; i<GAUSS_N; i++) {
const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 );
const double mu_proj = mu * sqrt(1.0-sigma*sigma);
// inner integral (with gauss points), integration limits = 0, 1
// corresponding to angles from 0 to pi/2.
double inner_total_F1 = 0.0;
double inner_total_F2 = 0.0;
for(int j=0; j<GAUSS_N; j++) {
const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 );
double sin_uu, cos_uu;
SINCOS(M_PI_2*uu, sin_uu, cos_uu);
const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
const double si2 = sas_sinx_x(mu_proj * cos_uu);
const double fq = si1 * si2;
inner_total_F1 += GAUSS_W[j] * fq;
inner_total_F2 += GAUSS_W[j] * fq * fq;
}
// now complete change of inner integration variable (1-0)/(1-(-1))= 0.5
inner_total_F1 *= 0.5;
inner_total_F2 *= 0.5;
const double si = sas_sinx_x(mu * c_scaled * sigma);
outer_total_F1 += GAUSS_W[i] * inner_total_F1 * si;
outer_total_F2 += GAUSS_W[i] * inner_total_F2 * si * si;
}
// now complete change of outer integration variable (1-0)/(1-(-1))= 0.5
outer_total_F1 *= 0.5;
outer_total_F2 *= 0.5;
// Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1]
const double V = form_volume(length_a, length_b, length_c);
const double contrast = (sld-solvent_sld);
const double s = contrast * V;
*F1 = 1.0e-2 * s * outer_total_F1;
*F2 = 1.0e-4 * s * s * outer_total_F2;
}
static double
Iqabc(double qa, double qb, double qc,
double sld,
double solvent_sld,
double length_a,
double length_b,
double length_c)
{
const double siA = sas_sinx_x(0.5*length_a*qa);
const double siB = sas_sinx_x(0.5*length_b*qb);
const double siC = sas_sinx_x(0.5*length_c*qc);
const double V = form_volume(length_a, length_b, length_c);
const double drho = (sld - solvent_sld);
const double form = V * drho * siA * siB * siC;
// Square and convert from [1e-12 A-1] to [cm-1]
return 1.0e-4 * form * form;
}
|
