source: sasview/sansmodels/src/c_models/onion.cpp @ a8d3b4f

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 a8d3b4f was 5617775, checked in by Jae Cho <jhjcho@…>, 13 years ago

Onion model: fixed polydispersity problem
  • Property mode set to 100644
File size: 21.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 <stdlib.h>
26using namespace std;
27#include "onion.h"
28
29// Convenience parameter structure
30typedef struct {
31  int n_shells;
32  double scale;
33  double rad_core0;
34  double sld_core0;
35  double sld_solv;
36  double background;
37  double sld_out_shell1;
38  double sld_out_shell2;
39  double sld_out_shell3;
40  double sld_out_shell4;
41  double sld_out_shell5;
42  double sld_out_shell6;
43  double sld_out_shell7;
44  double sld_out_shell8;
45  double sld_out_shell9;
46  double sld_out_shell10;
47  double sld_in_shell1;
48  double sld_in_shell2;
49  double sld_in_shell3;
50  double sld_in_shell4;
51  double sld_in_shell5;
52  double sld_in_shell6;
53  double sld_in_shell7;
54  double sld_in_shell8;
55  double sld_in_shell9;
56  double sld_in_shell10;
57
58  double A_shell1;
59  double A_shell2;
60  double A_shell3;
61  double A_shell4;
62  double A_shell5;
63  double A_shell6;
64  double A_shell7;
65  double A_shell8;
66  double A_shell9;
67  double A_shell10;
68
69  double thick_shell1;
70  double thick_shell2;
71  double thick_shell3;
72  double thick_shell4;
73  double thick_shell5;
74  double thick_shell6;
75  double thick_shell7;
76  double thick_shell8;
77  double thick_shell9;
78  double thick_shell10;
79
80  int func_shell1;
81  int func_shell2;
82  int func_shell3;
83  int func_shell4;
84  int func_shell5;
85  int func_shell6;
86  int func_shell7;
87  int func_shell8;
88  int func_shell9;
89  int func_shell10;
90} OnionParameters;
91
92// some details can be found in sld_cal.c
93static double so_kernel(double dp[], double q) {
94  int n = dp[0];
95  double scale = dp[1];
96  double rad_core0 = dp[2];
97  double sld_core0 = dp[3];
98  double sld_solv = dp[4];
99  double background = dp[5];
100  int i,j;
101  double bes,fun,alpha,f,vol,vol_pre,vol_sub,qr,r,contr,f2;
102  double sign;
103  double pi;
104  double r0 = 0.0;
105
106  double *sld_out;
107  double *slope;
108  double *sld_in;
109  double *thick;
110  double *A;
111  int *fun_type;
112
113  sld_out = (double*)malloc((n+2)*sizeof(double));
114  slope = (double*)malloc((n+2)*sizeof(double));
115  sld_in = (double*)malloc((n+2)*sizeof(double));
116  thick = (double*)malloc((n+2)*sizeof(double));
117  A = (double*)malloc((n+2)*sizeof(double));
118  fun_type = (int*)malloc((n+2)*sizeof(int));
119
120  for (i =1; i<=n; i++){
121    sld_out[i] = dp[i+5];
122    sld_in[i] = dp[i+15];
123    A[i] = dp[i+25];
124    thick[i] = dp[i+35];
125    fun_type[i] = dp[i+45];
126  }
127  sld_out[0] = sld_core0;
128  sld_out[n+1] = sld_solv;
129  sld_in[0] = sld_core0;
130  sld_in[n+1] = sld_solv;
131  thick[0] = rad_core0;
132  thick[n+1] = 1e+10;
133  A[0] = 0.0;
134  A[n+1] = 0.0;
135  fun_type[0] = 0;
136  fun_type[n+1] = 0;
137
138
139  pi = 4.0*atan(1.0);
140  f = 0.0;
141  r = 0.0;
142  vol = 0.0;
143  vol_pre = 0.0;
144  vol_sub = 0.0;
145
146  for (i =0; i<= n+1; i++){
147    if (thick[i] == 0.0){
148      continue;
149    }
150    if (fun_type[i]== 0 ){
151      slope[i] = 0.0;
152      A[i] = 0.0;
153    }
154    vol_pre = vol;
155    switch(fun_type[i]){
156    case 2 :
157      r0 = r;
158      if (A[i] == 0.0){
159        slope[i] = 0.0;
160      }
161      else{
162        slope[i] = (sld_out[i]-sld_in[i])/(exp(A[i])-1.0);
163      }
164      for (j=0; j<2; j++){
165        if ( i == 0 && j == 0){
166          continue;
167        }
168        if (i == n+1 && j == 1){
169          continue;
170        }
171        if ( j == 1){
172          sign = 1.0;
173          r += thick[i];
174        }
175        else{
176          sign = -1.0;
177        }
178        qr = q * r;
179        alpha = A[i] * r/thick[i];
180        fun = 0.0;
181        if(qr == 0.0){
182          fun = sign * 1.0;
183          bes = sign * 1.0;
184        }
185        else{
186          if (fabs(A[i]) > 0.0 ){
187            fun = 3.0 * ((alpha*alpha - qr * qr) * sin(qr) - 2.0 * alpha * qr * cos(qr))/ ((alpha * alpha + qr * qr) * (alpha * alpha + qr * qr) * qr);
188            fun = fun - 3.0 * (alpha * sin(qr) - qr * cos(qr)) / ((alpha * alpha + qr * qr) * qr);
189            fun = - sign *fun;
190            bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
191          }
192          else {
193            fun = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
194            bes = sign * 3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
195          }
196        }
197        contr = slope[i]*exp(A[i]*(r-r0)/thick[i]);
198
199        vol = 4.0 * pi / 3.0 * r * r * r;
200        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && A[i]==0.0){
201        //  vol_sub += (vol_pre - vol);
202        //}
203        f += vol * (contr * (fun) + (sld_in[i]-slope[i]) * bes);
204      }
205      break;
206    default :
207      if (fun_type[i]==0){
208        slope[i] = 0.0;
209      }
210      else{
211        slope[i]= (sld_out[i] -sld_in[i])/thick[i];
212      }
213      contr = sld_in[i]-slope[i]*r;
214      for (j=0; j<2; j++){
215        if ( i == 0 && j == 0){
216          continue;
217        }
218        if (i == n+1 && j == 1){
219          continue;
220        }
221        if ( j == 1){
222          sign = 1.0;
223          r += thick[i];
224        }
225        else{
226          sign = -1.0;
227        }
228
229        qr = q * r;
230        fun = 0.0;
231        if(qr == 0.0){
232          bes = sign * 1.0;
233        }
234        else{
235          bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr);
236          if (fabs(slope[i]) > 0.0 ){
237            fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr);
238          }
239        }
240        vol = 4.0 * pi / 3.0 * r * r * r;
241        //if (j == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){
242        //  vol_sub += (vol_pre - vol);
243        //}
244        f += vol * (bes * contr + fun * slope[i]);
245      }
246      break;
247    }
248
249  }
250  //vol += vol_sub;
251  f2 = f * f / vol * 1.0e8;
252  f2 *= scale;
253  f2 += background;
254
255  free(sld_out);
256  free(slope);
257  free(sld_in);
258  free(thick);
259  free(A);
260  free(fun_type);
261
262  return (f2);
263}
264
265
266OnionModel :: OnionModel() {
267  n_shells = Parameter(1.0);
268  scale = Parameter(1.0);
269  rad_core0 = Parameter(200.0);
270  rad_core0.set_min(0.0);
271  sld_core0 = Parameter(1e-06);
272  sld_solv = Parameter(6.4e-06);
273  background = Parameter(0.0);
274
275  sld_out_shell1 = Parameter(1.0e-06);
276  sld_out_shell2 = Parameter(1.0e-06);
277  sld_out_shell3 = Parameter(1.0e-06);
278  sld_out_shell4 = Parameter(1.0e-06);
279  sld_out_shell5 = Parameter(1.0e-06);
280  sld_out_shell6 = Parameter(1.0e-06);
281  sld_out_shell7 = Parameter(1.0e-06);
282  sld_out_shell8 = Parameter(1.0e-06);
283  sld_out_shell9 = Parameter(1.0e-06);
284  sld_out_shell10 = Parameter(1.0e-06);
285
286  sld_in_shell1 = Parameter(2.3e-06);
287  sld_in_shell2 = Parameter(2.6e-06);
288  sld_in_shell3 = Parameter(2.9e-06);
289  sld_in_shell4 = Parameter(3.2e-06);
290  sld_in_shell5 = Parameter(3.5e-06);
291  sld_in_shell6 = Parameter(3.8e-06);
292  sld_in_shell7 = Parameter(4.1e-06);
293  sld_in_shell8 = Parameter(4.4e-06);
294  sld_in_shell9 = Parameter(4.7e-06);
295  sld_in_shell10 = Parameter(5.0e-06);
296
297  A_shell1 = Parameter(1.0);
298  A_shell2 = Parameter(1.0);
299  A_shell3 = Parameter(1.0);
300  A_shell4 = Parameter(1.0);
301  A_shell5 = Parameter(1.0);
302  A_shell6 = Parameter(1.0);
303  A_shell7 = Parameter(1.0);
304  A_shell8 = Parameter(1.0);
305  A_shell9 = Parameter(1.0);
306  A_shell10 = Parameter(1.0);
307
308  thick_shell1 = Parameter(50.0);
309  thick_shell1.set_min(0.0);
310  thick_shell2 = Parameter(50.0);
311  thick_shell2.set_min(0.0);
312  thick_shell3 = Parameter(50.0);
313  thick_shell3.set_min(0.0);
314  thick_shell4 = Parameter(50.0);
315  thick_shell4.set_min(0.0);
316  thick_shell5 = Parameter(50.0);
317  thick_shell5.set_min(0.0);
318  thick_shell6 = Parameter(50.0);
319  thick_shell6.set_min(0.0);
320  thick_shell7 = Parameter(50.0);
321  thick_shell7.set_min(0.0);
322  thick_shell8 = Parameter(50.0);
323  thick_shell8.set_min(0.0);
324  thick_shell9 = Parameter(50.0);
325  thick_shell9.set_min(0.0);
326  thick_shell10 = Parameter(50.0);
327  thick_shell10.set_min(0.0);
328
329  func_shell1 = Parameter(2);
330  func_shell2 = Parameter(2);
331  func_shell3 = Parameter(2);
332  func_shell4 = Parameter(2);
333  func_shell5 = Parameter(2);
334  func_shell6 = Parameter(2);
335  func_shell7 = Parameter(2);
336  func_shell8 = Parameter(2);
337  func_shell9 = Parameter(2);
338  func_shell10 = Parameter(2);
339}
340
341/**
342 * Function to evaluate 1D scattering function
343 * The NIST IGOR library is used for the actual calculation.
344 * @param q: q-value
345 * @return: function value
346 */
347double OnionModel :: operator()(double q) {
348  double dp[56];
349  // Fill parameter array for IGOR library
350  // Add the background after averaging
351  dp[0] = n_shells();
352  dp[1] = scale();
353  dp[2] = rad_core0();
354  dp[3] = sld_core0();
355  dp[4] = sld_solv();
356  dp[5] = 0.0;
357
358  dp[6] = sld_out_shell1();
359  dp[7] = sld_out_shell2();
360  dp[8] = sld_out_shell3();
361  dp[9] = sld_out_shell4();
362  dp[10] = sld_out_shell5();
363  dp[11] = sld_out_shell6();
364  dp[12] = sld_out_shell7();
365  dp[13] = sld_out_shell8();
366  dp[14] = sld_out_shell9();
367  dp[15] = sld_out_shell10();
368
369  dp[16] = sld_in_shell1();
370  dp[17] = sld_in_shell2();
371  dp[18] = sld_in_shell3();
372  dp[19] = sld_in_shell4();
373  dp[20] = sld_in_shell5();
374  dp[21] = sld_in_shell6();
375  dp[22] = sld_in_shell7();
376  dp[23] = sld_in_shell8();
377  dp[24] = sld_in_shell9();
378  dp[25] = sld_in_shell10();
379
380  dp[26] = A_shell1();
381  dp[27] = A_shell2();
382  dp[28] = A_shell3();
383  dp[29] = A_shell4();
384  dp[30] = A_shell5();
385  dp[31] = A_shell6();
386  dp[32] = A_shell7();
387  dp[33] = A_shell8();
388  dp[34] = A_shell9();
389  dp[35] = A_shell10();
390
391  dp[36] = thick_shell1();
392  dp[37] = thick_shell2();
393  dp[38] = thick_shell3();
394  dp[39] = thick_shell4();
395  dp[40] = thick_shell5();
396  dp[41] = thick_shell6();
397  dp[42] = thick_shell7();
398  dp[43] = thick_shell8();
399  dp[44] = thick_shell9();
400  dp[45] = thick_shell10();
401
402  dp[46] = func_shell1();
403  dp[47] = func_shell2();
404  dp[48] = func_shell3();
405  dp[49] = func_shell4();
406  dp[50] = func_shell5();
407  dp[51] = func_shell6();
408  dp[52] = func_shell7();
409  dp[53] = func_shell8();
410  dp[54] = func_shell9();
411  dp[55] = func_shell10();
412
413
414  // Get the dispersion points for the radius
415  vector<WeightPoint> weights_rad;
416  rad_core0.get_weights(weights_rad);
417
418  // Get the dispersion points for the thick 1
419  vector<WeightPoint> weights_s1;
420  thick_shell1.get_weights(weights_s1);
421
422  // Get the dispersion points for the thick 2
423  vector<WeightPoint> weights_s2;
424  thick_shell2.get_weights(weights_s2);
425
426  // Get the dispersion points for the thick 3
427  vector<WeightPoint> weights_s3;
428  thick_shell3.get_weights(weights_s3);
429
430  // Get the dispersion points for the thick 4
431  vector<WeightPoint> weights_s4;
432  thick_shell4.get_weights(weights_s4);
433
434  // Get the dispersion points for the thick 5
435  vector<WeightPoint> weights_s5;
436  thick_shell5.get_weights(weights_s5);
437
438  // Get the dispersion points for the thick 6
439  vector<WeightPoint> weights_s6;
440  thick_shell6.get_weights(weights_s6);
441
442  // Get the dispersion points for the thick 7
443  vector<WeightPoint> weights_s7;
444  thick_shell7.get_weights(weights_s7);
445
446  // Get the dispersion points for the thick 8
447  vector<WeightPoint> weights_s8;
448  thick_shell8.get_weights(weights_s8);
449  // Get the dispersion points for the thick 9
450  vector<WeightPoint> weights_s9;
451  thick_shell9.get_weights(weights_s9);
452
453  // Get the dispersion points for the thick 10
454  vector<WeightPoint> weights_s10;
455  thick_shell10.get_weights(weights_s10);
456
457
458  // Perform the computation, with all weight points
459  double sum = 0.0;
460  double norm = 0.0;
461  double vol = 0.0;
462
463  // Loop over radius weight points
464  for(size_t i=0; i<weights_rad.size(); i++) {
465    dp[2] = weights_rad[i].value;
466    // Loop over radius weight points
467    for(size_t j=0; j<weights_s1.size(); j++) {
468      dp[36] = weights_s1[j].value;
469      // Loop over radius weight points
470      for(size_t k=0; k<weights_s2.size(); k++) {
471        dp[37] = weights_s2[k].value;
472        // Loop over radius weight points
473        for(size_t l=0; l<weights_s3.size(); l++) {
474          dp[38] = weights_s3[l].value;
475          // Loop over radius weight points
476          for(size_t m=0; m<weights_s4.size(); m++) {
477            dp[39] = weights_s4[m].value;
478            for(size_t n=0; n<weights_s5.size(); n++) {
479              dp[40] = weights_s5[n].value;
480              for(size_t o=0; o<weights_s6.size(); o++) {
481                dp[41] = weights_s6[o].value;
482                for(size_t p=0; p<weights_s7.size(); p++) {
483                  dp[42] = weights_s7[p].value;
484                  for(size_t t=0; t<weights_s8.size(); t++) {
485                    dp[43] = weights_s8[t].value;
486                    for(size_t r=0; r<weights_s9.size(); r++) {
487                      dp[44] = weights_s9[r].value;
488                      for(size_t s=0; s<weights_s10.size(); s++) {
489                        dp[45] = weights_s10[s].value;
490                        //Un-normalize Shells by volume
491                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
492                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
493                            *weights_s9[r].weight*weights_s10[s].weight
494                            * so_kernel(dp,q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
495                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
496                                +weights_s9[r].value+weights_s10[s].value),3.0);
497                        //Find average volume
498                        vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
499                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
500                            *weights_s9[r].weight*weights_s10[s].weight
501                            * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
502                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
503                                +weights_s9[r].value+weights_s10[s].value),3.0);
504                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
505                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
506                            *weights_s9[r].weight*weights_s10[s].weight;
507                      }
508                    }
509                  }
510                }
511              }
512            }
513          }
514        }
515      }
516    }
517  }
518
519  if (vol != 0.0 && norm != 0.0) {
520    //Re-normalize by avg volume
521    sum = sum/(vol/norm);}
522
523  return sum/norm + background();
524}
525
526/**
527 * Function to evaluate 2D scattering function
528 * @param q_x: value of Q along x
529 * @param q_y: value of Q along y
530 * @return: function value
531 */
532double OnionModel :: operator()(double qx, double qy) {
533  double q = sqrt(qx*qx + qy*qy);
534  return (*this).operator()(q);
535}
536
537/**
538 * Function to evaluate 2D scattering function
539 * @param pars: parameters of the sphere
540 * @param q: q-value
541 * @param phi: angle phi
542 * @return: function value
543 */
544double OnionModel :: evaluate_rphi(double q, double phi) {
545  return (*this).operator()(q);
546}
547
548/**
549 * Function to calculate effective radius
550 * @return: effective radius value
551 */
552double OnionModel :: calculate_ER() {
553  OnionParameters dp;
554  dp.rad_core0 = rad_core0();
555  dp.thick_shell1 = thick_shell1();
556  dp.thick_shell2 = thick_shell2();
557  dp.thick_shell3 = thick_shell3();
558  dp.thick_shell4 = thick_shell4();
559  dp.thick_shell5 = thick_shell5();
560  dp.thick_shell6 = thick_shell6();
561  dp.thick_shell7 = thick_shell7();
562  dp.thick_shell8 = thick_shell8();
563  dp.thick_shell9 = thick_shell9();
564  dp.thick_shell10 = thick_shell10();
565
566  double rad_out = 0.0;
567  // Perform the computation, with all weight points
568  double sum = 0.0;
569  double norm = 0.0;
570
571  // Get the dispersion points for the radius
572  vector<WeightPoint> weights_rad;
573  rad_core0.get_weights(weights_rad);
574
575  // Get the dispersion points for the thick 1
576  vector<WeightPoint> weights_s1;
577  thick_shell1.get_weights(weights_s1);
578
579  // Get the dispersion points for the thick 2
580  vector<WeightPoint> weights_s2;
581  thick_shell2.get_weights(weights_s2);
582
583  // Get the dispersion points for the thick 3
584  vector<WeightPoint> weights_s3;
585  thick_shell3.get_weights(weights_s3);
586
587  // Get the dispersion points for the thick 4
588  vector<WeightPoint> weights_s4;
589  thick_shell4.get_weights(weights_s4);
590  // Get the dispersion points for the thick 5
591  vector<WeightPoint> weights_s5;
592  thick_shell5.get_weights(weights_s5);
593
594  // Get the dispersion points for the thick 6
595  vector<WeightPoint> weights_s6;
596  thick_shell6.get_weights(weights_s6);
597
598  // Get the dispersion points for the thick 7
599  vector<WeightPoint> weights_s7;
600  thick_shell7.get_weights(weights_s7);
601
602  // Get the dispersion points for the thick 8
603  vector<WeightPoint> weights_s8;
604  thick_shell8.get_weights(weights_s8);
605  // Get the dispersion points for the thick 9
606  vector<WeightPoint> weights_s9;
607  thick_shell9.get_weights(weights_s9);
608
609  // Get the dispersion points for the thick 10
610  vector<WeightPoint> weights_s10;
611  thick_shell10.get_weights(weights_s10);
612
613
614  // Loop over radius weight points
615  for(size_t i=0; i<weights_rad.size(); i++) {
616    dp.rad_core0 = weights_rad[i].value;
617    // Loop over radius weight points
618    for(size_t j=0; j<weights_s1.size(); j++) {
619      dp.thick_shell1 = weights_s1[j].value;
620      // Loop over radius weight points
621      for(size_t k=0; k<weights_s2.size(); k++) {
622        dp.thick_shell2 = weights_s2[k].value;
623        // Loop over radius weight points
624        for(size_t l=0; l<weights_s3.size(); l++) {
625          dp.thick_shell3 = weights_s3[l].value;
626          // Loop over radius weight points
627          for(size_t m=0; m<weights_s4.size(); m++) {
628            dp.thick_shell4 = weights_s4[m].value;
629            // Loop over radius weight points
630            for(size_t n=0; j<weights_s5.size(); n++) {
631              dp.thick_shell5 = weights_s5[n].value;
632              // Loop over radius weight points
633              for(size_t o=0; k<weights_s6.size(); o++) {
634                dp.thick_shell6 = weights_s6[o].value;
635                // Loop over radius weight points
636                for(size_t p=0; l<weights_s7.size(); p++) {
637                  dp.thick_shell7 = weights_s7[p].value;
638                  // Loop over radius weight points
639                  for(size_t t=0; m<weights_s8.size(); t++) {
640                    dp.thick_shell8 = weights_s8[t].value;
641                    // Loop over radius weight points
642                    for(size_t r=0; l<weights_s9.size(); r++) {
643                      dp.thick_shell8 = weights_s9[r].value;
644                      // Loop over radius weight points
645                      for(size_t s=0; m<weights_s10.size(); s++) {
646                        dp.thick_shell10 = weights_s10[s].value;
647                        //Un-normalize FourShell by volume
648                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
649                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
650                            *weights_s9[r].weight*weights_s10[s].weight
651                            * (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
652                                +dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
653                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight
654                            *weights_s4[m].weight*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight
655                            *weights_s8[t].weight*weights_s9[r].weight*weights_s10[s].weight;
656                      }
657                    }
658                  }
659                }
660              }
661            }
662          }
663        }
664      }
665    }
666  }
667
668  if (norm != 0){
669    //return the averaged value
670    rad_out =  sum/norm;}
671  else{
672    //return normal value
673    rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
674        +dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
675  return rad_out;
676}
677double OnionModel :: calculate_VR() {
678  return 1.0;
679}
Note: See TracBrowser for help on using the repository browser.