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

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

refactor onion model
  • Property mode set to 100644
File size: 20.6 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  sld_core0 = Parameter(1e-06);
271  sld_solv = Parameter(6.4e-06);
272  background = Parameter(0.0);
273
274  sld_out_shell1 = Parameter(1.0e-06);
275  sld_out_shell2 = Parameter(1.0e-06);
276  sld_out_shell3 = Parameter(1.0e-06);
277  sld_out_shell4 = Parameter(1.0e-06);
278  sld_out_shell5 = Parameter(1.0e-06);
279  sld_out_shell6 = Parameter(1.0e-06);
280  sld_out_shell7 = Parameter(1.0e-06);
281  sld_out_shell8 = Parameter(1.0e-06);
282  sld_out_shell9 = Parameter(1.0e-06);
283  sld_out_shell10 = Parameter(1.0e-06);
284
285  sld_in_shell1 = Parameter(2.3e-06);
286  sld_in_shell2 = Parameter(2.6e-06);
287  sld_in_shell3 = Parameter(2.9e-06);
288  sld_in_shell4 = Parameter(3.2e-06);
289  sld_in_shell5 = Parameter(3.5e-06);
290  sld_in_shell6 = Parameter(3.8e-06);
291  sld_in_shell7 = Parameter(4.1e-06);
292  sld_in_shell8 = Parameter(4.4e-06);
293  sld_in_shell9 = Parameter(4.7e-06);
294  sld_in_shell10 = Parameter(5.0e-06);
295
296  A_shell1 = Parameter(1.0);
297  A_shell2 = Parameter(1.0);
298  A_shell3 = Parameter(1.0);
299  A_shell4 = Parameter(1.0);
300  A_shell5 = Parameter(1.0);
301  A_shell6 = Parameter(1.0);
302  A_shell7 = Parameter(1.0);
303  A_shell8 = Parameter(1.0);
304  A_shell9 = Parameter(1.0);
305  A_shell10 = Parameter(1.0);
306
307  thick_shell1 = Parameter(50.0);
308  thick_shell2 = Parameter(50.0);
309  thick_shell3 = Parameter(50.0);
310  thick_shell4 = Parameter(50.0);
311  thick_shell5 = Parameter(50.0);
312  thick_shell6 = Parameter(50.0);
313  thick_shell7 = Parameter(50.0);
314  thick_shell8 = Parameter(50.0);
315  thick_shell9 = Parameter(50.0);
316  thick_shell10 = Parameter(50.0);
317
318  func_shell1 = Parameter(2);
319  func_shell2 = Parameter(2);
320  func_shell3 = Parameter(2);
321  func_shell4 = Parameter(2);
322  func_shell5 = Parameter(2);
323  func_shell6 = Parameter(2);
324  func_shell7 = Parameter(2);
325  func_shell8 = Parameter(2);
326  func_shell9 = Parameter(2);
327  func_shell10 = Parameter(2);
328}
329
330/**
331 * Function to evaluate 1D scattering function
332 * The NIST IGOR library is used for the actual calculation.
333 * @param q: q-value
334 * @return: function value
335 */
336double OnionModel :: operator()(double q) {
337  double dp[56];
338  // Fill parameter array for IGOR library
339  // Add the background after averaging
340  dp[0] = n_shells();
341  dp[1] = scale();
342  dp[2] = rad_core0();
343  dp[3] = sld_core0();
344  dp[4] = sld_solv();
345  dp[5] = 0.0;
346
347  dp[6] = sld_out_shell1();
348  dp[7] = sld_out_shell2();
349  dp[8] = sld_out_shell3();
350  dp[9] = sld_out_shell4();
351  dp[10] = sld_out_shell5();
352  dp[11] = sld_out_shell6();
353  dp[12] = sld_out_shell7();
354  dp[13] = sld_out_shell8();
355  dp[14] = sld_out_shell9();
356  dp[15] = sld_out_shell10();
357
358  dp[16] = sld_in_shell1();
359  dp[17] = sld_in_shell2();
360  dp[18] = sld_in_shell3();
361  dp[19] = sld_in_shell4();
362  dp[20] = sld_in_shell5();
363  dp[21] = sld_in_shell6();
364  dp[22] = sld_in_shell7();
365  dp[23] = sld_in_shell8();
366  dp[24] = sld_in_shell9();
367  dp[25] = sld_in_shell10();
368
369  dp[26] = A_shell1();
370  dp[27] = A_shell2();
371  dp[28] = A_shell3();
372  dp[29] = A_shell4();
373  dp[30] = A_shell5();
374  dp[31] = A_shell6();
375  dp[32] = A_shell7();
376  dp[33] = A_shell8();
377  dp[34] = A_shell9();
378  dp[35] = A_shell10();
379
380  dp[36] = thick_shell1();
381  dp[37] = thick_shell2();
382  dp[38] = thick_shell3();
383  dp[39] = thick_shell4();
384  dp[40] = thick_shell5();
385  dp[41] = thick_shell6();
386  dp[42] = thick_shell7();
387  dp[43] = thick_shell8();
388  dp[44] = thick_shell9();
389  dp[45] = thick_shell10();
390
391  dp[46] = func_shell1();
392  dp[47] = func_shell2();
393  dp[48] = func_shell3();
394  dp[49] = func_shell4();
395  dp[50] = func_shell5();
396  dp[51] = func_shell6();
397  dp[52] = func_shell7();
398  dp[53] = func_shell8();
399  dp[54] = func_shell9();
400  dp[55] = func_shell10();
401
402
403  // Get the dispersion points for the radius
404  vector<WeightPoint> weights_rad;
405  rad_core0.get_weights(weights_rad);
406
407  // Get the dispersion points for the thick 1
408  vector<WeightPoint> weights_s1;
409  thick_shell1.get_weights(weights_s1);
410
411  // Get the dispersion points for the thick 2
412  vector<WeightPoint> weights_s2;
413  thick_shell2.get_weights(weights_s2);
414
415  // Get the dispersion points for the thick 3
416  vector<WeightPoint> weights_s3;
417  thick_shell3.get_weights(weights_s3);
418
419  // Get the dispersion points for the thick 4
420  vector<WeightPoint> weights_s4;
421  thick_shell4.get_weights(weights_s4);
422
423  // Get the dispersion points for the thick 5
424  vector<WeightPoint> weights_s5;
425  thick_shell5.get_weights(weights_s5);
426
427  // Get the dispersion points for the thick 6
428  vector<WeightPoint> weights_s6;
429  thick_shell6.get_weights(weights_s6);
430
431  // Get the dispersion points for the thick 7
432  vector<WeightPoint> weights_s7;
433  thick_shell7.get_weights(weights_s7);
434
435  // Get the dispersion points for the thick 8
436  vector<WeightPoint> weights_s8;
437  thick_shell8.get_weights(weights_s8);
438  // Get the dispersion points for the thick 9
439  vector<WeightPoint> weights_s9;
440  thick_shell9.get_weights(weights_s9);
441
442  // Get the dispersion points for the thick 10
443  vector<WeightPoint> weights_s10;
444  thick_shell10.get_weights(weights_s10);
445
446
447  // Perform the computation, with all weight points
448  double sum = 0.0;
449  double norm = 0.0;
450  double vol = 0.0;
451
452  // Loop over radius weight points
453  for(size_t i=0; i<weights_rad.size(); i++) {
454    dp[2] = weights_rad[i].value;
455    // Loop over radius weight points
456    for(size_t j=0; j<weights_s1.size(); j++) {
457      dp[36] = weights_s1[j].value;
458      // Loop over radius weight points
459      for(size_t k=0; k<weights_s2.size(); k++) {
460        dp[37] = weights_s2[k].value;
461        // Loop over radius weight points
462        for(size_t l=0; l<weights_s3.size(); l++) {
463          dp[38] = weights_s3[l].value;
464          // Loop over radius weight points
465          for(size_t m=0; m<weights_s4.size(); m++) {
466            dp[39] = weights_s4[m].value;
467            for(size_t n=0; n<weights_s5.size(); n++) {
468              dp[40] = weights_s5[n].value;
469              for(size_t o=0; o<weights_s6.size(); o++) {
470                dp[41] = weights_s6[o].value;
471                for(size_t p=0; p<weights_s7.size(); p++) {
472                  dp[42] = weights_s7[p].value;
473                  for(size_t t=0; t<weights_s8.size(); t++) {
474                    dp[43] = weights_s8[t].value;
475                    for(size_t r=0; r<weights_s9.size(); r++) {
476                      dp[44] = weights_s9[r].value;
477                      for(size_t s=0; s<weights_s10.size(); s++) {
478                        dp[45] = weights_s10[s].value;
479                        //Un-normalize Shells by volume
480                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
481                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
482                            *weights_s9[r].weight*weights_s10[s].weight
483                            * 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
484                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
485                                +weights_s9[r].value+weights_s10[s].value),3.0);
486                        //Find average volume
487                        vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
488                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
489                            *weights_s9[r].weight*weights_s10[s].weight
490                            * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
491                                +weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
492                                +weights_s9[r].value+weights_s10[s].value),3.0);
493                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
494                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
495                            *weights_s9[r].weight*weights_s10[s].weight;
496                      }
497                    }
498                  }
499                }
500              }
501            }
502          }
503        }
504      }
505    }
506  }
507
508  if (vol != 0.0 && norm != 0.0) {
509    //Re-normalize by avg volume
510    sum = sum/(vol/norm);}
511
512  return sum/norm + background();
513}
514
515/**
516 * Function to evaluate 2D scattering function
517 * @param q_x: value of Q along x
518 * @param q_y: value of Q along y
519 * @return: function value
520 */
521double OnionModel :: operator()(double qx, double qy) {
522  double q = sqrt(qx*qx + qy*qy);
523  return (*this).operator()(q);
524}
525
526/**
527 * Function to evaluate 2D scattering function
528 * @param pars: parameters of the sphere
529 * @param q: q-value
530 * @param phi: angle phi
531 * @return: function value
532 */
533double OnionModel :: evaluate_rphi(double q, double phi) {
534  return (*this).operator()(q);
535}
536
537/**
538 * Function to calculate effective radius
539 * @return: effective radius value
540 */
541double OnionModel :: calculate_ER() {
542  OnionParameters dp;
543  dp.rad_core0 = rad_core0();
544  dp.thick_shell1 = thick_shell1();
545  dp.thick_shell2 = thick_shell2();
546  dp.thick_shell3 = thick_shell3();
547  dp.thick_shell4 = thick_shell4();
548  dp.thick_shell5 = thick_shell5();
549  dp.thick_shell6 = thick_shell6();
550  dp.thick_shell7 = thick_shell7();
551  dp.thick_shell8 = thick_shell8();
552  dp.thick_shell9 = thick_shell9();
553  dp.thick_shell10 = thick_shell10();
554
555  double rad_out = 0.0;
556  // Perform the computation, with all weight points
557  double sum = 0.0;
558  double norm = 0.0;
559
560  // Get the dispersion points for the radius
561  vector<WeightPoint> weights_rad;
562  rad_core0.get_weights(weights_rad);
563
564  // Get the dispersion points for the thick 1
565  vector<WeightPoint> weights_s1;
566  thick_shell1.get_weights(weights_s1);
567
568  // Get the dispersion points for the thick 2
569  vector<WeightPoint> weights_s2;
570  thick_shell2.get_weights(weights_s2);
571
572  // Get the dispersion points for the thick 3
573  vector<WeightPoint> weights_s3;
574  thick_shell3.get_weights(weights_s3);
575
576  // Get the dispersion points for the thick 4
577  vector<WeightPoint> weights_s4;
578  thick_shell4.get_weights(weights_s4);
579  // Get the dispersion points for the thick 5
580  vector<WeightPoint> weights_s5;
581  thick_shell5.get_weights(weights_s5);
582
583  // Get the dispersion points for the thick 6
584  vector<WeightPoint> weights_s6;
585  thick_shell6.get_weights(weights_s6);
586
587  // Get the dispersion points for the thick 7
588  vector<WeightPoint> weights_s7;
589  thick_shell7.get_weights(weights_s7);
590
591  // Get the dispersion points for the thick 8
592  vector<WeightPoint> weights_s8;
593  thick_shell8.get_weights(weights_s8);
594  // Get the dispersion points for the thick 9
595  vector<WeightPoint> weights_s9;
596  thick_shell9.get_weights(weights_s9);
597
598  // Get the dispersion points for the thick 10
599  vector<WeightPoint> weights_s10;
600  thick_shell10.get_weights(weights_s10);
601
602
603  // Loop over radius weight points
604  for(size_t i=0; i<weights_rad.size(); i++) {
605    dp.rad_core0 = weights_rad[i].value;
606    // Loop over radius weight points
607    for(size_t j=0; j<weights_s1.size(); j++) {
608      dp.thick_shell1 = weights_s1[j].value;
609      // Loop over radius weight points
610      for(size_t k=0; k<weights_s2.size(); k++) {
611        dp.thick_shell2 = weights_s2[k].value;
612        // Loop over radius weight points
613        for(size_t l=0; l<weights_s3.size(); l++) {
614          dp.thick_shell3 = weights_s3[l].value;
615          // Loop over radius weight points
616          for(size_t m=0; m<weights_s4.size(); m++) {
617            dp.thick_shell4 = weights_s4[m].value;
618            // Loop over radius weight points
619            for(size_t n=0; j<weights_s5.size(); n++) {
620              dp.thick_shell5 = weights_s5[n].value;
621              // Loop over radius weight points
622              for(size_t o=0; k<weights_s6.size(); o++) {
623                dp.thick_shell6 = weights_s6[o].value;
624                // Loop over radius weight points
625                for(size_t p=0; l<weights_s7.size(); p++) {
626                  dp.thick_shell7 = weights_s7[p].value;
627                  // Loop over radius weight points
628                  for(size_t t=0; m<weights_s8.size(); t++) {
629                    dp.thick_shell8 = weights_s8[t].value;
630                    // Loop over radius weight points
631                    for(size_t r=0; l<weights_s9.size(); r++) {
632                      dp.thick_shell8 = weights_s9[r].value;
633                      // Loop over radius weight points
634                      for(size_t s=0; m<weights_s10.size(); s++) {
635                        dp.thick_shell10 = weights_s10[s].value;
636                        //Un-normalize FourShell by volume
637                        sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
638                            *weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight*weights_s8[t].weight
639                            *weights_s9[r].weight*weights_s10[s].weight
640                            * (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
641                                +dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
642                        norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight
643                            *weights_s4[m].weight*weights_s5[n].weight*weights_s6[o].weight*weights_s7[p].weight
644                            *weights_s8[t].weight*weights_s9[r].weight*weights_s10[s].weight;
645                      }
646                    }
647                  }
648                }
649              }
650            }
651          }
652        }
653      }
654    }
655  }
656
657  if (norm != 0){
658    //return the averaged value
659    rad_out =  sum/norm;}
660  else{
661    //return normal value
662    rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
663        +dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
664  return rad_out;
665}
Note: See TracBrowser for help on using the repository browser.