source: sasview/sansmodels/src/c_models/csparallelepiped.cpp @ df88829

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 df88829 was df88829, checked in by Mathieu Doucet <doucetm@…>, 12 years ago

refactor csparallelepiped
  • Property mode set to 100644
File size: 15.9 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 2010, 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#include <math.h>
22#include "models.hh"
23#include "parameters.hh"
24#include <stdio.h>
25using namespace std;
26
27extern "C" {
28        #include "libCylinder.h"
29        #include "libStructureFactor.h"
30        #include "csparallelepiped.h"
31}
32
33// Convenience parameter structure
34typedef struct {
35    double scale;
36    double shortA;
37    double midB;
38    double longC;
39    double rimA;
40    double rimB;
41    double rimC;
42    double sld_rimA;
43    double sld_rimB;
44    double sld_rimC;
45    double sld_pcore;
46    double sld_solv;
47    double background;
48    double parallel_theta;
49    double parallel_phi;
50    double parallel_psi;
51} CSParallelepipedParameters;
52
53static double cspkernel(double dp[],double q, double ala, double alb, double alc){
54    // mu passed in is really mu*sqrt(1-sig^2)
55    double argA,argB,argC,argtA,argtB,argtC,tmp1,tmp2,tmp3,tmpt1,tmpt2,tmpt3;     //local variables
56
57  double aa,bb,cc, ta,tb,tc;
58  double Vin,Vot,V1,V2,V3;
59  double rhoA,rhoB,rhoC, rhoP, rhosolv;
60  double dr0, drA,drB, drC;
61  double retVal;
62
63  aa = dp[1];
64  bb = dp[2];
65  cc = dp[3];
66  ta = dp[4];
67  tb = dp[5];
68  tc = dp[6];
69  rhoA=dp[7];
70  rhoB=dp[8];
71  rhoC=dp[9];
72  rhoP=dp[10];
73  rhosolv=dp[11];
74  dr0=rhoP-rhosolv;
75  drA=rhoA-rhosolv;
76  drB=rhoB-rhosolv;
77  drC=rhoC-rhosolv;
78  Vin=(aa*bb*cc);
79  Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
80  V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
81  V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
82  V3=(2.0*aa*bb*tc);
83  //aa = aa/bb;
84  ta=(aa+2.0*ta);///bb;
85  tb=(aa+2.0*tb);///bb;
86  tc=(aa+2.0*tc);
87    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
88    argA = q*aa*ala/2.0;
89    argB = q*bb*alb/2.0;
90    argC = q*cc*alc/2.0;
91    argtA = q*ta*ala/2.0;
92  argtB = q*tb*alb/2.0;
93  argtC = q*tc*alc/2.0;
94
95    if(argA==0.0) {
96    tmp1 = 1.0;
97  } else {
98    tmp1 = sin(argA)/argA;
99    }
100    if (argB==0.0) {
101    tmp2 = 1.0;
102  } else {
103    tmp2 = sin(argB)/argB;
104    }
105
106    if (argC==0.0) {
107    tmp3 = 1.0;
108  } else {
109    tmp3 = sin(argC)/argC;
110    }
111    if(argtA==0.0) {
112    tmpt1 = 1.0;
113  } else {
114    tmpt1 = sin(argtA)/argtA;
115    }
116    if (argtB==0.0) {
117    tmpt2 = 1.0;
118  } else {
119    tmpt2 = sin(argtB)/argtB;
120    }
121    if (argtC==0.0) {
122    tmpt3 = 1.0;
123  } else {
124    tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC;
125    }
126    // This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
127    retVal =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
128        ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
129    //retVal *= (tmp3*tmp3);
130    retVal /= Vot;
131
132    return (retVal);
133
134}//Function cspkernel()
135
136/**
137 * Function to evaluate 2D scattering function
138 * @param pars: parameters of the CSparallelepiped
139 * @param q: q-value
140 * @param q_x: q_x / q
141 * @param q_y: q_y / q
142 * @return: function value
143 */
144static double csparallelepiped_analytical_2D_scaled(CSParallelepipedParameters *pars, double q, double q_x, double q_y) {
145  double dp[13];
146  double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y;
147  double q_z;
148  double alpha, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
149
150  double answer;
151  //convert angle degree to radian
152  double pi = 4.0*atan(1.0);
153  double theta = pars->parallel_theta * pi/180.0;
154  double phi = pars->parallel_phi * pi/180.0;
155  double psi = pars->parallel_psi* pi/180.0;
156
157  // Fill paramater array
158  dp[0] = 1.0;
159  dp[1] = pars->shortA;
160  dp[2] = pars->midB;
161  dp[3] = pars->longC;
162  dp[4] = pars->rimA;
163  dp[5] = pars->rimB;
164  dp[6] = pars->rimC;
165  dp[7] = pars->sld_rimA;
166  dp[8] = pars->sld_rimB;
167  dp[9] = pars->sld_rimC;
168  dp[10] = pars->sld_pcore;
169  dp[11] = pars->sld_solv;
170  dp[12] = 0.0;
171
172
173  edgeA = pars->shortA;
174  edgeB = pars->midB;
175  edgeC = pars->longC;
176
177
178    // parallelepiped c axis orientation
179    cparallel_x = sin(theta) * cos(phi);
180    cparallel_y = sin(theta) * sin(phi);
181    cparallel_z = cos(theta);
182
183    // q vector
184    q_z = 0.0;
185
186    // Compute the angle btw vector q and the
187    // axis of the parallelepiped
188    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
189    alpha = acos(cos_val_c);
190
191    // parallelepiped a axis orientation
192    parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
193    parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
194
195    cos_val_a = parallel_x*q_x + parallel_y*q_y;
196
197
198
199    // parallelepiped b axis orientation
200    bparallel_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
201    bparallel_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
202    // axis of the parallelepiped
203    cos_val_b = sin(acos(cos_val_a)) ;
204
205
206
207    // The following test should always pass
208    if (fabs(cos_val_c)>1.0) {
209      printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
210      return 0;
211    }
212
213  // Call the IGOR library function to get the kernel
214  answer = cspkernel( dp,q, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
215
216  //convert to [cm-1]
217  answer *= 1.0e8;
218
219  //Scale
220  answer *= pars->scale;
221
222  // add in the background
223  answer += pars->background;
224
225  return answer;
226}
227
228/**
229 * Function to evaluate 2D scattering function
230 * @param pars: parameters of the CSparallelepiped
231 * @param q: q-value
232 * @return: function value
233 */
234static double csparallelepiped_analytical_2DXY(CSParallelepipedParameters *pars, double qx, double qy) {
235  double q;
236  q = sqrt(qx*qx+qy*qy);
237    return csparallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
238}
239
240
241
242
243CSParallelepipedModel :: CSParallelepipedModel() {
244        scale      = Parameter(1.0);
245        shortA     = Parameter(35.0, true);
246        shortA.set_min(1.0);
247        midB     = Parameter(75.0, true);
248        midB.set_min(1.0);
249        longC    = Parameter(400.0, true);
250        longC.set_min(1.0);
251        rimA     = Parameter(10.0, true);
252        rimB     = Parameter(10.0, true);
253        rimC     = Parameter(10.0, true);
254        sld_rimA     = Parameter(2.0e-6, true);
255        sld_rimB     = Parameter(4.0e-6, true);
256        sld_rimC    = Parameter(2.0e-6, true);
257        sld_pcore   = Parameter(1.0e-6);
258        sld_solv   = Parameter(6.0e-6);
259        background = Parameter(0.06);
260        parallel_theta  = Parameter(0.0, true);
261        parallel_phi    = Parameter(0.0, true);
262        parallel_psi    = Parameter(0.0, true);
263}
264
265/**
266 * Function to evaluate 1D scattering function
267 * The NIST IGOR library is used for the actual calculation.
268 * @param q: q-value
269 * @return: function value
270 */
271double CSParallelepipedModel :: operator()(double q) {
272        double dp[13];
273
274        // Fill parameter array for IGOR library
275        // Add the background after averaging
276        dp[0] = scale();
277        dp[1] = shortA();
278        dp[2] = midB();
279        dp[3] = longC();
280        dp[4] = rimA();
281        dp[5] = rimB();
282        dp[6] = rimC();
283        dp[7] = sld_rimA();
284        dp[8] = sld_rimB();
285        dp[9] = sld_rimC();
286        dp[10] = sld_pcore();
287        dp[11] = sld_solv();
288        dp[12] = 0.0;
289
290        // Get the dispersion points for the short_edgeA
291        vector<WeightPoint> weights_shortA;
292        shortA.get_weights(weights_shortA);
293
294        // Get the dispersion points for the longer_edgeB
295        vector<WeightPoint> weights_midB;
296        midB.get_weights(weights_midB);
297
298        // Get the dispersion points for the longuest_edgeC
299        vector<WeightPoint> weights_longC;
300        longC.get_weights(weights_longC);
301
302
303
304        // Perform the computation, with all weight points
305        double sum = 0.0;
306        double norm = 0.0;
307        double vol = 0.0;
308
309        // Loop over short_edgeA weight points
310        for(int i=0; i< (int)weights_shortA.size(); i++) {
311                dp[1] = weights_shortA[i].value;
312
313                // Loop over longer_edgeB weight points
314                for(int j=0; j< (int)weights_midB.size(); j++) {
315                        dp[2] = weights_midB[j].value;
316
317                        // Loop over longuest_edgeC weight points
318                        for(int k=0; k< (int)weights_longC.size(); k++) {
319                                dp[3] = weights_longC[k].value;
320                                //Un-normalize  by volume
321                                sum += weights_shortA[i].weight * weights_midB[j].weight
322                                        * weights_longC[k].weight * CSParallelepiped(dp, q)
323                                        * weights_shortA[i].value*weights_midB[j].value
324                                        * weights_longC[k].value;
325                                //Find average volume
326                                vol += weights_shortA[i].weight * weights_midB[j].weight
327                                        * weights_longC[k].weight
328                                        * weights_shortA[i].value * weights_midB[j].value
329                                        * weights_longC[k].value;
330
331                                norm += weights_shortA[i].weight
332                                         * weights_midB[j].weight * weights_longC[k].weight;
333                        }
334                }
335        }
336        if (vol != 0.0 && norm != 0.0) {
337                //Re-normalize by avg volume
338                sum = sum/(vol/norm);}
339
340        return sum/norm + background();
341}
342/**
343 * Function to evaluate 2D scattering function
344 * @param q_x: value of Q along x
345 * @param q_y: value of Q along y
346 * @return: function value
347 */
348double CSParallelepipedModel :: operator()(double qx, double qy) {
349        CSParallelepipedParameters dp;
350        // Fill parameter array
351        dp.scale      = scale();
352        dp.shortA   = shortA();
353        dp.midB   = midB();
354        dp.longC  = longC();
355        dp.rimA   = rimA();
356        dp.rimB   = rimB();
357        dp.rimC  = rimC();
358        dp.sld_rimA   = sld_rimA();
359        dp.sld_rimB   = sld_rimB();
360        dp.sld_rimC  = sld_rimC();
361        dp.sld_pcore   = sld_pcore();
362        dp.sld_solv   = sld_solv();
363        dp.background = 0.0;
364        //dp.background = background();
365        dp.parallel_theta  = parallel_theta();
366        dp.parallel_phi    = parallel_phi();
367        dp.parallel_psi    = parallel_psi();
368
369
370
371        // Get the dispersion points for the short_edgeA
372        vector<WeightPoint> weights_shortA;
373        shortA.get_weights(weights_shortA);
374
375        // Get the dispersion points for the longer_edgeB
376        vector<WeightPoint> weights_midB;
377        midB.get_weights(weights_midB);
378
379        // Get the dispersion points for the longuest_edgeC
380        vector<WeightPoint> weights_longC;
381        longC.get_weights(weights_longC);
382
383        // Get angular averaging for theta
384        vector<WeightPoint> weights_parallel_theta;
385        parallel_theta.get_weights(weights_parallel_theta);
386
387        // Get angular averaging for phi
388        vector<WeightPoint> weights_parallel_phi;
389        parallel_phi.get_weights(weights_parallel_phi);
390
391        // Get angular averaging for psi
392        vector<WeightPoint> weights_parallel_psi;
393        parallel_psi.get_weights(weights_parallel_psi);
394
395        // Perform the computation, with all weight points
396        double sum = 0.0;
397        double norm = 0.0;
398        double norm_vol = 0.0;
399        double vol = 0.0;
400        double pi = 4.0*atan(1.0);
401
402        // Loop over radius weight points
403        for(int i=0; i< (int)weights_shortA.size(); i++) {
404                dp.shortA = weights_shortA[i].value;
405
406                // Loop over longer_edgeB weight points
407                for(int j=0; j< (int)weights_midB.size(); j++) {
408                        dp.midB = weights_midB[j].value;
409
410                        // Average over longuest_edgeC distribution
411                        for(int k=0; k< (int)weights_longC.size(); k++) {
412                                dp.longC = weights_longC[k].value;
413
414                                // Average over theta distribution
415                                for(int l=0; l< (int)weights_parallel_theta.size(); l++) {
416                                dp.parallel_theta = weights_parallel_theta[l].value;
417
418                                        // Average over phi distribution
419                                        for(int m=0; m< (int)weights_parallel_phi.size(); m++) {
420                                                dp.parallel_phi = weights_parallel_phi[m].value;
421
422                                                // Average over phi distribution
423                                                for(int n=0; n< (int)weights_parallel_psi.size(); n++) {
424                                                        dp.parallel_psi = weights_parallel_psi[n].value;
425                                                        //Un-normalize by volume
426                                                        double _ptvalue = weights_shortA[i].weight
427                                                                * weights_midB[j].weight
428                                                                * weights_longC[k].weight
429                                                                * weights_parallel_theta[l].weight
430                                                                * weights_parallel_phi[m].weight
431                                                                * weights_parallel_psi[n].weight
432                                                                * csparallelepiped_analytical_2DXY(&dp, qx, qy)
433                                                                * weights_shortA[i].value*weights_midB[j].value
434                                                                * weights_longC[k].value;
435
436                                                        if (weights_parallel_theta.size()>1) {
437                                                                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
438                                                        }
439                                                        sum += _ptvalue;
440                                                        //Find average volume
441                                                        vol += weights_shortA[i].weight
442                                                                * weights_midB[j].weight
443                                                                * weights_longC[k].weight
444                                                                * weights_shortA[i].value*weights_midB[j].value
445                                                                * weights_longC[k].value;
446                                                        //Find norm for volume
447                                                        norm_vol += weights_shortA[i].weight
448                                                                * weights_midB[j].weight
449                                                                * weights_longC[k].weight;
450
451                                                        norm += weights_shortA[i].weight
452                                                                * weights_midB[j].weight
453                                                                * weights_longC[k].weight
454                                                                * weights_parallel_theta[l].weight
455                                                                * weights_parallel_phi[m].weight
456                                                                * weights_parallel_psi[n].weight;
457                                                }
458                                        }
459
460                                }
461                        }
462                }
463        }
464        // Averaging in theta needs an extra normalization
465        // factor to account for the sin(theta) term in the
466        // integration (see documentation).
467        if (weights_parallel_theta.size()>1) norm = norm / asin(1.0);
468
469        if (vol != 0.0 && norm_vol != 0.0) {
470                //Re-normalize by avg volume
471                sum = sum/(vol/norm_vol);}
472
473        return sum/norm + background();
474}
475
476
477/**
478 * Function to evaluate 2D scattering function
479 * @param pars: parameters of the cylinder
480 * @param q: q-value
481 * @param phi: angle phi
482 * @return: function value
483 */
484double CSParallelepipedModel :: evaluate_rphi(double q, double phi) {
485        double qx = q*cos(phi);
486        double qy = q*sin(phi);
487        return (*this).operator()(qx, qy);
488}
489/**
490 * Function to calculate effective radius
491 * @return: effective radius value
492 */
493double CSParallelepipedModel :: calculate_ER() {
494        CSParallelepipedParameters dp;
495        dp.shortA   = shortA();
496        dp.midB   = midB();
497        dp.longC  = longC();
498        dp.rimA   = rimA();
499        dp.rimB   = rimB();
500        dp.rimC  = rimC();
501
502        double rad_out = 0.0;
503        double pi = 4.0*atan(1.0);
504        double suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
505        double height =(dp.longC + 2.0*dp.rimC);
506        // Perform the computation, with all weight points
507        double sum = 0.0;
508        double norm = 0.0;
509
510        // Get the dispersion points for the short_edgeA
511        vector<WeightPoint> weights_shortA;
512        shortA.get_weights(weights_shortA);
513
514        // Get the dispersion points for the longer_edgeB
515        vector<WeightPoint> weights_midB;
516        midB.get_weights(weights_midB);
517
518        // Get angular averaging for the longuest_edgeC
519        vector<WeightPoint> weights_longC;
520        longC.get_weights(weights_longC);
521
522        // Loop over radius weight points
523        for(int i=0; i< (int)weights_shortA.size(); i++) {
524                dp.shortA = weights_shortA[i].value;
525
526                // Loop over longer_edgeB weight points
527                for(int j=0; j< (int)weights_midB.size(); j++) {
528                        dp.midB = weights_midB[j].value;
529
530                        // Average over longuest_edgeC distribution
531                        for(int k=0; k< (int)weights_longC.size(); k++) {
532                                dp.longC = weights_longC[k].value;
533                                //Calculate surface averaged radius
534                                //This is rough approximation.
535                                suf_rad = sqrt((dp.shortA*dp.midB+2.0*dp.rimA*dp.midB+2.0*dp.rimA*dp.shortA)/pi);
536                                height =(dp.longC + 2.0*dp.rimC);
537                                //Note: output of "DiamCyl(dp.length,dp.radius)" is a DIAMETER.
538                                sum +=weights_shortA[i].weight* weights_midB[j].weight
539                                        * weights_longC[k].weight*DiamCyl(height, suf_rad)/2.0;
540                                norm += weights_shortA[i].weight* weights_midB[j].weight*weights_longC[k].weight;
541                        }
542                }
543        }
544
545        if (norm != 0){
546                //return the averaged value
547                rad_out =  sum/norm;}
548        else{
549                //return normal value
550                //Note: output of "DiamCyl(length,radius)" is DIAMETER.
551                rad_out = DiamCyl(dp.longC, suf_rad)/2.0;}
552        return rad_out;
553
554}
Note: See TracBrowser for help on using the repository browser.