1 | FourShell(double dp[], double q) |
---|
2 | { |
---|
3 | // variables are: |
---|
4 | //[0] scale factor |
---|
5 | //[1] radius of core [ᅵ] |
---|
6 | //[2] SLD of the core [ᅵ-2] |
---|
7 | //[3] thickness of shell 1 [ᅵ] |
---|
8 | //[4] SLD of shell 1 |
---|
9 | //[5] thickness of shell 2 [ᅵ] |
---|
10 | //[6] SLD of shell 2 |
---|
11 | //[7] thickness of shell 3 |
---|
12 | //[8] SLD of shell 3 |
---|
13 | //[9] thickness of shell 3 |
---|
14 | //[10] SLD of shell 3 |
---|
15 | //[11] SLD of solvent |
---|
16 | //[12] background [cm-1] |
---|
17 | |
---|
18 | double x,pi; |
---|
19 | double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg; //my local names |
---|
20 | double bes,f,vol,qr,contr,f2; |
---|
21 | double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; |
---|
22 | |
---|
23 | pi = 4.0*atan(1.0); |
---|
24 | x=q; |
---|
25 | |
---|
26 | scale = dp[0]; |
---|
27 | rcore = dp[1]; |
---|
28 | rhocore = dp[2]; |
---|
29 | thick1 = dp[3]; |
---|
30 | rhoshel1 = dp[4]; |
---|
31 | thick2 = dp[5]; |
---|
32 | rhoshel2 = dp[6]; |
---|
33 | thick3 = dp[7]; |
---|
34 | rhoshel3 = dp[8]; |
---|
35 | thick4 = dp[9]; |
---|
36 | rhoshel4 = dp[10]; |
---|
37 | rhosolv = dp[11]; |
---|
38 | bkg = dp[12]; |
---|
39 | |
---|
40 | // core first, then add in shells |
---|
41 | qr=x*rcore; |
---|
42 | contr = rhocore-rhoshel1; |
---|
43 | if(qr == 0){ |
---|
44 | bes = 1.0; |
---|
45 | }else{ |
---|
46 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
47 | } |
---|
48 | vol = 4.0*pi/3.0*rcore*rcore*rcore; |
---|
49 | f = vol*bes*contr; |
---|
50 | //now the shell (1) |
---|
51 | qr=x*(rcore+thick1); |
---|
52 | contr = rhoshel1-rhoshel2; |
---|
53 | if(qr == 0){ |
---|
54 | bes = 1.0; |
---|
55 | }else{ |
---|
56 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
57 | } |
---|
58 | vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); |
---|
59 | f += vol*bes*contr; |
---|
60 | //now the shell (2) |
---|
61 | qr=x*(rcore+thick1+thick2); |
---|
62 | contr = rhoshel2-rhoshel3; |
---|
63 | if(qr == 0){ |
---|
64 | bes = 1.0; |
---|
65 | }else{ |
---|
66 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
67 | } |
---|
68 | vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); |
---|
69 | f += vol*bes*contr; |
---|
70 | //now the shell (3) |
---|
71 | qr=x*(rcore+thick1+thick2+thick3); |
---|
72 | contr = rhoshel3-rhoshel4; |
---|
73 | if(qr == 0){ |
---|
74 | bes = 1.0; |
---|
75 | }else{ |
---|
76 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
77 | } |
---|
78 | vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); |
---|
79 | f += vol*bes*contr; |
---|
80 | //now the shell (4) |
---|
81 | qr=x*(rcore+thick1+thick2+thick3+thick4); |
---|
82 | contr = rhoshel4-rhosolv; |
---|
83 | if(qr == 0){ |
---|
84 | bes = 1.0; |
---|
85 | }else{ |
---|
86 | bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); |
---|
87 | } |
---|
88 | vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); |
---|
89 | f += vol*bes*contr; |
---|
90 | |
---|
91 | |
---|
92 | // normalize to particle volume and rescale from [ᅵ-1] to [cm-1] |
---|
93 | f2 = f*f/vol*1.0e8; |
---|
94 | |
---|
95 | //scale if desired |
---|
96 | f2 *= scale; |
---|
97 | // then add in the background |
---|
98 | f2 += bkg; |
---|
99 | |
---|
100 | return(f2); |
---|
101 | } |
---|
102 | |
---|
103 | // above is cut and past with no changes from libigor four shell |
---|
104 | static double |
---|
105 | f_exp(double q, double r, double sld_in, double sld_out, |
---|
106 | double thickness, double A) |
---|
107 | { |
---|
108 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
109 | const double qr = q * r; |
---|
110 | const double alpha = A * r/thickness; |
---|
111 | const double bes = sph_j1c(qr); |
---|
112 | const double B = (sld_out - sld_in)/expm1(A); |
---|
113 | const double C = sld_in - B; |
---|
114 | double fun; |
---|
115 | if (qr == 0.0) { |
---|
116 | fun = 1.0; |
---|
117 | } else if (fabs(A) > 0.0) { |
---|
118 | const double qrsq = qr*qr; |
---|
119 | const double alphasq = alpha*alpha; |
---|
120 | const double sumsq = alphasq + qrsq; |
---|
121 | double sinqr, cosqr; |
---|
122 | SINCOS(qr, sinqr, cosqr); |
---|
123 | fun = -3.0( |
---|
124 | ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq |
---|
125 | - (alpha*sinqr/qr - cosqr) |
---|
126 | ) / sumsq; |
---|
127 | } else { |
---|
128 | fun = bes; |
---|
129 | } |
---|
130 | return vol * (B*fun + C*bes); |
---|
131 | } |
---|
132 | |
---|
133 | static double |
---|
134 | f_linear(double q, double r, double sld, double slope) |
---|
135 | { |
---|
136 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
137 | const double qr = q * r; |
---|
138 | const double bes = sph_j1c(qr); |
---|
139 | double fun = 0.0; |
---|
140 | if (qr > 0.0) { |
---|
141 | const double qrsq = qr*qr; |
---|
142 | double sinqr, cosqr; |
---|
143 | SINCOS(qr, sinqr, cosqr); |
---|
144 | // Jae-He's code seems to simplify to this |
---|
145 | // fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); |
---|
146 | // Rederiving the math, we get the following instead: |
---|
147 | fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); |
---|
148 | } |
---|
149 | return vol * (sld*bes + fun); |
---|
150 | } |
---|
151 | |
---|
152 | static double |
---|
153 | f_constant(double q, double r, double sld) |
---|
154 | { |
---|
155 | const double bes = sph_j1c(q * r); |
---|
156 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
157 | return sld * vol * bes; |
---|
158 | } |
---|
159 | |
---|
160 | static double |
---|
161 | form_volume(double core_radius, double n, double thickness[]) |
---|
162 | { |
---|
163 | int i; |
---|
164 | double r = core_radius; |
---|
165 | for (i=0; i < n; i++) { |
---|
166 | r += thickness[i]; |
---|
167 | } |
---|
168 | return 4.0/3.0 * M_PI * r * r * r; |
---|
169 | } |
---|
170 | |
---|
171 | static double |
---|
172 | Iq(double q, double core_sld, double core_radius, double solvent_sld, |
---|
173 | double n, double in_sld[], double out_sld[], double thickness[], |
---|
174 | double A[]) |
---|
175 | { |
---|
176 | int i; |
---|
177 | r = core_radius; |
---|
178 | f = f_constant(q, r, core_sld); |
---|
179 | for (i=0; i<n; i++){ |
---|
180 | const double r0 = r; |
---|
181 | r += thickness[i]; |
---|
182 | if (r == r0) { |
---|
183 | // no thickness, so nothing to add |
---|
184 | } else if (fabs(A[i]) < 1e-16 || sld_out[i] == sld_in[i]) { |
---|
185 | f -= f_constant(q, r0, sld_in[i]); |
---|
186 | f += f_constant(q, r, sld_in[i]); |
---|
187 | } else if (fabs(A[i]) < 1e-4) { |
---|
188 | const double slope = (sld_out[i] - sld_in[i])/thickness[i]; |
---|
189 | f -= f_linear(q, r0, sld_in[i], slope); |
---|
190 | f += f_linear(q, r, sld_out[i], slope); |
---|
191 | } else { |
---|
192 | f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
193 | f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
194 | } |
---|
195 | } |
---|
196 | f -= f_constant(q, r, solvent_sld); |
---|
197 | f2 = f * f * 1.0e-4; |
---|
198 | |
---|
199 | return f2; |
---|
200 | } |
---|