source: sasmodels/sasmodels/models/rpa.c @ 20c856a

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 20c856a was 20c856a, checked in by butler, 7 years ago

Final rpa model corrections (correct scale and units) and verified with
one of original authors on paper. Fixes #593.

while documentation was updated minimally to be usable it could use more
thorough work as can many which is the subject of other tickets. Also
annoying is that the example plot is useless. It is however
autogenerated from "default parameters" which is not helpful in this
particular case. Will submit new ticket for that.

  • Property mode set to 100644
File size: 18.1 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[],    // DEGREE OF POLYMERIZATION
9    double Phi[],  // VOL FRACTION
10    double v[],    // SPECIFIC VOLUME
11    double L[],    // SCATT. LENGTH
12    double b[],    // SEGMENT LENGTH
13    double Kab, double Kac, double Kad,  // CHI PARAM
14    double Kbc, double Kbd, double Kcd
15    )
16{
17  int icase = (int)case_num;
18
19  double Nab,Nac,Nad,Nbc,Nbd,Ncd;
20  double Phiab,Phiac,Phiad,Phibc,Phibd,Phicd;
21  double vab,vac,vad,vbc,vbd,vcd;
22  double m;
23  double Xa,Xb,Xc,Xd;
24  double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad;
25  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd;
26  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd;
27  double S0da,S0db,S0dc;
28  double Pdd,S0dd;
29  double Kaa,Kbb,Kcc;
30  double Kba,Kca,Kcb;
31  double Kda,Kdb,Kdc,Kdd;
32  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc;
33  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33;
34  double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33;
35  double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33;
36  double N11,N12,N13,N21,N22,N23,N31,N32,N33;
37  double M11,M12,M13,M21,M22,M23,M31,M32,M33;
38  double S11,S12,S13,S14,S21,S22,S23,S24;
39  double S31,S32,S33,S34,S41,S42,S43,S44;
40  double Lad,Lbd,Lcd,Nav,Intg;
41 
42  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc)
43  //icase was shifted to N-1 from the original code
44  if (icase <= 1){
45    Phi[0] = Phi[1] = 0.0000001;
46    N[0] = N[1] = 1000.0;
47    L[0] = L[1] = 1.e-12;
48    v[0] = v[1] = 100.0;
49    b[0] = b[1] = 5.0;
50    Kab = Kac = Kad = Kbc = Kbd = -0.0004;
51  }
52  else if ((icase > 1) && (icase <= 4)){
53    Phi[0] = 0.0000001;
54    N[0] = 1000.0;
55    L[0] = 1.e-12;
56    v[0] = 100.0;
57    b[0] = 5.0;
58    Kab = Kac = Kad = -0.0004;
59  }
60 
61  // Set volume fraction of component D based on constraint that sum of vol frac =1
62  Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2];
63
64  //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9)
65  Nab=sqrt(N[0]*N[1]);
66  Nac=sqrt(N[0]*N[2]);
67  Nad=sqrt(N[0]*N[3]);
68  Nbc=sqrt(N[1]*N[2]);
69  Nbd=sqrt(N[1]*N[3]);
70  Ncd=sqrt(N[2]*N[3]);
71
72  vab=sqrt(v[0]*v[1]);
73  vac=sqrt(v[0]*v[2]);
74  vad=sqrt(v[0]*v[3]);
75  vbc=sqrt(v[1]*v[2]);
76  vbd=sqrt(v[1]*v[3]);
77  vcd=sqrt(v[2]*v[3]);
78
79  Phiab=sqrt(Phi[0]*Phi[1]);
80  Phiac=sqrt(Phi[0]*Phi[2]);
81  Phiad=sqrt(Phi[0]*Phi[3]);
82  Phibc=sqrt(Phi[1]*Phi[2]);
83  Phibd=sqrt(Phi[1]*Phi[3]);
84  Phicd=sqrt(Phi[2]*Phi[3]);
85
86  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk
87  Xa=q*q*b[0]*b[0]*N[0]/6.0;
88  Xb=q*q*b[1]*b[1]*N[1]/6.0;
89  Xc=q*q*b[2]*b[2]*N[2]/6.0;
90  Xd=q*q*b[3]*b[3]*N[3]/6.0;
91
92  //calculate all partial structure factors Pij and normalize n^2
93  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor
94  S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2
95  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor
96  S0ab=(Phiab*vab*Nab)*Pab;
97  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor
98  S0ac=(Phiac*vac*Nac)*Pac;
99  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block
100  S0ad=(Phiad*vad*Nad)*Pad;
101
102  S0ba=S0ab;
103  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain
104  S0bb=N[1]*Phi[1]*v[1]*Pbb;
105  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock
106  S0bc=(Phibc*vbc*Nbc)*Pbc;
107  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock
108  S0bd=(Phibd*vbd*Nbd)*Pbd;
109
110  S0ca=S0ac;
111  S0cb=S0bc;
112  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain
113  S0cc=N[2]*Phi[2]*v[2]*Pcc;
114  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock
115  S0cd=(Phicd*vcd*Ncd)*Pcd;
116
117  S0da=S0ad;
118  S0db=S0bd;
119  S0dc=S0cd;
120  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain
121  S0dd=N[3]*Phi[3]*v[3]*Pdd;
122
123  // Reset all unused partial structure factors to 0 (depends on case)
124  //icase was shifted to N-1 from the original code
125  switch(icase){
126  case 0:
127    S0aa=0.000001;
128    S0ab=0.000002;
129    S0ac=0.000003;
130    S0ad=0.000004;
131    S0bb=0.000005;
132    S0bc=0.000006;
133    S0bd=0.000007;
134    S0cd=0.000008;
135    break;
136  case 1:
137    S0aa=0.000001;
138    S0ab=0.000002;
139    S0ac=0.000003;
140    S0ad=0.000004;
141    S0bb=0.000005;
142    S0bc=0.000006;
143    S0bd=0.000007;
144    break;
145  case 2:
146    S0aa=0.000001;
147    S0ab=0.000002;
148    S0ac=0.000003;
149    S0ad=0.000004;
150    S0bc=0.000005;
151    S0bd=0.000006;
152    S0cd=0.000007;
153    break;
154  case 3:
155    S0aa=0.000001;
156    S0ab=0.000002;
157    S0ac=0.000003;
158    S0ad=0.000004;
159    S0bc=0.000005;
160    S0bd=0.000006;
161    break;
162  case 4:
163    S0aa=0.000001;
164    S0ab=0.000002;
165    S0ac=0.000003;
166    S0ad=0.000004;
167    break;
168  case 5:
169    S0ab=0.000001;
170    S0ac=0.000002;
171    S0ad=0.000003;
172    S0bc=0.000004;
173    S0bd=0.000005;
174    S0cd=0.000006;
175    break;
176  case 6:
177    S0ab=0.000001;
178    S0ac=0.000002;
179    S0ad=0.000003;
180    S0bc=0.000004;
181    S0bd=0.000005;
182    break;
183  case 7:
184    S0ab=0.000001;
185    S0ac=0.000002;
186    S0ad=0.000003;
187    break;
188  case 8:
189    S0ac=0.000001;
190    S0ad=0.000002;
191    S0bc=0.000003;
192    S0bd=0.000004;
193    break;
194  default : //case 9:
195    break;
196  }
197  S0ba=S0ab;
198  S0ca=S0ac;
199  S0cb=S0bc;
200  S0da=S0ad;
201  S0db=S0bd;
202  S0dc=S0cd;
203
204  // self chi parameter is 0 ... of course
205  Kaa=0.0;
206  Kbb=0.0;
207  Kcc=0.0;
208  Kdd=0.0;
209
210  Kba=Kab;
211  Kca=Kac;
212  Kcb=Kbc;
213  Kda=Kad;
214  Kdb=Kbd;
215  Kdc=Kcd;
216
217  Zaa=Kaa-Kad-Kad;
218  Zab=Kab-Kad-Kbd;
219  Zac=Kac-Kad-Kcd;
220  Zba=Kba-Kbd-Kad;
221  Zbb=Kbb-Kbd-Kbd;
222  Zbc=Kbc-Kbd-Kcd;
223  Zca=Kca-Kcd-Kad;
224  Zcb=Kcb-Kcd-Kbd;
225  Zcc=Kcc-Kcd-Kcd;
226
227  DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc);
228
229  T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT;
230  T12= (S0ac*S0cb - S0ab*S0cc)/DenT;
231  T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT;
232  T21= (S0bc*S0ca - S0ba*S0cc)/DenT;
233  T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT;
234  T23= (S0ac*S0ba - S0aa*S0bc)/DenT;
235  T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT;
236  T32= (S0ab*S0ca - S0aa*S0cb)/DenT;
237  T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT;
238
239  Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0;
240  Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0;
241  Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0;
242
243  X11=Y1*Y1;
244  X12=Y1*Y2;
245  X13=Y1*Y3;
246  X21=Y2*Y1;
247  X22=Y2*Y2;
248  X23=Y2*Y3;
249  X31=Y3*Y1;
250  X32=Y3*Y2;
251  X33=Y3*Y3;
252
253  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd);
254
255  // D is considered the matrix or background component so enters here
256  m=1.0/(S0dd-ZZ);
257
258  N11=m*X11+Zaa;
259  N12=m*X12+Zab;
260  N13=m*X13+Zac;
261  N21=m*X21+Zba;
262  N22=m*X22+Zbb;
263  N23=m*X23+Zbc;
264  N31=m*X31+Zca;
265  N32=m*X32+Zcb;
266  N33=m*X33+Zcc;
267
268  M11= N11*S0aa + N12*S0ab + N13*S0ac;
269  M12= N11*S0ab + N12*S0bb + N13*S0bc;
270  M13= N11*S0ac + N12*S0bc + N13*S0cc;
271  M21= N21*S0aa + N22*S0ab + N23*S0ac;
272  M22= N21*S0ab + N22*S0bb + N23*S0bc;
273  M23= N21*S0ac + N22*S0bc + N23*S0cc;
274  M31= N31*S0aa + N32*S0ab + N33*S0ac;
275  M32= N31*S0ab + N32*S0bb + N33*S0bc;
276  M33= N31*S0ac + N32*S0bc + N33*S0cc;
277
278  DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31;
279  DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33;
280  DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33;
281  DenQ=DenQ1+DenQ2+DenQ3;
282
283  Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ;
284  Q12= (-M12 + M13*M32 - M12*M33)/DenQ;
285  Q13= (-M13 - M13*M22 + M12*M23)/DenQ;
286  Q21= (-M21 + M23*M31 - M21*M33)/DenQ;
287  Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ;
288  Q23= (M13*M21 - M23 - M11*M23)/DenQ;
289  Q31= (-M31 - M22*M31 + M21*M32)/DenQ;
290  Q32= (M12*M31 - M32 - M11*M32)/DenQ;
291  Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ;
292
293  S11= Q11*S0aa + Q21*S0ab + Q31*S0ac;
294  S12= Q12*S0aa + Q22*S0ab + Q32*S0ac;
295  S13= Q13*S0aa + Q23*S0ab + Q33*S0ac;
296  S14=-S11-S12-S13;
297  S21= Q11*S0ba + Q21*S0bb + Q31*S0bc;
298  S22= Q12*S0ba + Q22*S0bb + Q32*S0bc;
299  S23= Q13*S0ba + Q23*S0bb + Q33*S0bc;
300  S24=-S21-S22-S23;
301  S31= Q11*S0ca + Q21*S0cb + Q31*S0cc;
302  S32= Q12*S0ca + Q22*S0cb + Q32*S0cc;
303  S33= Q13*S0ca + Q23*S0cb + Q33*S0cc;
304  S34=-S31-S32-S33;
305  S41=S14;
306  S42=S24;
307  S43=S34;
308  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23;
309
310  //calculate contrast where L[i] is the scattering length of i and D is the matrix
311  //Note that should multiply by Nav to get units of SLD which will become
312  // Nav*2 in the next line (SLD^2) but then normalization in that line would
313  //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.
314  Nav=6.022045e+23;
315  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav);
316  Lbd=(L[1]/v[1]-L[3]/v[3])*sqrt(Nav);
317  Lcd=(L[2]/v[2]-L[3]/v[3])*sqrt(Nav);
318
319  Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13;
320
321  //rescale for units of Lij^2 (fm^2 to cm^2)
322  Intg *= 1.0e-26;   
323
324  return Intg;
325
326
327/*  Attempts at a new implementation --- supressed for now
328#if 1  // Sasview defaults
329  if (icase <= 1) {
330    N[0]=N[1]=1000.0;
331    Phi[0]=Phi[1]=0.0000001;
332    Kab=Kac=Kad=Kbc=Kbd=-0.0004;
333    L[0]=L[1]=1.0e-12;
334    v[0]=v[1]=100.0;
335    b[0]=b[1]=5.0;
336  } else if (icase <= 4) {
337    Phi[0]=0.0000001;
338    Kab=Kac=Kad=-0.0004;
339    L[0]=1.0e-12;
340    v[0]=100.0;
341    b[0]=5.0;
342  }
343#else
344  if (icase <= 1) {
345    N[0]=N[1]=0.0;
346    Phi[0]=Phi[1]=0.0;
347    Kab=Kac=Kad=Kbc=Kbd=0.0;
348    L[0]=L[1]=L[3];
349    v[0]=v[1]=v[3];
350    b[0]=b[1]=0.0;
351  } else if (icase <= 4) {
352    N[0] = 0.0;
353    Phi[0]=0.0;
354    Kab=Kac=Kad=0.0;
355    L[0]=L[3];
356    v[0]=v[3];
357    b[0]=0.0;
358  }
359#endif
360
361  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
362  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
363  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
364  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;
365
366  // limit as Xa goes to 0 is 1
367  const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa;
368  const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb;
369  const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc;
370  const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd;
371
372  // limit as Xa goes to 0 is 1
373  const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa;
374  const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb;
375  const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc;
376  const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd;
377
378
379  // Note: S0ij only defined for copolymers; otherwise set to zero
380  // 0: C/D     binary mixture
381  // 1: C-D     diblock copolymer
382  // 2: B/C/D   ternery mixture
383  // 3: B/C-D   binary mixture,1 homopolymer, 1 diblock copolymer
384  // 4: B-C-D   triblock copolymer
385  // 5: A/B/C/D quaternary mixture
386  // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer
387  // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer
388  // 8: A-B/C-D binary mixture, 2 diblock copolymer
389  // 9: A-B-C-D tetra-block copolymer
390#if 0
391  const double S0aa = icase<5
392                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
393  const double S0bb = icase<2
394                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
395  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
396  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
397  const double S0ab = icase<8
398                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
399  const double S0ac = icase<9
400                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
401  const double S0ad = icase<9
402                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
403  const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
404                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
405  const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
406                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
407  const double S0cd = (icase==0 || icase==2 || icase==5)
408                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
409#else  // sasview equivalent
410//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
411  double S0aa = N[0]*Phi[0]*v[0]*Paa;
412  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
413  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
414  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
415  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
416  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
417  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
418  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
419  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
420  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
421switch(icase){
422  case 0:
423    S0aa=0.000001;
424    S0ab=0.000002;
425    S0ac=0.000003;
426    S0ad=0.000004;
427    S0bb=0.000005;
428    S0bc=0.000006;
429    S0bd=0.000007;
430    S0cd=0.000008;
431    break;
432  case 1:
433    S0aa=0.000001;
434    S0ab=0.000002;
435    S0ac=0.000003;
436    S0ad=0.000004;
437    S0bb=0.000005;
438    S0bc=0.000006;
439    S0bd=0.000007;
440    break;
441  case 2:
442    S0aa=0.000001;
443    S0ab=0.000002;
444    S0ac=0.000003;
445    S0ad=0.000004;
446    S0bc=0.000005;
447    S0bd=0.000006;
448    S0cd=0.000007;
449    break;
450  case 3:
451    S0aa=0.000001;
452    S0ab=0.000002;
453    S0ac=0.000003;
454    S0ad=0.000004;
455    S0bc=0.000005;
456    S0bd=0.000006;
457    break;
458  case 4:
459    S0aa=0.000001;
460    S0ab=0.000002;
461    S0ac=0.000003;
462    S0ad=0.000004;
463    break;
464  case 5:
465    S0ab=0.000001;
466    S0ac=0.000002;
467    S0ad=0.000003;
468    S0bc=0.000004;
469    S0bd=0.000005;
470    S0cd=0.000006;
471    break;
472  case 6:
473    S0ab=0.000001;
474    S0ac=0.000002;
475    S0ad=0.000003;
476    S0bc=0.000004;
477    S0bd=0.000005;
478    break;
479  case 7:
480    S0ab=0.000001;
481    S0ac=0.000002;
482    S0ad=0.000003;
483    break;
484  case 8:
485    S0ac=0.000001;
486    S0ad=0.000002;
487    S0bc=0.000003;
488    S0bd=0.000004;
489    break;
490  default : //case 9:
491    break;
492  }
493#endif
494
495  // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F
496  const double Kaa = 0.0;
497  const double Kbb = 0.0;
498  const double Kcc = 0.0;
499  //const double Kdd = 0.0;
500  const double Zaa = Kaa - Kad - Kad;
501  const double Zab = Kab - Kad - Kbd;
502  const double Zac = Kac - Kad - Kcd;
503  const double Zbb = Kbb - Kbd - Kbd;
504  const double Zbc = Kbc - Kbd - Kcd;
505  const double Zcc = Kcc - Kcd - Kcd;
506//printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac);
507//printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc);
508//printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc);
509
510  // T = inv(S0)
511  const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc
512                       - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc);
513  const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT;
514  const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT;
515  const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT;
516  const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT;
517  const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT;
518  const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT;
519
520//printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13);
521//printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23);
522//printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33);
523
524  // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0)
525  const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd)
526                  + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd)
527                  + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd);
528
529  const double m=1.0/(S0dd-ZZ);
530
531  // eq 18d: Y = inv(S0)s0 + e
532  const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0;
533  const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0;
534  const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0;
535
536  // N = mYY^T + \kappa^F
537  const double N11 = m*Y1*Y1 + Zaa;
538  const double N12 = m*Y1*Y2 + Zab;
539  const double N13 = m*Y1*Y3 + Zac;
540  const double N22 = m*Y2*Y2 + Zbb;
541  const double N23 = m*Y2*Y3 + Zbc;
542  const double N33 = m*Y3*Y3 + Zcc;
543
544//printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13);
545//printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23);
546//printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33);
547//printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac);
548//printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc);
549//printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc);
550
551  // M = I + S0 N
552  const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0;
553  const double Mab = N11*S0ab + N12*S0bb + N13*S0bc;
554  const double Mac = N11*S0ac + N12*S0bc + N13*S0cc;
555  const double Mba = N12*S0aa + N22*S0ab + N23*S0ac;
556  const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0;
557  const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc;
558  const double Mca = N13*S0aa + N23*S0ab + N33*S0ac;
559  const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc;
560  const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0;
561//printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac);
562//printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc);
563//printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc);
564
565  // Q = inv(M) = inv(I + S0 N)
566  const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc
567                       + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca);
568
569  const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ;
570  const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ;
571  const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ;
572  //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
573  const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ;
574  const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ;
575  //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
576  //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
577  const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ;
578
579//printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13);
580//printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23);
581//printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33);
582  // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F
583  // eq A1 in the appendix
584  // To solve for S, use:
585  //      S = inv(inv(S^0) + N) inv(S^0) S^0
586  //        = inv(S^0 inv(S^0) + N) S^0
587  //        = inv(I + S^0 N) S^0
588  //        = Q S^0
589  const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac;
590  const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac;
591  const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac;
592  const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc;
593  const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc;
594  const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc;
595  // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below
596  //const double S14=-S11-S12-S13;
597  //const double S24=-S12-S22-S23;
598  //const double S34=-S13-S23-S33;
599  //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23);
600
601  // eq 12 of Akcasu, 1990: I(q) = L^T S L
602  // Note: eliminate cases without A and B polymers by setting Lij to 0
603  // Note: 1e-13 to convert from fm to cm for scattering length
604  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
605  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
606  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
607  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;
608
609  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
610                    + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13);
611
612  return result;
613*/
614}
Note: See TracBrowser for help on using the repository browser.