source: sasview/sansmodels/src/c_models/rpa.cpp @ bbbed8c

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since bbbed8c was 7ffa8196, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

refactor RPA model
  • Property mode set to 100644
File size: 12.0 KB
Line 
1/**
2        This software was developed by the University of Tennessee as part of the
3        Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4        project funded by the US National Science Foundation.
5
6        If you use DANSE applications to do scientific research that leads to
7        publication, we ask that you acknowledge the use of the software with the
8        following sentence:
9
10        "This work benefited from DANSE software developed under NSF award DMR-0520547."
11
12        copyright 2008, University of Tennessee
13 */
14
15/**
16 * Scattering model classes
17 * The classes use the IGOR library found in
18 *   sansmodels/src/libigor
19 *
20 */
21
22#include <math.h>
23#include "parameters.hh"
24#include <stdio.h>
25#include "rpa.h"
26
27using namespace std;
28
29static double rpa_kernel(double dp[], double q) {
30  int lCASE = dp[0];
31
32  double Na,Nb,Nc,Nd,Nab,Nac,Nad,Nbc,Nbd,Ncd;
33  double Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad;
34  double Phibc,Phibd,Phicd;
35  double va,vb,vc,vd,vab,vac,vad,vbc,vbd,vcd;
36  double m;
37  double ba,bb,bc,bd;
38  double Xa,Xb,Xc,Xd;
39  double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad;
40  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd;
41  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd;
42  double S0da,S0db,S0dc,Pdd,S0dd;
43  double Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd;
44  double Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd;
45  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc;
46  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33;
47  double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33;
48  double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33;
49  double N11,N12,N13,N21,N22,N23,N31,N32,N33;
50  double M11,M12,M13,M21,M22,M23,M31,M32,M33;
51  double S11,S12,S13,S14,S21,S22,S23,S24;
52  double S31,S32,S33,S34,S41,S42,S43,S44;
53  double La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Intg;
54  double scale, background;
55
56  int ii=13;  //dp[ii<13] = fittable, else fixed
57
58  Na=1000.0;
59  Nb=1000.0;
60  Nc=1000.0;
61  Nd=1000.0;  //DEGREE OF POLYMERIZATION
62  Phia=0.25;
63  Phib=0.25;
64  Phic=0.25;
65  Phid=0.25 ; //VOL FRACTION
66  Kab=-0.0004;
67  Kac=-0.0004;
68  Kad=-0.0004;  //CHI PARAM
69  Kbc=-0.0004;
70  Kbd=-0.0004;
71  Kcd=-0.0004;
72  La=1.0e-12;
73  Lb=1.0e-12;
74  Lc=1.0e-12;
75  Ld=0.0e-12;   //SCATT. LENGTH
76  va=100.0;
77  vb=100.0;
78  vc=100.0;
79  vd=100.0  ; //SPECIFIC VOLUME
80  ba=5.0;
81  bb=5.0;
82  bc=5.0;
83  bd=5.0;   //SEGMENT LENGTH
84
85  //lCASE was shifted to N-1 from the original code
86  if (lCASE <= 1){
87    Phia=0.0000001;
88    Phib=0.0000001;
89    Phic=0.5;
90    Phid=0.5;
91    Nc=dp[ii+8];
92    Phic=dp[ii+9];
93    vc=dp[ii+10];
94    Lc=dp[ii+11];
95    bc=dp[3];
96    Nd=dp[ii+12];
97    Phid=dp[ii+13];
98    vd=dp[ii+14];
99    Ld=dp[ii+15];
100    bd=dp[4];
101    Kcd=dp[10];
102    scale=dp[11];
103    background=dp[12];
104  }
105  else if ((lCASE > 1) && (lCASE <= 4)){
106    Phia=0.0000001;
107    Phib=0.333333;
108    Phic=0.333333;
109    Phid=0.333333;
110    Nb=dp[ii+4];
111    Phib=dp[ii+5];
112    vb=dp[ii+6];
113    Lb=dp[ii+7];
114    bb=dp[2];
115    Nc=dp[ii+8];
116    Phic=dp[ii+9];
117    vc=dp[ii+10];
118    Lc=dp[ii+11];
119    bc=dp[3];
120    Nd=dp[ii+12];
121    Phid=dp[ii+13];
122    vd=dp[ii+14];
123    Ld=dp[ii+15];
124    bd=dp[4];
125    Kbc=dp[8];
126    Kbd=dp[9];
127    Kcd=dp[10];
128    scale=dp[11];
129    background=dp[12];
130  }
131  else {
132    Phia=0.25;
133    Phib=0.25;
134    Phic=0.25;
135    Phid=0.25;
136    Na=dp[ii+0];
137    Phia=dp[ii+1];
138    va=dp[ii+2];
139    La=dp[ii+3];
140    ba=dp[1];
141    Nb=dp[ii+4];
142    Phib=dp[ii+5];
143    vb=dp[ii+6];
144    Lb=dp[ii+7];
145    bb=dp[2];
146    Nc=dp[ii+8];
147    Phic=dp[ii+9];
148    vc=dp[ii+10];
149    Lc=dp[ii+11];
150    bc=dp[3];
151    Nd=dp[ii+12];
152    Phid=dp[ii+13];
153    vd=dp[ii+14];
154    Ld=dp[ii+15];
155    bd=dp[4];
156    Kab=dp[5];
157    Kac=dp[6];
158    Kad=dp[7];
159    Kbc=dp[8];
160    Kbd=dp[9];
161    Kcd=dp[10];
162    scale=dp[11];
163    background=dp[12];
164  }
165
166  Nab=pow((Na*Nb),(0.5));
167  Nac=pow((Na*Nc),(0.5));
168  Nad=pow((Na*Nd),(0.5));
169  Nbc=pow((Nb*Nc),(0.5));
170  Nbd=pow((Nb*Nd),(0.5));
171  Ncd=pow((Nc*Nd),(0.5));
172
173  vab=pow((va*vb),(0.5));
174  vac=pow((va*vc),(0.5));
175  vad=pow((va*vd),(0.5));
176  vbc=pow((vb*vc),(0.5));
177  vbd=pow((vb*vd),(0.5));
178  vcd=pow((vc*vd),(0.5));
179
180  Phiab=pow((Phia*Phib),(0.5));
181  Phiac=pow((Phia*Phic),(0.5));
182  Phiad=pow((Phia*Phid),(0.5));
183  Phibc=pow((Phib*Phic),(0.5));
184  Phibd=pow((Phib*Phid),(0.5));
185  Phicd=pow((Phic*Phid),(0.5));
186
187  Xa=q*q*ba*ba*Na/6.0;
188  Xb=q*q*bb*bb*Nb/6.0;
189  Xc=q*q*bc*bb*Nc/6.0;
190  Xd=q*q*bd*bb*Nd/6.0;
191
192  Paa=2.0*(exp(-Xa)-1.0+Xa)/pow(Xa,2.0);
193  S0aa=Na*Phia*va*Paa;
194  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb);
195  S0ab=(Phiab*vab*Nab)*Pab;
196  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc);
197  S0ac=(Phiac*vac*Nac)*Pac;
198  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd);
199  S0ad=(Phiad*vad*Nad)*Pad;
200
201  S0ba=S0ab;
202  Pbb=2.0*(exp(-Xb)-1.0+Xb)/pow(Xb,2.0);
203  S0bb=Nb*Phib*vb*Pbb;
204  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc);
205  S0bc=(Phibc*vbc*Nbc)*Pbc;
206  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd);
207  S0bd=(Phibd*vbd*Nbd)*Pbd;
208
209  S0ca=S0ac;
210  S0cb=S0bc;
211  Pcc=2.0*(exp(-Xc)-1.0+Xc)/pow(Xc,2.0);
212  S0cc=Nc*Phic*vc*Pcc;
213  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd);
214  S0cd=(Phicd*vcd*Ncd)*Pcd;
215
216  S0da=S0ad;
217  S0db=S0bd;
218  S0dc=S0cd;
219  Pdd=2.0*(exp(-Xd)-1.0+Xd)/pow(Xd,2.0);
220  S0dd=Nd*Phid*vd*Pdd;
221  //lCASE was shifted to N-1 from the original code
222  switch(lCASE){
223  case 0:
224    S0aa=0.000001;
225    S0ab=0.000002;
226    S0ac=0.000003;
227    S0ad=0.000004;
228    S0bb=0.000005;
229    S0bc=0.000006;
230    S0bd=0.000007;
231    S0cd=0.000008;
232    break;
233  case 1:
234    S0aa=0.000001;
235    S0ab=0.000002;
236    S0ac=0.000003;
237    S0ad=0.000004;
238    S0bb=0.000005;
239    S0bc=0.000006;
240    S0bd=0.000007;
241    break;
242  case 2:
243    S0aa=0.000001;
244    S0ab=0.000002;
245    S0ac=0.000003;
246    S0ad=0.000004;
247    S0bc=0.000005;
248    S0bd=0.000006;
249    S0cd=0.000007;
250    break;
251  case 3:
252    S0aa=0.000001;
253    S0ab=0.000002;
254    S0ac=0.000003;
255    S0ad=0.000004;
256    S0bc=0.000005;
257    S0bd=0.000006;
258    break;
259  case 4:
260    S0aa=0.000001;
261    S0ab=0.000002;
262    S0ac=0.000003;
263    S0ad=0.000004;
264    break;
265  case 5:
266    S0ab=0.000001;
267    S0ac=0.000002;
268    S0ad=0.000003;
269    S0bc=0.000004;
270    S0bd=0.000005;
271    S0cd=0.000006;
272    break;
273  case 6:
274    S0ab=0.000001;
275    S0ac=0.000002;
276    S0ad=0.000003;
277    S0bc=0.000004;
278    S0bd=0.000005;
279    break;
280  case 7:
281    S0ab=0.000001;
282    S0ac=0.000002;
283    S0ad=0.000003;
284    break;
285  case 8:
286    S0ac=0.000001;
287    S0ad=0.000002;
288    S0bc=0.000003;
289    S0bd=0.000004;
290    break;
291  default : //case 9:
292    break;
293  }
294  S0ba=S0ab;
295  S0ca=S0ac;
296  S0cb=S0bc;
297  S0da=S0ad;
298  S0db=S0bd;
299  S0dc=S0cd;
300
301  Kaa=0.0;
302  Kbb=0.0;
303  Kcc=0.0;
304  Kdd=0.0;
305
306  Kba=Kab;
307  Kca=Kac;
308  Kcb=Kbc;
309  Kda=Kad;
310  Kdb=Kbd;
311  Kdc=Kcd;
312
313  Zaa=Kaa-Kad-Kad;
314  Zab=Kab-Kad-Kbd;
315  Zac=Kac-Kad-Kcd;
316  Zba=Kba-Kbd-Kad;
317  Zbb=Kbb-Kbd-Kbd;
318  Zbc=Kbc-Kbd-Kcd;
319  Zca=Kca-Kcd-Kad;
320  Zcb=Kcb-Kcd-Kbd;
321  Zcc=Kcc-Kcd-Kcd;
322
323  DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc);
324
325  T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT;
326  T12= (S0ac*S0cb - S0ab*S0cc)/DenT;
327  T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT;
328  T21= (S0bc*S0ca - S0ba*S0cc)/DenT;
329  T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT;
330  T23= (S0ac*S0ba - S0aa*S0bc)/DenT;
331  T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT;
332  T32= (S0ab*S0ca - S0aa*S0cb)/DenT;
333  T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT;
334
335  Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0;
336  Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0;
337  Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0;
338
339  X11=Y1*Y1;
340  X12=Y1*Y2;
341  X13=Y1*Y3;
342  X21=Y2*Y1;
343  X22=Y2*Y2;
344  X23=Y2*Y3;
345  X31=Y3*Y1;
346  X32=Y3*Y2;
347  X33=Y3*Y3;
348
349  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd);
350
351  m=1.0/(S0dd-ZZ);
352
353  N11=m*X11+Zaa;
354  N12=m*X12+Zab;
355  N13=m*X13+Zac;
356  N21=m*X21+Zba;
357  N22=m*X22+Zbb;
358  N23=m*X23+Zbc;
359  N31=m*X31+Zca;
360  N32=m*X32+Zcb;
361  N33=m*X33+Zcc;
362
363  M11= N11*S0aa + N12*S0ab + N13*S0ac;
364  M12= N11*S0ab + N12*S0bb + N13*S0bc;
365  M13= N11*S0ac + N12*S0bc + N13*S0cc;
366  M21= N21*S0aa + N22*S0ab + N23*S0ac;
367  M22= N21*S0ab + N22*S0bb + N23*S0bc;
368  M23= N21*S0ac + N22*S0bc + N23*S0cc;
369  M31= N31*S0aa + N32*S0ab + N33*S0ac;
370  M32= N31*S0ab + N32*S0bb + N33*S0bc;
371  M33= N31*S0ac + N32*S0bc + N33*S0cc;
372
373  DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31;
374  DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33;
375  DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33;
376  DenQ=DenQ1+DenQ2+DenQ3;
377
378  Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ;
379  Q12= (-M12 + M13*M32 - M12*M33)/DenQ;
380  Q13= (-M13 - M13*M22 + M12*M23)/DenQ;
381  Q21= (-M21 + M23*M31 - M21*M33)/DenQ;
382  Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ;
383  Q23= (M13*M21 - M23 - M11*M23)/DenQ;
384  Q31= (-M31 - M22*M31 + M21*M32)/DenQ;
385  Q32= (M12*M31 - M32 - M11*M32)/DenQ;
386  Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ;
387
388  S11= Q11*S0aa + Q21*S0ab + Q31*S0ac;
389  S12= Q12*S0aa + Q22*S0ab + Q32*S0ac;
390  S13= Q13*S0aa + Q23*S0ab + Q33*S0ac;
391  S14=-S11-S12-S13;
392  S21= Q11*S0ba + Q21*S0bb + Q31*S0bc;
393  S22= Q12*S0ba + Q22*S0bb + Q32*S0bc;
394  S23= Q13*S0ba + Q23*S0bb + Q33*S0bc;
395  S24=-S21-S22-S23;
396  S31= Q11*S0ca + Q21*S0cb + Q31*S0cc;
397  S32= Q12*S0ca + Q22*S0cb + Q32*S0cc;
398  S33= Q13*S0ca + Q23*S0cb + Q33*S0cc;
399  S34=-S31-S32-S33;
400  S41=S14;
401  S42=S24;
402  S43=S34;
403  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23;
404
405  Nav=6.022045e+23;
406  Lad=(La/va-Ld/vd)*sqrt(Nav);
407  Lbd=(Lb/vb-Ld/vd)*sqrt(Nav);
408  Lcd=(Lc/vc-Ld/vd)*sqrt(Nav);
409
410  Intg=pow(Lad,2.0)*S11+pow(Lbd,2.0)*S22+pow(Lcd,2.0)*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13;
411
412  Intg *= scale;
413  Intg += background;
414
415  return Intg;
416
417}
418
419RPAModel :: RPAModel() {
420  lcase_n = Parameter(0.0);
421  ba = Parameter(5.0);
422  bb = Parameter(5.0);
423  bc = Parameter(5.0);
424  bd = Parameter(5.0);
425
426  Kab = Parameter(-0.0004);
427  Kac = Parameter(-0.0004);
428  Kad = Parameter(-0.0004);
429  Kbc = Parameter(-0.0004);
430  Kbd = Parameter(-0.0004);
431  Kcd = Parameter(-0.0004);
432
433  scale = Parameter(1.0);
434  background = Parameter(0.0);
435
436  Na = Parameter(1000.0);
437  Phia = Parameter(0.25);
438  va = Parameter(100.0);
439  La = Parameter(1.0e-012);
440
441  Nb = Parameter(1000.0);
442  Phib = Parameter(0.25);
443  vb = Parameter(100.0);
444  Lb = Parameter(1.0e-012);
445
446  Nc = Parameter(1000.0);
447  Phic = Parameter(0.25);
448  vc = Parameter(100.0);
449  Lc = Parameter(1.0e-012);
450
451  Nd = Parameter(1000.0);
452  Phid = Parameter(0.25);
453  vd = Parameter(100.0);
454  Ld = Parameter(0.0e-012);
455
456}
457
458/**
459 * Function to evaluate 1D scattering function
460 * The NIST IGOR library is used for the actual calculation.
461 * @param q: q-value
462 * @return: function value
463 */
464double RPAModel :: operator()(double q) {
465  double dp[29];
466  // Fill parameter array for IGOR library
467  // Add the background after averaging
468  //Fittable parameters
469  dp[0] = lcase_n();
470
471  dp[1] = ba();
472  dp[2] = bb();
473  dp[3] = bc();
474  dp[4] = bd();
475
476  dp[5] = Kab();
477  dp[6] = Kac();
478  dp[7] = Kad();
479  dp[8] = Kbc();
480  dp[9] = Kbd();
481  dp[10] = Kcd();
482
483  dp[11] = scale();
484  dp[12] = background();
485
486  //Fixed parameters
487  dp[13] = Na();
488  dp[14] = Phia();
489  dp[15] = va();
490  dp[16] = La();
491
492  dp[17] = Nb();
493  dp[18] = Phib();
494  dp[19] = vb();
495  dp[20] = Lb();
496
497  dp[21] = Nc();
498  dp[22] = Phic();
499  dp[23] = vc();
500  dp[24] = Lc();
501
502  dp[25] = Nd();
503  dp[26] = Phid();
504  dp[27] = vd();
505  dp[28] = Ld();
506
507  return rpa_kernel(dp,q);
508}
509
510/**
511 * Function to evaluate 2D scattering function
512 * @param q_x: value of Q along x
513 * @param q_y: value of Q along y
514 * @return: function value
515 */
516double RPAModel :: operator()(double qx, double qy) {
517  double q = sqrt(qx*qx + qy*qy);
518  return (*this).operator()(q);
519}
520
521/**
522 * Function to evaluate 2D scattering function
523 * @param pars: parameters of the sphere
524 * @param q: q-value
525 * @param phi: angle phi
526 * @return: function value
527 */
528double RPAModel :: evaluate_rphi(double q, double phi) {
529  return (*this).operator()(q);
530}
531
532/**
533 * Function to calculate effective radius
534 * @return: effective radius value
535 */
536double RPAModel :: calculate_ER() {
537  //NOT implemented!!!
538  return 0.0;
539}
Note: See TracBrowser for help on using the repository browser.