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

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

scale intensity returned from rpa.c code to be close to correct
(currently using 1e-24). However th exact power of 10 will depend on
the units of Li (the scattering lengths) that are assumed in the
computation (we are now telling users it is in fm but clearly the code
is assuming something else — need to verify with rpa authors. Once done
this should complete work on #593. The documentation could use a lot of
TLC but that is a different ticket.
  • Property mode set to 100644
File size: 18.0 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  //need to verify why the sqrt of Nav rather than just Nav (assuming v is molar volume)
312  Nav=6.022045e+23;
313  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav);
314  Lbd=(L[1]/v[1]-L[3]/v[3])*sqrt(Nav);
315  Lcd=(L[2]/v[2]-L[3]/v[3])*sqrt(Nav);
316
317  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;
318
319  //rescale for units of Lij^2 (in 10e-12 m^2 to m^2 ?)
320  Intg *= 1.0e-24;   
321
322  return Intg;
323
324
325/*  Attempts at a new implementation --- supressed for now
326#if 1  // Sasview defaults
327  if (icase <= 1) {
328    N[0]=N[1]=1000.0;
329    Phi[0]=Phi[1]=0.0000001;
330    Kab=Kac=Kad=Kbc=Kbd=-0.0004;
331    L[0]=L[1]=1.0e-12;
332    v[0]=v[1]=100.0;
333    b[0]=b[1]=5.0;
334  } else if (icase <= 4) {
335    Phi[0]=0.0000001;
336    Kab=Kac=Kad=-0.0004;
337    L[0]=1.0e-12;
338    v[0]=100.0;
339    b[0]=5.0;
340  }
341#else
342  if (icase <= 1) {
343    N[0]=N[1]=0.0;
344    Phi[0]=Phi[1]=0.0;
345    Kab=Kac=Kad=Kbc=Kbd=0.0;
346    L[0]=L[1]=L[3];
347    v[0]=v[1]=v[3];
348    b[0]=b[1]=0.0;
349  } else if (icase <= 4) {
350    N[0] = 0.0;
351    Phi[0]=0.0;
352    Kab=Kac=Kad=0.0;
353    L[0]=L[3];
354    v[0]=v[3];
355    b[0]=0.0;
356  }
357#endif
358
359  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
360  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
361  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
362  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;
363
364  // limit as Xa goes to 0 is 1
365  const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa;
366  const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb;
367  const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc;
368  const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd;
369
370  // limit as Xa goes to 0 is 1
371  const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa;
372  const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb;
373  const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc;
374  const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd;
375
376
377  // Note: S0ij only defined for copolymers; otherwise set to zero
378  // 0: C/D     binary mixture
379  // 1: C-D     diblock copolymer
380  // 2: B/C/D   ternery mixture
381  // 3: B/C-D   binary mixture,1 homopolymer, 1 diblock copolymer
382  // 4: B-C-D   triblock copolymer
383  // 5: A/B/C/D quaternary mixture
384  // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer
385  // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer
386  // 8: A-B/C-D binary mixture, 2 diblock copolymer
387  // 9: A-B-C-D tetra-block copolymer
388#if 0
389  const double S0aa = icase<5
390                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
391  const double S0bb = icase<2
392                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
393  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
394  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
395  const double S0ab = icase<8
396                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
397  const double S0ac = icase<9
398                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
399  const double S0ad = icase<9
400                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
401  const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
402                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
403  const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
404                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
405  const double S0cd = (icase==0 || icase==2 || icase==5)
406                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
407#else  // sasview equivalent
408//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
409  double S0aa = N[0]*Phi[0]*v[0]*Paa;
410  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
411  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
412  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
413  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
414  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
415  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
416  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
417  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
418  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
419switch(icase){
420  case 0:
421    S0aa=0.000001;
422    S0ab=0.000002;
423    S0ac=0.000003;
424    S0ad=0.000004;
425    S0bb=0.000005;
426    S0bc=0.000006;
427    S0bd=0.000007;
428    S0cd=0.000008;
429    break;
430  case 1:
431    S0aa=0.000001;
432    S0ab=0.000002;
433    S0ac=0.000003;
434    S0ad=0.000004;
435    S0bb=0.000005;
436    S0bc=0.000006;
437    S0bd=0.000007;
438    break;
439  case 2:
440    S0aa=0.000001;
441    S0ab=0.000002;
442    S0ac=0.000003;
443    S0ad=0.000004;
444    S0bc=0.000005;
445    S0bd=0.000006;
446    S0cd=0.000007;
447    break;
448  case 3:
449    S0aa=0.000001;
450    S0ab=0.000002;
451    S0ac=0.000003;
452    S0ad=0.000004;
453    S0bc=0.000005;
454    S0bd=0.000006;
455    break;
456  case 4:
457    S0aa=0.000001;
458    S0ab=0.000002;
459    S0ac=0.000003;
460    S0ad=0.000004;
461    break;
462  case 5:
463    S0ab=0.000001;
464    S0ac=0.000002;
465    S0ad=0.000003;
466    S0bc=0.000004;
467    S0bd=0.000005;
468    S0cd=0.000006;
469    break;
470  case 6:
471    S0ab=0.000001;
472    S0ac=0.000002;
473    S0ad=0.000003;
474    S0bc=0.000004;
475    S0bd=0.000005;
476    break;
477  case 7:
478    S0ab=0.000001;
479    S0ac=0.000002;
480    S0ad=0.000003;
481    break;
482  case 8:
483    S0ac=0.000001;
484    S0ad=0.000002;
485    S0bc=0.000003;
486    S0bd=0.000004;
487    break;
488  default : //case 9:
489    break;
490  }
491#endif
492
493  // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F
494  const double Kaa = 0.0;
495  const double Kbb = 0.0;
496  const double Kcc = 0.0;
497  //const double Kdd = 0.0;
498  const double Zaa = Kaa - Kad - Kad;
499  const double Zab = Kab - Kad - Kbd;
500  const double Zac = Kac - Kad - Kcd;
501  const double Zbb = Kbb - Kbd - Kbd;
502  const double Zbc = Kbc - Kbd - Kcd;
503  const double Zcc = Kcc - Kcd - Kcd;
504//printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac);
505//printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc);
506//printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc);
507
508  // T = inv(S0)
509  const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc
510                       - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc);
511  const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT;
512  const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT;
513  const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT;
514  const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT;
515  const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT;
516  const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT;
517
518//printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13);
519//printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23);
520//printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33);
521
522  // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0)
523  const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd)
524                  + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd)
525                  + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd);
526
527  const double m=1.0/(S0dd-ZZ);
528
529  // eq 18d: Y = inv(S0)s0 + e
530  const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0;
531  const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0;
532  const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0;
533
534  // N = mYY^T + \kappa^F
535  const double N11 = m*Y1*Y1 + Zaa;
536  const double N12 = m*Y1*Y2 + Zab;
537  const double N13 = m*Y1*Y3 + Zac;
538  const double N22 = m*Y2*Y2 + Zbb;
539  const double N23 = m*Y2*Y3 + Zbc;
540  const double N33 = m*Y3*Y3 + Zcc;
541
542//printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13);
543//printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23);
544//printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33);
545//printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac);
546//printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc);
547//printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc);
548
549  // M = I + S0 N
550  const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0;
551  const double Mab = N11*S0ab + N12*S0bb + N13*S0bc;
552  const double Mac = N11*S0ac + N12*S0bc + N13*S0cc;
553  const double Mba = N12*S0aa + N22*S0ab + N23*S0ac;
554  const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0;
555  const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc;
556  const double Mca = N13*S0aa + N23*S0ab + N33*S0ac;
557  const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc;
558  const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0;
559//printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac);
560//printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc);
561//printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc);
562
563  // Q = inv(M) = inv(I + S0 N)
564  const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc
565                       + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca);
566
567  const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ;
568  const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ;
569  const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ;
570  //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
571  const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ;
572  const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ;
573  //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
574  //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
575  const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ;
576
577//printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13);
578//printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23);
579//printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33);
580  // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F
581  // eq A1 in the appendix
582  // To solve for S, use:
583  //      S = inv(inv(S^0) + N) inv(S^0) S^0
584  //        = inv(S^0 inv(S^0) + N) S^0
585  //        = inv(I + S^0 N) S^0
586  //        = Q S^0
587  const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac;
588  const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac;
589  const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac;
590  const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc;
591  const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc;
592  const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc;
593  // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below
594  //const double S14=-S11-S12-S13;
595  //const double S24=-S12-S22-S23;
596  //const double S34=-S13-S23-S33;
597  //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23);
598
599  // eq 12 of Akcasu, 1990: I(q) = L^T S L
600  // Note: eliminate cases without A and B polymers by setting Lij to 0
601  // Note: 1e-13 to convert from fm to cm for scattering length
602  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
603  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
604  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
605  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;
606
607  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
608                    + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13);
609
610  return result;
611*/
612}
Note: See TracBrowser for help on using the repository browser.