1 | double Iq(double q, double case_num, |
---|
2 | double N[], double Phi[], double v[], double L[], double b[], |
---|
3 | double Kab, double Kac, double Kad, |
---|
4 | double Kbc, double Kbd, double Kcd |
---|
5 | ); |
---|
6 | |
---|
7 | double Iq(double q, double case_num, |
---|
8 | double N[], double Phi[], double v[], double L[], double b[], |
---|
9 | double Kab, double Kac, double Kad, |
---|
10 | double Kbc, double Kbd, double Kcd |
---|
11 | ) { |
---|
12 | int icase = (int)case_num; |
---|
13 | |
---|
14 | #if 0 // Sasview defaults |
---|
15 | if (icase <= 1) { |
---|
16 | N[0]=N[1]=1000.0; |
---|
17 | Phi[0]=Phi[1]=0.0000001; |
---|
18 | Kab=Kac=Kad=Kbc=Kbd=-0.0004; |
---|
19 | La=Lb=1.0e-12; |
---|
20 | va=vb=100.0; |
---|
21 | ba=bb=5.0; |
---|
22 | } else if (icase <= 4) { |
---|
23 | Phi[0]=0.0000001; |
---|
24 | Kab=Kac=Kad=-0.0004; |
---|
25 | La=1.0e-12; |
---|
26 | va=100.0; |
---|
27 | ba=5.0; |
---|
28 | } |
---|
29 | #else |
---|
30 | if (icase <= 1) { |
---|
31 | N[0]=N[1]=0.0; |
---|
32 | Phi[0]=Phi[1]=0.0; |
---|
33 | Kab=Kac=Kad=Kbc=Kbd=0.0; |
---|
34 | L[0]=L[1]=L[3]; |
---|
35 | v[0]=v[1]=v[3]; |
---|
36 | b[0]=b[1]=0.0; |
---|
37 | } else if (icase <= 4) { |
---|
38 | N[0] = 0.0; |
---|
39 | Phi[0]=0.0; |
---|
40 | Kab=Kac=Kad=0.0; |
---|
41 | L[0]=L[3]; |
---|
42 | v[0]=v[3]; |
---|
43 | b[0]=0.0; |
---|
44 | } |
---|
45 | #endif |
---|
46 | |
---|
47 | const double Xa = q*q*b[0]*b[0]*N[0]/6.0; |
---|
48 | const double Xb = q*q*b[1]*b[1]*N[1]/6.0; |
---|
49 | const double Xc = q*q*b[2]*b[2]*N[2]/6.0; |
---|
50 | const double Xd = q*q*b[3]*b[3]*N[3]/6.0; |
---|
51 | |
---|
52 | // limit as Xa goes to 0 is 1 |
---|
53 | const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa; |
---|
54 | const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb; |
---|
55 | const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc; |
---|
56 | const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd; |
---|
57 | |
---|
58 | // limit as Xa goes to 0 is 1 |
---|
59 | const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa; |
---|
60 | const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb; |
---|
61 | const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc; |
---|
62 | const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd; |
---|
63 | |
---|
64 | |
---|
65 | // Note: S0ij only defined for copolymers; otherwise set to zero |
---|
66 | // 0: C/D binary mixture |
---|
67 | // 1: C-D diblock copolymer |
---|
68 | // 2: B/C/D ternery mixture |
---|
69 | // 3: B/C-D binary mixture,1 homopolymer, 1 diblock copolymer |
---|
70 | // 4: B-C-D triblock copolymer |
---|
71 | // 5: A/B/C/D quaternary mixture |
---|
72 | // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer |
---|
73 | // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer |
---|
74 | // 8: A-B/C-D binary mixture, 2 diblock copolymer |
---|
75 | // 9: A-B-C-D tetra-block copolymer |
---|
76 | #if 0 |
---|
77 | const double S0aa = icase<5 |
---|
78 | ? 1.0 : N[0]*Phi[0]*v[0]*Paa; |
---|
79 | const double S0bb = icase<2 |
---|
80 | ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; |
---|
81 | const double S0cc = N[2]*Phi[2]*v[2]*Pcc; |
---|
82 | const double S0dd = N[3]*Phi[3]*v[3]*Pdd; |
---|
83 | const double S0ab = icase<8 |
---|
84 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; |
---|
85 | const double S0ac = icase<9 |
---|
86 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); |
---|
87 | const double S0ad = icase<9 |
---|
88 | ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); |
---|
89 | const double S0bc = (icase!=4 && icase!=7 && icase!= 9) |
---|
90 | ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; |
---|
91 | const double S0bd = (icase!=4 && icase!=7 && icase!= 9) |
---|
92 | ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); |
---|
93 | const double S0cd = (icase==0 || icase==2 || icase==5) |
---|
94 | ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; |
---|
95 | #else // sasview equivalent |
---|
96 | //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); |
---|
97 | double S0aa = N[0]*Phi[0]*v[0]*Paa; |
---|
98 | double S0bb = N[1]*Phi[1]*v[1]*Pbb; |
---|
99 | double S0cc = N[2]*Phi[2]*v[2]*Pcc; |
---|
100 | double S0dd = N[3]*Phi[3]*v[3]*Pdd; |
---|
101 | double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; |
---|
102 | double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); |
---|
103 | double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); |
---|
104 | double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; |
---|
105 | double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); |
---|
106 | double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; |
---|
107 | switch(icase){ |
---|
108 | case 0: |
---|
109 | S0aa=0.000001; |
---|
110 | S0ab=0.000002; |
---|
111 | S0ac=0.000003; |
---|
112 | S0ad=0.000004; |
---|
113 | S0bb=0.000005; |
---|
114 | S0bc=0.000006; |
---|
115 | S0bd=0.000007; |
---|
116 | S0cd=0.000008; |
---|
117 | break; |
---|
118 | case 1: |
---|
119 | S0aa=0.000001; |
---|
120 | S0ab=0.000002; |
---|
121 | S0ac=0.000003; |
---|
122 | S0ad=0.000004; |
---|
123 | S0bb=0.000005; |
---|
124 | S0bc=0.000006; |
---|
125 | S0bd=0.000007; |
---|
126 | break; |
---|
127 | case 2: |
---|
128 | S0aa=0.000001; |
---|
129 | S0ab=0.000002; |
---|
130 | S0ac=0.000003; |
---|
131 | S0ad=0.000004; |
---|
132 | S0bc=0.000005; |
---|
133 | S0bd=0.000006; |
---|
134 | S0cd=0.000007; |
---|
135 | break; |
---|
136 | case 3: |
---|
137 | S0aa=0.000001; |
---|
138 | S0ab=0.000002; |
---|
139 | S0ac=0.000003; |
---|
140 | S0ad=0.000004; |
---|
141 | S0bc=0.000005; |
---|
142 | S0bd=0.000006; |
---|
143 | break; |
---|
144 | case 4: |
---|
145 | S0aa=0.000001; |
---|
146 | S0ab=0.000002; |
---|
147 | S0ac=0.000003; |
---|
148 | S0ad=0.000004; |
---|
149 | break; |
---|
150 | case 5: |
---|
151 | S0ab=0.000001; |
---|
152 | S0ac=0.000002; |
---|
153 | S0ad=0.000003; |
---|
154 | S0bc=0.000004; |
---|
155 | S0bd=0.000005; |
---|
156 | S0cd=0.000006; |
---|
157 | break; |
---|
158 | case 6: |
---|
159 | S0ab=0.000001; |
---|
160 | S0ac=0.000002; |
---|
161 | S0ad=0.000003; |
---|
162 | S0bc=0.000004; |
---|
163 | S0bd=0.000005; |
---|
164 | break; |
---|
165 | case 7: |
---|
166 | S0ab=0.000001; |
---|
167 | S0ac=0.000002; |
---|
168 | S0ad=0.000003; |
---|
169 | break; |
---|
170 | case 8: |
---|
171 | S0ac=0.000001; |
---|
172 | S0ad=0.000002; |
---|
173 | S0bc=0.000003; |
---|
174 | S0bd=0.000004; |
---|
175 | break; |
---|
176 | default : //case 9: |
---|
177 | break; |
---|
178 | } |
---|
179 | #endif |
---|
180 | |
---|
181 | // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F |
---|
182 | const double Kaa = 0.0; |
---|
183 | const double Kbb = 0.0; |
---|
184 | const double Kcc = 0.0; |
---|
185 | //const double Kdd = 0.0; |
---|
186 | const double Zaa = Kaa - Kad - Kad; |
---|
187 | const double Zab = Kab - Kad - Kbd; |
---|
188 | const double Zac = Kac - Kad - Kcd; |
---|
189 | const double Zbb = Kbb - Kbd - Kbd; |
---|
190 | const double Zbc = Kbc - Kbd - Kcd; |
---|
191 | const double Zcc = Kcc - Kcd - Kcd; |
---|
192 | //printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac); |
---|
193 | //printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc); |
---|
194 | //printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc); |
---|
195 | |
---|
196 | // T = inv(S0) |
---|
197 | const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc |
---|
198 | - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc); |
---|
199 | const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT; |
---|
200 | const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT; |
---|
201 | const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT; |
---|
202 | const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT; |
---|
203 | const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT; |
---|
204 | const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT; |
---|
205 | |
---|
206 | //printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13); |
---|
207 | //printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23); |
---|
208 | //printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33); |
---|
209 | |
---|
210 | // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0) |
---|
211 | const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd) |
---|
212 | + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd) |
---|
213 | + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd); |
---|
214 | |
---|
215 | const double m=1.0/(S0dd-ZZ); |
---|
216 | |
---|
217 | // eq 18d: Y = inv(S0)s0 + e |
---|
218 | const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0; |
---|
219 | const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0; |
---|
220 | const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0; |
---|
221 | |
---|
222 | // N = mYY^T + \kappa^F |
---|
223 | const double N11 = m*Y1*Y1 + Zaa; |
---|
224 | const double N12 = m*Y1*Y2 + Zab; |
---|
225 | const double N13 = m*Y1*Y3 + Zac; |
---|
226 | const double N22 = m*Y2*Y2 + Zbb; |
---|
227 | const double N23 = m*Y2*Y3 + Zbc; |
---|
228 | const double N33 = m*Y3*Y3 + Zcc; |
---|
229 | |
---|
230 | //printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13); |
---|
231 | //printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23); |
---|
232 | //printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33); |
---|
233 | //printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac); |
---|
234 | //printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc); |
---|
235 | //printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc); |
---|
236 | |
---|
237 | // M = I + S0 N |
---|
238 | const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0; |
---|
239 | const double Mab = N11*S0ab + N12*S0bb + N13*S0bc; |
---|
240 | const double Mac = N11*S0ac + N12*S0bc + N13*S0cc; |
---|
241 | const double Mba = N12*S0aa + N22*S0ab + N23*S0ac; |
---|
242 | const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0; |
---|
243 | const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc; |
---|
244 | const double Mca = N13*S0aa + N23*S0ab + N33*S0ac; |
---|
245 | const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc; |
---|
246 | const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0; |
---|
247 | //printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac); |
---|
248 | //printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc); |
---|
249 | //printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc); |
---|
250 | |
---|
251 | // Q = inv(M) = inv(I + S0 N) |
---|
252 | const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc |
---|
253 | + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca); |
---|
254 | |
---|
255 | const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ; |
---|
256 | const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ; |
---|
257 | const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ; |
---|
258 | //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ; |
---|
259 | const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ; |
---|
260 | const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ; |
---|
261 | //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ; |
---|
262 | //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ; |
---|
263 | const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ; |
---|
264 | |
---|
265 | //printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13); |
---|
266 | //printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23); |
---|
267 | //printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33); |
---|
268 | // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F |
---|
269 | // eq A1 in the appendix |
---|
270 | // To solve for S, use: |
---|
271 | // S = inv(inv(S^0) + N) inv(S^0) S^0 |
---|
272 | // = inv(S^0 inv(S^0) + N) S^0 |
---|
273 | // = inv(I + S^0 N) S^0 |
---|
274 | // = Q S^0 |
---|
275 | const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac; |
---|
276 | const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac; |
---|
277 | const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac; |
---|
278 | const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc; |
---|
279 | const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc; |
---|
280 | const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc; |
---|
281 | // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below |
---|
282 | //const double S14=-S11-S12-S13; |
---|
283 | //const double S24=-S12-S22-S23; |
---|
284 | //const double S34=-S13-S23-S33; |
---|
285 | //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23); |
---|
286 | |
---|
287 | // eq 12 of Akcasu, 1990: I(q) = L^T S L |
---|
288 | // Note: eliminate cases without A and B polymers by setting Lij to 0 |
---|
289 | // Note: 1e-13 to convert from fm to cm for scattering length |
---|
290 | const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; |
---|
291 | const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; |
---|
292 | const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; |
---|
293 | const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; |
---|
294 | |
---|
295 | const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 |
---|
296 | + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13); |
---|
297 | |
---|
298 | return result; |
---|
299 | |
---|
300 | } |
---|