source: sasmodels/sasmodels/models/rpa.c @ d507c3a

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

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