source: sasmodels/sasmodels/models/rpa.c @ 830cf6b

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since 830cf6b was 19dc29e7, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

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