source: sasmodels/sasmodels/models/rpa.c @ 639c4e3

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 639c4e3 was 639c4e3, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

Merge branch 'master' into polydisp

Conflicts:

sasmodels/kerneldll.py
sasmodels/models/rpa.c
  • Property mode set to 100644
File size: 9.9 KB
Line 
1double 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
7double 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;
107switch(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}
Note: See TracBrowser for help on using the repository browser.