source: sasview/sansmodels/src/c_extensions/rpa.c @ 31af069

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 31af069 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Reorganizing models in preparation of cpp cleanup
  • Property mode set to 100644
File size: 10.0 KB
Line 
1/**
2 * Scattering model for multiblock copolymer
3 * The below rpa_kernel is modified to fit on 'C'
4 * originated from a IgorProcedure by B. HAMMOUDA, NIST, JULY 1998,
5 * and it was tried to keep it as original as possible: JHCho (10/05/10)
6 */
7#include <math.h>
8#include "rpa.h"
9#include <stdio.h>
10#include <stdlib.h>
11
12double rpa_kernel(double dp[], double q) {
13        int lCASE = dp[0];
14
15        double Na,Nb,Nc,Nd,Nab,Nac,Nad,Nbc,Nbd,Ncd;
16        double Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad;
17        double Phibc,Phibd,Phicd;
18        double va,vb,vc,vd,vab,vac,vad,vbc,vbd,vcd;
19        double m;
20        double ba,bb,bc,bd;
21        double Xa,Xb,Xc,Xd;
22        double Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad;
23        double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd;
24        double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd;
25        double S0da,S0db,S0dc,Pdd,S0dd;
26        double Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd;
27        double Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd;
28        double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc;
29        double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33;
30        double Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33;
31        double ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33;
32        double N11,N12,N13,N21,N22,N23,N31,N32,N33;
33        double M11,M12,M13,M21,M22,M23,M31,M32,M33;
34        double S11,S12,S13,S14,S21,S22,S23,S24;
35        double S31,S32,S33,S34,S41,S42,S43,S44;
36        double La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Intg;
37        double scale, background;
38
39        int ii=13;      //dp[ii<13] = fittable, else fixed
40
41        Na=1000.0;
42        Nb=1000.0;
43        Nc=1000.0;
44        Nd=1000.0;      //DEGREE OF POLYMERIZATION
45        Phia=0.25;
46        Phib=0.25;
47        Phic=0.25;
48        Phid=0.25 ;     //VOL FRACTION
49        Kab=-0.0004;
50        Kac=-0.0004;
51        Kad=-0.0004;    //CHI PARAM
52        Kbc=-0.0004;
53        Kbd=-0.0004;
54        Kcd=-0.0004;
55        La=1.0e-12;
56        Lb=1.0e-12;
57        Lc=1.0e-12;
58        Ld=0.0e-12;     //SCATT. LENGTH
59        va=100.0;
60        vb=100.0;
61        vc=100.0;
62        vd=100.0        ;       //SPECIFIC VOLUME
63        ba=5.0;
64        bb=5.0;
65        bc=5.0;
66        bd=5.0;         //SEGMENT LENGTH
67
68        //lCASE was shifted to N-1 from the original code
69        if (lCASE <= 1){
70                Phia=0.0000001;
71                Phib=0.0000001;
72                Phic=0.5;
73                Phid=0.5;
74                Nc=dp[ii+8];
75                Phic=dp[ii+9];
76                vc=dp[ii+10];
77                Lc=dp[ii+11];
78                bc=dp[3];
79                Nd=dp[ii+12];
80                Phid=dp[ii+13];
81                vd=dp[ii+14];
82                Ld=dp[ii+15];
83                bd=dp[4];
84                Kcd=dp[10];
85                scale=dp[11];
86                background=dp[12];
87        }
88        else if ((lCASE > 1) && (lCASE <= 4)){
89                Phia=0.0000001;
90                Phib=0.333333;
91                Phic=0.333333;
92                Phid=0.333333;
93                Nb=dp[ii+4];
94                Phib=dp[ii+5];
95                vb=dp[ii+6];
96                Lb=dp[ii+7];
97                bb=dp[2];
98                Nc=dp[ii+8];
99                Phic=dp[ii+9];
100                vc=dp[ii+10];
101                Lc=dp[ii+11];
102                bc=dp[3];
103                Nd=dp[ii+12];
104                Phid=dp[ii+13];
105                vd=dp[ii+14];
106                Ld=dp[ii+15];
107                bd=dp[4];
108                Kbc=dp[8];
109                Kbd=dp[9];
110                Kcd=dp[10];
111                scale=dp[11];
112                background=dp[12];
113        }
114        else {
115                Phia=0.25;
116                Phib=0.25;
117                Phic=0.25;
118                Phid=0.25;
119                Na=dp[ii+0];
120                Phia=dp[ii+1];
121                va=dp[ii+2];
122                La=dp[ii+3];
123                ba=dp[1];
124                Nb=dp[ii+4];
125                Phib=dp[ii+5];
126                vb=dp[ii+6];
127                Lb=dp[ii+7];
128                bb=dp[2];
129                Nc=dp[ii+8];
130                Phic=dp[ii+9];
131                vc=dp[ii+10];
132                Lc=dp[ii+11];
133                bc=dp[3];
134                Nd=dp[ii+12];
135                Phid=dp[ii+13];
136                vd=dp[ii+14];
137                Ld=dp[ii+15];
138                bd=dp[4];
139                Kab=dp[5];
140                Kac=dp[6];
141                Kad=dp[7];
142                Kbc=dp[8];
143                Kbd=dp[9];
144                Kcd=dp[10];
145                scale=dp[11];
146                background=dp[12];
147        }
148
149        Nab=pow((Na*Nb),(0.5));
150        Nac=pow((Na*Nc),(0.5));
151        Nad=pow((Na*Nd),(0.5));
152        Nbc=pow((Nb*Nc),(0.5));
153        Nbd=pow((Nb*Nd),(0.5));
154        Ncd=pow((Nc*Nd),(0.5));
155
156        vab=pow((va*vb),(0.5));
157        vac=pow((va*vc),(0.5));
158        vad=pow((va*vd),(0.5));
159        vbc=pow((vb*vc),(0.5));
160        vbd=pow((vb*vd),(0.5));
161        vcd=pow((vc*vd),(0.5));
162
163        Phiab=pow((Phia*Phib),(0.5));
164        Phiac=pow((Phia*Phic),(0.5));
165        Phiad=pow((Phia*Phid),(0.5));
166        Phibc=pow((Phib*Phic),(0.5));
167        Phibd=pow((Phib*Phid),(0.5));
168        Phicd=pow((Phic*Phid),(0.5));
169
170        Xa=q*q*ba*ba*Na/6.0;
171        Xb=q*q*bb*bb*Nb/6.0;
172        Xc=q*q*bc*bb*Nc/6.0;
173        Xd=q*q*bd*bb*Nd/6.0;
174
175        Paa=2.0*(exp(-Xa)-1.0+Xa)/pow(Xa,2.0);
176        S0aa=Na*Phia*va*Paa;
177        Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb);
178        S0ab=(Phiab*vab*Nab)*Pab;
179        Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc);
180        S0ac=(Phiac*vac*Nac)*Pac;
181        Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd);
182        S0ad=(Phiad*vad*Nad)*Pad;
183
184        S0ba=S0ab;
185        Pbb=2.0*(exp(-Xb)-1.0+Xb)/pow(Xb,2.0);
186        S0bb=Nb*Phib*vb*Pbb;
187        Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc);
188        S0bc=(Phibc*vbc*Nbc)*Pbc;
189        Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd);
190        S0bd=(Phibd*vbd*Nbd)*Pbd;
191
192        S0ca=S0ac;
193        S0cb=S0bc;
194        Pcc=2.0*(exp(-Xc)-1.0+Xc)/pow(Xc,2.0);
195        S0cc=Nc*Phic*vc*Pcc;
196        Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd);
197        S0cd=(Phicd*vcd*Ncd)*Pcd;
198
199        S0da=S0ad;
200        S0db=S0bd;
201        S0dc=S0cd;
202        Pdd=2.0*(exp(-Xd)-1.0+Xd)/pow(Xd,2.0);
203        S0dd=Nd*Phid*vd*Pdd;
204        //lCASE was shifted to N-1 from the original code
205        switch(lCASE){
206                case 0:
207                        S0aa=0.000001;
208                        S0ab=0.000002;
209                        S0ac=0.000003;
210                        S0ad=0.000004;
211                        S0bb=0.000005;
212                        S0bc=0.000006;
213                        S0bd=0.000007;
214                        S0cd=0.000008;
215                        break;
216                case 1:
217                        S0aa=0.000001;
218                        S0ab=0.000002;
219                        S0ac=0.000003;
220                        S0ad=0.000004;
221                        S0bb=0.000005;
222                        S0bc=0.000006;
223                        S0bd=0.000007;
224                        break;
225                case 2:
226                        S0aa=0.000001;
227                        S0ab=0.000002;
228                        S0ac=0.000003;
229                        S0ad=0.000004;
230                        S0bc=0.000005;
231                        S0bd=0.000006;
232                        S0cd=0.000007;
233                        break;
234                case 3:
235                        S0aa=0.000001;
236                        S0ab=0.000002;
237                        S0ac=0.000003;
238                        S0ad=0.000004;
239                        S0bc=0.000005;
240                        S0bd=0.000006;
241                        break;
242                case 4:
243                        S0aa=0.000001;
244                        S0ab=0.000002;
245                        S0ac=0.000003;
246                        S0ad=0.000004;
247                        break;
248                case 5:
249                        S0ab=0.000001;
250                        S0ac=0.000002;
251                        S0ad=0.000003;
252                        S0bc=0.000004;
253                        S0bd=0.000005;
254                        S0cd=0.000006;
255                        break;
256                case 6:
257                        S0ab=0.000001;
258                        S0ac=0.000002;
259                        S0ad=0.000003;
260                        S0bc=0.000004;
261                        S0bd=0.000005;
262                        break;
263                case 7:
264                        S0ab=0.000001;
265                        S0ac=0.000002;
266                        S0ad=0.000003;
267                        break;
268                case 8:
269                        S0ac=0.000001;
270                        S0ad=0.000002;
271                        S0bc=0.000003;
272                        S0bd=0.000004;
273                        break;
274                default : //case 9:
275                        break;
276        }
277        S0ba=S0ab;
278        S0ca=S0ac;
279        S0cb=S0bc;
280        S0da=S0ad;
281        S0db=S0bd;
282        S0dc=S0cd;
283
284        Kaa=0.0;
285        Kbb=0.0;
286        Kcc=0.0;
287        Kdd=0.0;
288
289        Kba=Kab;
290        Kca=Kac;
291        Kcb=Kbc;
292        Kda=Kad;
293        Kdb=Kbd;
294        Kdc=Kcd;
295
296        Zaa=Kaa-Kad-Kad;
297        Zab=Kab-Kad-Kbd;
298        Zac=Kac-Kad-Kcd;
299        Zba=Kba-Kbd-Kad;
300        Zbb=Kbb-Kbd-Kbd;
301        Zbc=Kbc-Kbd-Kcd;
302        Zca=Kca-Kcd-Kad;
303        Zcb=Kcb-Kcd-Kbd;
304        Zcc=Kcc-Kcd-Kcd;
305
306        DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc);
307
308        T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT;
309        T12= (S0ac*S0cb - S0ab*S0cc)/DenT;
310        T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT;
311        T21= (S0bc*S0ca - S0ba*S0cc)/DenT;
312        T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT;
313        T23= (S0ac*S0ba - S0aa*S0bc)/DenT;
314        T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT;
315        T32= (S0ab*S0ca - S0aa*S0cb)/DenT;
316        T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT;
317
318        Y1=T11*S0ad+T12*S0bd+T13*S0cd+1.0;
319        Y2=T21*S0ad+T22*S0bd+T23*S0cd+1.0;
320        Y3=T31*S0ad+T32*S0bd+T33*S0cd+1.0;
321
322        X11=Y1*Y1;
323        X12=Y1*Y2;
324        X13=Y1*Y3;
325        X21=Y2*Y1;
326        X22=Y2*Y2;
327        X23=Y2*Y3;
328        X31=Y3*Y1;
329        X32=Y3*Y2;
330        X33=Y3*Y3;
331
332        ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd);
333
334        m=1.0/(S0dd-ZZ);
335
336        N11=m*X11+Zaa;
337        N12=m*X12+Zab;
338        N13=m*X13+Zac;
339        N21=m*X21+Zba;
340        N22=m*X22+Zbb;
341        N23=m*X23+Zbc;
342        N31=m*X31+Zca;
343        N32=m*X32+Zcb;
344        N33=m*X33+Zcc;
345
346        M11= N11*S0aa + N12*S0ab + N13*S0ac;
347        M12= N11*S0ab + N12*S0bb + N13*S0bc;
348        M13= N11*S0ac + N12*S0bc + N13*S0cc;
349        M21= N21*S0aa + N22*S0ab + N23*S0ac;
350        M22= N21*S0ab + N22*S0bb + N23*S0bc;
351        M23= N21*S0ac + N22*S0bc + N23*S0cc;
352        M31= N31*S0aa + N32*S0ab + N33*S0ac;
353        M32= N31*S0ab + N32*S0bb + N33*S0bc;
354        M33= N31*S0ac + N32*S0bc + N33*S0cc;
355
356        DenQ1=1.0+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31;
357        DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33;
358        DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33;
359        DenQ=DenQ1+DenQ2+DenQ3;
360
361        Q11= (1.0 + M22-M23*M32 + M33 + M22*M33)/DenQ;
362        Q12= (-M12 + M13*M32 - M12*M33)/DenQ;
363        Q13= (-M13 - M13*M22 + M12*M23)/DenQ;
364        Q21= (-M21 + M23*M31 - M21*M33)/DenQ;
365        Q22= (1.0 + M11 - M13*M31 + M33 + M11*M33)/DenQ;
366        Q23= (M13*M21 - M23 - M11*M23)/DenQ;
367        Q31= (-M31 - M22*M31 + M21*M32)/DenQ;
368        Q32= (M12*M31 - M32 - M11*M32)/DenQ;
369        Q33= (1.0 + M11 - M12*M21 + M22 + M11*M22)/DenQ;
370
371        S11= Q11*S0aa + Q21*S0ab + Q31*S0ac;
372        S12= Q12*S0aa + Q22*S0ab + Q32*S0ac;
373        S13= Q13*S0aa + Q23*S0ab + Q33*S0ac;
374        S14=-S11-S12-S13;
375        S21= Q11*S0ba + Q21*S0bb + Q31*S0bc;
376        S22= Q12*S0ba + Q22*S0bb + Q32*S0bc;
377        S23= Q13*S0ba + Q23*S0bb + Q33*S0bc;
378        S24=-S21-S22-S23;
379        S31= Q11*S0ca + Q21*S0cb + Q31*S0cc;
380        S32= Q12*S0ca + Q22*S0cb + Q32*S0cc;
381        S33= Q13*S0ca + Q23*S0cb + Q33*S0cc;
382        S34=-S31-S32-S33;
383        S41=S14;
384        S42=S24;
385        S43=S34;
386        S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23;
387
388        Nav=6.022045e+23;
389        Lad=(La/va-Ld/vd)*sqrt(Nav);
390        Lbd=(Lb/vb-Ld/vd)*sqrt(Nav);
391        Lcd=(Lc/vc-Ld/vd)*sqrt(Nav);
392
393        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;
394
395        Intg *= scale;
396        Intg += background;
397
398    return Intg;
399
400}
401/**
402 * Function to evaluate 1D scattering function
403 * @param pars: parameters of the sphere
404 * @param q: q-value
405 * @return: function value
406 */
407double rpa_analytical_1D(RPAParameters *pars, double q) {
408        double dp[29];
409        //Fittable parameters
410        dp[0] = pars->lcase_n;
411
412        dp[1] = pars->ba;
413        dp[2] = pars->bb;
414        dp[3] = pars->bc;
415        dp[4] = pars->bd;
416
417        dp[5] = pars->Kab;
418        dp[6] = pars->Kac;
419        dp[7] = pars->Kad;
420        dp[8] = pars->Kbc;
421        dp[9] = pars->Kbd;
422        dp[10] = pars->Kcd;
423
424        dp[11] = pars->scale;
425        dp[12] = pars->background;
426
427        //Fixed parameters
428        dp[13] = pars->Na;
429        dp[14] = pars->Phia;
430        dp[15] = pars->va;
431        dp[16] = pars->La;
432
433        dp[17] = pars->Nb;
434        dp[18] = pars->Phib;
435        dp[19] = pars->vb;
436        dp[20] = pars->Lb;
437
438        dp[21] = pars->Nc;
439        dp[22] = pars->Phic;
440        dp[23] = pars->vc;
441        dp[24] = pars->Lc;
442
443        dp[25] = pars->Nd;
444        dp[26] = pars->Phid;
445        dp[27] = pars->vd;
446        dp[28] = pars->Ld;
447
448        return rpa_kernel(dp, q);
449}
450
451/**
452 * Function to evaluate 2D scattering function
453 * @param pars: parameters of the sphere
454 * @param q: q-value
455 * @return: function value
456 */
457double rpa_analytical_2D(RPAParameters *pars, double q, double phi) {
458        return rpa_analytical_1D(pars,q);
459}
460
461double rpa_analytical_2DXY(RPAParameters *pars, double qx, double qy){
462        return rpa_analytical_1D(pars,sqrt(qx*qx+qy*qy));
463}
Note: See TracBrowser for help on using the repository browser.