1 | static double effective_radius(double radius_core, double epsilon, double alpha) |
---|
2 | { |
---|
3 | return radius_core |
---|
4 | * sqrt(square(sin(alpha)) + square(epsilon * cos(alpha))); |
---|
5 | } |
---|
6 | |
---|
7 | static double bes_core(double q, double or_radius) |
---|
8 | { |
---|
9 | return sph_j1c(q*or_radius); |
---|
10 | } |
---|
11 | |
---|
12 | static double bes_corona(double q, double or_radius, double radius_gyr, double d_penetration) |
---|
13 | { |
---|
14 | return sinc(q * (or_radius + d_penetration * radius_gyr)); |
---|
15 | } |
---|
16 | |
---|
17 | static double Fs(double q, double radius_core, double epsilon) |
---|
18 | { |
---|
19 | const int bins = 100; |
---|
20 | const double dalpha = M_PI_2 / bins; |
---|
21 | double total = 0.0; |
---|
22 | for(int i = 0; i != bins; ++i) { |
---|
23 | const double alpha = i * dalpha; |
---|
24 | const double or_radius=effective_radius(radius_core,epsilon,alpha); |
---|
25 | total += square(bes_core(q, or_radius)) * sin(alpha); |
---|
26 | } |
---|
27 | return total * dalpha; |
---|
28 | } |
---|
29 | |
---|
30 | static double Ssc(double q, double radius_core, double radius_gyr, double d_penetration, double epsilon, double chain_ampl) |
---|
31 | { |
---|
32 | const int bins = 100; |
---|
33 | const double dalpha = M_PI_2 / bins; |
---|
34 | double total = 0.0; |
---|
35 | for(int i = 0; i != bins; ++i) { |
---|
36 | const double alpha = i * dalpha; |
---|
37 | const double or_radius=effective_radius(radius_core, epsilon, alpha); |
---|
38 | total += bes_core(q, or_radius) * bes_corona(q, or_radius, radius_gyr, d_penetration) * sin(alpha); |
---|
39 | } |
---|
40 | return total * dalpha * chain_ampl; |
---|
41 | } |
---|
42 | |
---|
43 | static double Scc(double q, double radius_core, double radius_gyr, double d_penetration, double epsilon, double chain_ampl) |
---|
44 | { |
---|
45 | const int bins = 100; |
---|
46 | const double dalpha = M_PI_2 / bins; |
---|
47 | double total = 0.0; |
---|
48 | for(int i = 0; i != bins; ++i) { |
---|
49 | const double alpha = i * dalpha; |
---|
50 | const double or_radius=effective_radius(radius_core,epsilon,alpha); |
---|
51 | total += square(bes_corona(q, or_radius, radius_gyr, d_penetration)) * sin(alpha); |
---|
52 | } |
---|
53 | return total * dalpha * chain_ampl * chain_ampl; |
---|
54 | } |
---|
55 | |
---|
56 | static double Iq( |
---|
57 | double q, |
---|
58 | double ndensity, |
---|
59 | double v_core, |
---|
60 | double v_corona, |
---|
61 | double solvent_sld, |
---|
62 | double core_sld, |
---|
63 | double corona_sld, |
---|
64 | double radius_core, |
---|
65 | double radius_gyr, |
---|
66 | double d_penetration, |
---|
67 | double n_aggreg, |
---|
68 | double epsilon) |
---|
69 | { |
---|
70 | const double beta_core = v_core * (core_sld - solvent_sld); |
---|
71 | const double beta_corona = v_corona * (corona_sld - solvent_sld); |
---|
72 | |
---|
73 | const double term1 = square(n_aggreg*beta_core)*Fs(q, radius_core, epsilon); |
---|
74 | |
---|
75 | // Self-correlation term of the chains |
---|
76 | const double qrg2 = square(q*radius_gyr); |
---|
77 | const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2); |
---|
78 | const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain; |
---|
79 | |
---|
80 | |
---|
81 | const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; |
---|
82 | |
---|
83 | // Interference cross-term between core and chains |
---|
84 | const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona * Ssc(q, radius_core, radius_gyr, d_penetration, epsilon, chain_ampl); |
---|
85 | |
---|
86 | // Interference cross-term between chains |
---|
87 | const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona * Scc(q, radius_core, radius_gyr, d_penetration, epsilon, chain_ampl); |
---|
88 | |
---|
89 | // I(q)_micelle : Sum of 4 terms computed above |
---|
90 | double i_micelle = term1 + term2 + term3 + term4; |
---|
91 | |
---|
92 | // rescale from [A^2] to [cm^2] |
---|
93 | i_micelle *= 1.0e-13; |
---|
94 | |
---|
95 | // "normalize" by number density --> intensity in [cm-1] |
---|
96 | i_micelle *= ndensity; |
---|
97 | |
---|
98 | return(i_micelle); |
---|
99 | } |
---|
100 | |
---|