Changeset 37805e9 in sasview for sansmodels


Ignore:
Timestamp:
Jan 5, 2012 11:14:43 AM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
6e10cff
Parents:
a8eab1c
Message:

spheresld refactor

Location:
sansmodels/src
Files:
2 deleted
4 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_extensions/spheresld.h

    r67424cd r37805e9  
    11#if !defined(o_h) 
    22#define sphere_sld_h 
     3#include "parameters.hh" 
    34 
    45/** 
    … …  
    3031 //[ORIENTATION_PARAMS]= <text> </text> 
    3132 
    32 typedef struct { 
    33         /// number of shells 
    34         //  [DEFAULT]=n_shells=1 
    35         int n_shells; 
     33class SphereSLDModel{ 
     34public: 
     35  // Model parameters 
     36  /// number of shells 
     37  //  [DEFAULT]=n_shells=1 
     38  Parameter n_shells; 
    3639    /// Scale factor 
    3740    //  [DEFAULT]=scale= 1.0 
    38         double scale; 
    39     /// thick_inter0 [A] 
     41  Parameter scale; 
     42    /// thick_inter0 [A] 
    4043    //  [DEFAULT]=thick_inter0=50.0 [A] 
    41         double thick_inter0; 
    42         ///     func_inter0 
    43         //  [DEFAULT]=func_inter0= 0 
    44         double func_inter0; 
    45         ///     sld_core0 [1/A^(2)] 
    46         //  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)] 
    47         double sld_core0; 
    48         ///     sld_solv [1/A^(2)] 
    49         //  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)] 
    50         double sld_solv; 
    51         /// Background 
    52         //  [DEFAULT]=background=0 
    53         double background; 
     44  Parameter thick_inter0; 
     45  /// func_inter0 
     46  //  [DEFAULT]=func_inter0= 0 
     47  Parameter func_inter0; 
     48  /// sld_core0 [1/A^(2)] 
     49  //  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)] 
     50  Parameter sld_core0; 
     51  /// sld_solv [1/A^(2)] 
     52  //  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)] 
     53  Parameter sld_solv; 
     54  /// Background 
     55  //  [DEFAULT]=background=0 
     56  Parameter background; 
    5457 
    5558    //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)] 
    56     double sld_flat1; 
     59    Parameter sld_flat1; 
    5760    //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)] 
    58     double sld_flat2; 
     61    Parameter sld_flat2; 
    5962    //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)] 
    60     double sld_flat3; 
     63    Parameter sld_flat3; 
    6164    //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)] 
    62     double sld_flat4; 
     65    Parameter sld_flat4; 
    6366    //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)] 
    64     double sld_flat5; 
     67    Parameter sld_flat5; 
    6568    //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)] 
    66     double sld_flat6; 
     69    Parameter sld_flat6; 
    6770    //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)] 
    68     double sld_flat7; 
     71    Parameter sld_flat7; 
    6972    //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)] 
    70     double sld_flat8; 
     73    Parameter sld_flat8; 
    7174    //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)] 
    72     double sld_flat9; 
     75    Parameter sld_flat9; 
    7376    //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)] 
    74     double sld_flat10; 
     77    Parameter sld_flat10; 
    7578 
    7679    //  [DEFAULT]=thick_inter1=50.0 [A] 
    77     double thick_inter1; 
     80    Parameter thick_inter1; 
    7881    //  [DEFAULT]=thick_inter2=50.0 [A] 
    79     double thick_inter2; 
     82    Parameter thick_inter2; 
    8083    //  [DEFAULT]=thick_inter3=50.0 [A] 
    81     double thick_inter3; 
     84    Parameter thick_inter3; 
    8285    //  [DEFAULT]=thick_inter4=50.0 [A] 
    83     double thick_inter4; 
     86    Parameter thick_inter4; 
    8487    //  [DEFAULT]=thick_inter5=50.0 [A] 
    85     double thick_inter5; 
     88    Parameter thick_inter5; 
    8689    //  [DEFAULT]=thick_inter6=50.0 [A] 
    87     double thick_inter6; 
     90    Parameter thick_inter6; 
    8891    //  [DEFAULT]=thick_inter7=50.0 [A] 
    89     double thick_inter7; 
     92    Parameter thick_inter7; 
    9093    //  [DEFAULT]=thick_inter8=50.0 [A] 
    91     double thick_inter8; 
     94    Parameter thick_inter8; 
    9295    //  [DEFAULT]=thick_inter9=50.0 [A] 
    93     double thick_inter9; 
     96    Parameter thick_inter9; 
    9497    //  [DEFAULT]=thick_inter10=50.0 [A] 
    95     double thick_inter10; 
     98    Parameter thick_inter10; 
    9699 
    97100    //  [DEFAULT]=thick_flat1=100 [A] 
    98     double thick_flat1; 
     101    Parameter thick_flat1; 
    99102    //  [DEFAULT]=thick_flat2=100 [A] 
    100     double thick_flat2; 
     103    Parameter thick_flat2; 
    101104    //  [DEFAULT]=thick_flat3=100 [A] 
    102     double thick_flat3; 
     105    Parameter thick_flat3; 
    103106    //  [DEFAULT]=thick_flat4=100 [A] 
    104     double thick_flat4; 
     107    Parameter thick_flat4; 
    105108    //  [DEFAULT]=thick_flat5=100 [A] 
    106     double thick_flat5; 
     109    Parameter thick_flat5; 
    107110    //  [DEFAULT]=thick_flat6=100 [A] 
    108     double thick_flat6; 
     111    Parameter thick_flat6; 
    109112    //  [DEFAULT]=thick_flat7=100 [A] 
    110     double thick_flat7; 
     113    Parameter thick_flat7; 
    111114    //  [DEFAULT]=thick_flat8=100 [A] 
    112     double thick_flat8; 
     115    Parameter thick_flat8; 
    113116    //  [DEFAULT]=thick_flat9=100 [A] 
    114     double thick_flat9; 
     117    Parameter thick_flat9; 
    115118    //  [DEFAULT]=thick_flat10=100 [A] 
    116     double thick_flat10; 
     119    Parameter thick_flat10; 
    117120 
    118121    //  [DEFAULT]=func_inter1=0 
    119     double func_inter1; 
     122    Parameter func_inter1; 
    120123    //  [DEFAULT]=func_inter2=0 
    121     double func_inter2; 
     124    Parameter func_inter2; 
    122125    //  [DEFAULT]=func_inter3=0 
    123     double func_inter3; 
     126    Parameter func_inter3; 
    124127    //  [DEFAULT]=func_inter4=0 
    125     double func_inter4; 
     128    Parameter func_inter4; 
    126129    //  [DEFAULT]=func_inter5=0 
    127     double func_inter5; 
     130    Parameter func_inter5; 
    128131    //  [DEFAULT]=func_inter6=0 
    129     double func_inter6; 
     132    Parameter func_inter6; 
    130133    //  [DEFAULT]=func_inter7=0 
    131     double func_inter7; 
     134    Parameter func_inter7; 
    132135    //  [DEFAULT]=func_inter8=0 
    133     double func_inter8; 
     136    Parameter func_inter8; 
    134137    //  [DEFAULT]=func_inter9=0 
    135     double func_inter9; 
     138    Parameter func_inter9; 
    136139    //  [DEFAULT]=func_inter10=0 
    137     double func_inter10; 
     140    Parameter func_inter10; 
    138141 
    139142    //  [DEFAULT]=nu_inter1=2.5 
    140     double nu_inter1; 
     143    Parameter nu_inter1; 
    141144    //  [DEFAULT]=nu_inter2=2.5 
    142     double nu_inter2; 
     145    Parameter nu_inter2; 
    143146    //  [DEFAULT]=nu_inter3=2.5 
    144     double nu_inter3; 
     147    Parameter nu_inter3; 
    145148    //  [DEFAULT]=nu_inter4=2.5 
    146     double nu_inter4; 
     149    Parameter nu_inter4; 
    147150    //  [DEFAULT]=nu_inter5=2.5 
    148     double nu_inter5; 
     151    Parameter nu_inter5; 
    149152    //  [DEFAULT]=nu_inter6=2.5 
    150     double nu_inter6; 
     153    Parameter nu_inter6; 
    151154    //  [DEFAULT]=nu_inter7=2.5 
    152     double nu_inter7; 
     155    Parameter nu_inter7; 
    153156    //  [DEFAULT]=nu_inter8=2.5 
    154     double nu_inter8; 
     157    Parameter nu_inter8; 
    155158    //  [DEFAULT]=nu_inter9=2.5 
    156     double nu_inter9; 
     159    Parameter nu_inter9; 
    157160    //  [DEFAULT]=nu_inter10=2.5 
    158     double nu_inter10; 
     161    Parameter nu_inter10; 
    159162 
    160163    //  [DEFAULT]=npts_inter=35.0 
    161     double npts_inter; 
     164    Parameter npts_inter; 
    162165    //  [DEFAULT]=nu_inter0=2.5 
    163     double nu_inter0; 
     166    Parameter nu_inter0; 
    164167    //  [DEFAULT]=rad_core0=50.0 [A] 
    165     double rad_core0; 
    166 } SphereSLDParameters; 
     168    Parameter rad_core0; 
    167169 
    168 double sphere_sld_kernel(double dq[], double q); 
     170  // Constructor 
     171  SphereSLDModel(); 
    169172 
    170 /// 1D scattering function 
    171 double sphere_sld_analytical_1D(SphereSLDParameters *pars, double q); 
     173  // Operators to get I(Q) 
     174  double operator()(double q); 
     175  double operator()(double qx, double qy); 
     176  double calculate_ER(); 
     177  double evaluate_rphi(double q, double phi); 
     178}; 
    172179 
    173 /// 2D scattering function 
    174 double sphere_sld_analytical_2D(SphereSLDParameters *pars, double q, double phi); 
    175 double sphere_sld_analytical_2DXY(SphereSLDParameters *pars, double qx, double qy); 
    176180 
    177181#endif 

  • sansmodels/src/c_models/models.hh

    ra8eab1c r37805e9  
    3030 
    3131 
    32  
    33 class SphereSLDModel{ 
    34 public: 
    35         // Model parameters 
    36         Parameter n_shells; 
    37         Parameter scale; 
    38         Parameter thick_inter0; 
    39         Parameter func_inter0; 
    40         Parameter sld_core0; 
    41         Parameter sld_solv; 
    42         Parameter background; 
    43  
    44         Parameter sld_flat1; 
    45         Parameter sld_flat2; 
    46         Parameter sld_flat3; 
    47         Parameter sld_flat4; 
    48         Parameter sld_flat5; 
    49         Parameter sld_flat6; 
    50         Parameter sld_flat7; 
    51         Parameter sld_flat8; 
    52         Parameter sld_flat9; 
    53         Parameter sld_flat10; 
    54  
    55         Parameter thick_inter1; 
    56         Parameter thick_inter2; 
    57         Parameter thick_inter3; 
    58         Parameter thick_inter4; 
    59         Parameter thick_inter5; 
    60         Parameter thick_inter6; 
    61         Parameter thick_inter7; 
    62         Parameter thick_inter8; 
    63         Parameter thick_inter9; 
    64         Parameter thick_inter10; 
    65  
    66         Parameter thick_flat1; 
    67         Parameter thick_flat2; 
    68         Parameter thick_flat3; 
    69         Parameter thick_flat4; 
    70         Parameter thick_flat5; 
    71         Parameter thick_flat6; 
    72         Parameter thick_flat7; 
    73         Parameter thick_flat8; 
    74         Parameter thick_flat9; 
    75         Parameter thick_flat10; 
    76  
    77         Parameter func_inter1; 
    78         Parameter func_inter2; 
    79         Parameter func_inter3; 
    80         Parameter func_inter4; 
    81         Parameter func_inter5; 
    82         Parameter func_inter6; 
    83         Parameter func_inter7; 
    84         Parameter func_inter8; 
    85         Parameter func_inter9; 
    86         Parameter func_inter10; 
    87  
    88         Parameter nu_inter1; 
    89         Parameter nu_inter2; 
    90         Parameter nu_inter3; 
    91         Parameter nu_inter4; 
    92         Parameter nu_inter5; 
    93         Parameter nu_inter6; 
    94         Parameter nu_inter7; 
    95         Parameter nu_inter8; 
    96         Parameter nu_inter9; 
    97         Parameter nu_inter10; 
    98  
    99         Parameter npts_inter; 
    100         Parameter nu_inter0; 
    101         Parameter rad_core0; 
    102  
    103         // Constructor 
    104         SphereSLDModel(); 
    105  
    106         // Operators to get I(Q) 
    107         double operator()(double q); 
    108         double operator()(double qx, double qy); 
    109         double calculate_ER(); 
    110         double evaluate_rphi(double q, double phi); 
    111 }; 
    11232 
    11333 

  • sansmodels/src/c_models/refl.cpp

    r2d1b700 r37805e9  
    1313#define lamda 4.62 
    1414 
    15 double re_kernel(double dp[], double q) { 
     15static double re_kernel(double dp[], double q) { 
    1616  int n = dp[0]; 
    1717  int i,j; 

  • sansmodels/src/c_models/spheresld.cpp

    r67424cd r37805e9  
    11 
    22#include <math.h> 
    3 #include "models.hh" 
    43#include "parameters.hh" 
    54#include <stdio.h> 
     5#include <stdlib.h> 
     6 
     7extern "C" { 
     8#include "libmultifunc/librefl.h" 
     9} 
     10 
     11#include "spheresld.h" 
     12 
    613using namespace std; 
    714 
    8 extern "C" { 
    9         #include "spheresld.h" 
     15// Convenience structure 
     16typedef struct { 
     17  /// number of shells 
     18  //  [DEFAULT]=n_shells=1 
     19  int n_shells; 
     20  /// Scale factor 
     21  //  [DEFAULT]=scale= 1.0 
     22  double scale; 
     23  /// thick_inter0 [A] 
     24  //  [DEFAULT]=thick_inter0=50.0 [A] 
     25  double thick_inter0; 
     26  /// func_inter0 
     27  //  [DEFAULT]=func_inter0= 0 
     28  double func_inter0; 
     29  /// sld_core0 [1/A^(2)] 
     30  //  [DEFAULT]=sld_core0= 2.07e-6 [1/A^(2)] 
     31  double sld_core0; 
     32  /// sld_solv [1/A^(2)] 
     33  //  [DEFAULT]=sld_solv= 1.0e-6 [1/A^(2)] 
     34  double sld_solv; 
     35  /// Background 
     36  //  [DEFAULT]=background=0 
     37  double background; 
     38 
     39  //  [DEFAULT]=sld_flat1=4.0e-06 [1/A^(2)] 
     40  double sld_flat1; 
     41  //  [DEFAULT]=sld_flat2=3.5e-06 [1/A^(2)] 
     42  double sld_flat2; 
     43  //  [DEFAULT]=sld_flat3=4.0e-06 [1/A^(2)] 
     44  double sld_flat3; 
     45  //  [DEFAULT]=sld_flat4=3.5e-06 [1/A^(2)] 
     46  double sld_flat4; 
     47  //  [DEFAULT]=sld_flat5=4.0e-06 [1/A^(2)] 
     48  double sld_flat5; 
     49  //  [DEFAULT]=sld_flat6=3.5e-06 [1/A^(2)] 
     50  double sld_flat6; 
     51  //  [DEFAULT]=sld_flat7=4.0e-06 [1/A^(2)] 
     52  double sld_flat7; 
     53  //  [DEFAULT]=sld_flat8=3.5e-06 [1/A^(2)] 
     54  double sld_flat8; 
     55  //  [DEFAULT]=sld_flat9=4.0e-06 [1/A^(2)] 
     56  double sld_flat9; 
     57  //  [DEFAULT]=sld_flat10=3.5e-06 [1/A^(2)] 
     58  double sld_flat10; 
     59 
     60  //  [DEFAULT]=thick_inter1=50.0 [A] 
     61  double thick_inter1; 
     62  //  [DEFAULT]=thick_inter2=50.0 [A] 
     63  double thick_inter2; 
     64  //  [DEFAULT]=thick_inter3=50.0 [A] 
     65  double thick_inter3; 
     66  //  [DEFAULT]=thick_inter4=50.0 [A] 
     67  double thick_inter4; 
     68  //  [DEFAULT]=thick_inter5=50.0 [A] 
     69  double thick_inter5; 
     70  //  [DEFAULT]=thick_inter6=50.0 [A] 
     71  double thick_inter6; 
     72  //  [DEFAULT]=thick_inter7=50.0 [A] 
     73  double thick_inter7; 
     74  //  [DEFAULT]=thick_inter8=50.0 [A] 
     75  double thick_inter8; 
     76  //  [DEFAULT]=thick_inter9=50.0 [A] 
     77  double thick_inter9; 
     78  //  [DEFAULT]=thick_inter10=50.0 [A] 
     79  double thick_inter10; 
     80 
     81  //  [DEFAULT]=thick_flat1=100 [A] 
     82  double thick_flat1; 
     83  //  [DEFAULT]=thick_flat2=100 [A] 
     84  double thick_flat2; 
     85  //  [DEFAULT]=thick_flat3=100 [A] 
     86  double thick_flat3; 
     87  //  [DEFAULT]=thick_flat4=100 [A] 
     88  double thick_flat4; 
     89  //  [DEFAULT]=thick_flat5=100 [A] 
     90  double thick_flat5; 
     91  //  [DEFAULT]=thick_flat6=100 [A] 
     92  double thick_flat6; 
     93  //  [DEFAULT]=thick_flat7=100 [A] 
     94  double thick_flat7; 
     95  //  [DEFAULT]=thick_flat8=100 [A] 
     96  double thick_flat8; 
     97  //  [DEFAULT]=thick_flat9=100 [A] 
     98  double thick_flat9; 
     99  //  [DEFAULT]=thick_flat10=100 [A] 
     100  double thick_flat10; 
     101 
     102  //  [DEFAULT]=func_inter1=0 
     103  double func_inter1; 
     104  //  [DEFAULT]=func_inter2=0 
     105  double func_inter2; 
     106  //  [DEFAULT]=func_inter3=0 
     107  double func_inter3; 
     108  //  [DEFAULT]=func_inter4=0 
     109  double func_inter4; 
     110  //  [DEFAULT]=func_inter5=0 
     111  double func_inter5; 
     112  //  [DEFAULT]=func_inter6=0 
     113  double func_inter6; 
     114  //  [DEFAULT]=func_inter7=0 
     115  double func_inter7; 
     116  //  [DEFAULT]=func_inter8=0 
     117  double func_inter8; 
     118  //  [DEFAULT]=func_inter9=0 
     119  double func_inter9; 
     120  //  [DEFAULT]=func_inter10=0 
     121  double func_inter10; 
     122 
     123  //  [DEFAULT]=nu_inter1=2.5 
     124  double nu_inter1; 
     125  //  [DEFAULT]=nu_inter2=2.5 
     126  double nu_inter2; 
     127  //  [DEFAULT]=nu_inter3=2.5 
     128  double nu_inter3; 
     129  //  [DEFAULT]=nu_inter4=2.5 
     130  double nu_inter4; 
     131  //  [DEFAULT]=nu_inter5=2.5 
     132  double nu_inter5; 
     133  //  [DEFAULT]=nu_inter6=2.5 
     134  double nu_inter6; 
     135  //  [DEFAULT]=nu_inter7=2.5 
     136  double nu_inter7; 
     137  //  [DEFAULT]=nu_inter8=2.5 
     138  double nu_inter8; 
     139  //  [DEFAULT]=nu_inter9=2.5 
     140  double nu_inter9; 
     141  //  [DEFAULT]=nu_inter10=2.5 
     142  double nu_inter10; 
     143 
     144  //  [DEFAULT]=npts_inter=35.0 
     145  double npts_inter; 
     146  //  [DEFAULT]=nu_inter0=2.5 
     147  double nu_inter0; 
     148  //  [DEFAULT]=rad_core0=50.0 [A] 
     149  double rad_core0; 
     150} SphereSLDParameters; 
     151 
     152 
     153static double sphere_sld_kernel(double dp[], double q) { 
     154  int n = dp[0]; 
     155  int i,j,k; 
     156 
     157  double scale = dp[1]; 
     158  double thick_inter_core = dp[2]; 
     159  double sld_core = dp[4]; 
     160  double sld_solv = dp[5]; 
     161  double background = dp[6]; 
     162  double npts = dp[57]; //number of sub_layers in each interface 
     163  double nsl=npts;//21.0; //nsl = Num_sub_layer:  MUST ODD number in double //no other number works now 
     164  int n_s; 
     165 
     166  double sld_i,sld_f,dz,bes,fun,f,vol,vol_pre,vol_sub,qr,r,contr,f2; 
     167  double sign,slope=0.0; 
     168  double pi; 
     169 
     170  int* fun_type; 
     171  double* sld; 
     172  double* thick_inter; 
     173  double* thick; 
     174  double* fun_coef; 
     175 
     176  double total_thick=0.0; 
     177 
     178  fun_type = (int*)malloc((n+2)*sizeof(int)); 
     179  sld = (double*)malloc((n+2)*sizeof(double)); 
     180  thick_inter = (double*)malloc((n+2)*sizeof(double)); 
     181  thick = (double*)malloc((n+2)*sizeof(double)); 
     182  fun_coef = (double*)malloc((n+2)*sizeof(double)); 
     183 
     184  fun_type[0] = dp[3]; 
     185  fun_coef[0] = fabs(dp[58]); 
     186  for (i =1; i<=n; i++){ 
     187    sld[i] = dp[i+6]; 
     188    thick_inter[i]= dp[i+16]; 
     189    thick[i] = dp[i+26]; 
     190    fun_type[i] = dp[i+36]; 
     191    fun_coef[i] = fabs(dp[i+46]); 
     192    total_thick += thick[i] + thick_inter[i]; 
     193  } 
     194  sld[0] = sld_core; 
     195  sld[n+1] = sld_solv; 
     196  thick[0] = dp[59]; 
     197  thick[n+1] = total_thick/5.0; 
     198  thick_inter[0] = thick_inter_core; 
     199  thick_inter[n+1] = 0.0; 
     200  fun_coef[n+1] = 0.0; 
     201 
     202  pi = 4.0*atan(1.0); 
     203  f = 0.0; 
     204  r = 0.0; 
     205  vol = 0.0; 
     206  vol_pre = 0.0; 
     207  vol_sub = 0.0; 
     208  sld_f = sld_core; 
     209 
     210  //floor_nsl = floor(nsl/2.0); 
     211 
     212  dz = 0.0; 
     213  // iteration for # of shells + core + solvent 
     214  for (i=0;i<=n+1; i++){ 
     215    //iteration for N sub-layers 
     216    //if (fabs(thick[i]) <= 1e-24){ 
     217    //   continue; 
     218    //} 
     219    // iteration for flat and interface 
     220    for (j=0;j<2;j++){ 
     221      // iteration for sub_shells in the interface 
     222      // starts from #1 sub-layer 
     223      for (n_s=1;n_s<=nsl; n_s++){ 
     224        // for solvent, it doesn't have an interface 
     225        if (i==n+1 && j==1) 
     226          break; 
     227        // for flat layers 
     228        if (j==0){ 
     229          dz = thick[i]; 
     230          sld_i = sld[i]; 
     231          slope = 0.0; 
     232        } 
     233        // for interfacial sub_shells 
     234        else{ 
     235          dz = thick_inter[i]/nsl; 
     236          // find sld_i at the outer boundary of sub-layer #n_s 
     237          sld_i = intersldfunc(fun_type[i],nsl, n_s, fun_coef[i], sld[i], sld[i+1]); 
     238          // calculate slope 
     239          slope= (sld_i -sld_f)/dz; 
     240        } 
     241        contr = sld_f-slope*r; 
     242        // iteration for the left and right boundary of the shells(or sub_shells) 
     243        for (k=0; k<2; k++){ 
     244          // At r=0, the contribution to I is zero, so skip it. 
     245          if ( i == 0 && j == 0 && k == 0){ 
     246            continue; 
     247          } 
     248          // On the top of slovent there is no interface; skip it. 
     249          if (i == n+1 && k == 1){ 
     250            continue; 
     251          } 
     252          // At the right side (outer) boundary 
     253          if ( k == 1){ 
     254            sign = 1.0; 
     255            r += dz; 
     256          } 
     257          // At the left side (inner) boundary 
     258          else{ 
     259            sign = -1.0; 
     260          } 
     261          qr = q * r; 
     262          fun = 0.0; 
     263          if(qr == 0.0){ 
     264            // sigular point 
     265            bes = sign * 1.0; 
     266          } 
     267          else{ 
     268            // for flat sub-layer 
     269            bes = sign *  3.0 * (sin(qr) - qr * cos(qr)) / (qr * qr * qr); 
     270            // with linear slope 
     271            if (fabs(slope) > 0.0 ){ 
     272              fun = sign * 3.0 * r * (2.0*qr*sin(qr)-((qr*qr)-2.0)*cos(qr))/(qr * qr * qr * qr); 
     273            } 
     274          } 
     275          // update total volume 
     276          vol = 4.0 * pi / 3.0 * r * r * r; 
     277          // we won't do the following volume correction for now. 
     278          // substrate empty area of volume 
     279          //if (k == 1 && fabs(sld_in[i]-sld_solv) < 1e-04*fabs(sld_solv) && fun_type[i]==0){ 
     280          //  vol_sub += (vol_pre - vol); 
     281          //} 
     282          f += vol * (bes * contr + fun * slope); 
     283        } 
     284        // remember this sld as sld_f 
     285        sld_f =sld_i; 
     286        // no sub-layer iteration (n_s loop) for the flat layer 
     287        if (j==0) 
     288          break; 
     289      } 
     290    } 
     291  } 
     292  //vol += vol_sub; 
     293  f2 = f * f / vol * 1.0e8; 
     294  f2 *= scale; 
     295  f2 += background; 
     296 
     297  free(fun_type); 
     298  free(sld); 
     299  free(thick_inter); 
     300  free(thick); 
     301  free(fun_coef); 
     302 
     303  return (f2); 
    10304} 
    11305 
     306 
    12307SphereSLDModel :: SphereSLDModel() { 
    13         n_shells = Parameter(1.0); 
    14         scale = Parameter(1.0); 
    15         thick_inter0 = Parameter(1.0, true); 
    16         thick_inter0.set_min(0.0); 
    17         func_inter0 = Parameter(0); 
    18         sld_core0 = Parameter(2.07e-06); 
    19         sld_solv = Parameter(1.0e-06); 
    20     background = Parameter(0.0); 
    21  
    22  
    23     sld_flat1 = Parameter(2.7e-06); 
    24     sld_flat2 = Parameter(3.5e-06); 
    25     sld_flat3 = Parameter(4.0e-06); 
    26     sld_flat4 = Parameter(3.5e-06); 
    27     sld_flat5 = Parameter(4.0e-06); 
    28     sld_flat6 = Parameter(3.5e-06); 
    29     sld_flat7 = Parameter(4.0e-06); 
    30     sld_flat8 = Parameter(3.5e-06); 
    31     sld_flat9 = Parameter(4.0e-06); 
    32     sld_flat10 = Parameter(3.5e-06); 
    33  
    34  
    35     thick_inter1 = Parameter(1.0); 
    36     thick_inter2 = Parameter(1.0); 
    37     thick_inter3 = Parameter(1.0); 
    38     thick_inter4 = Parameter(1.0); 
    39     thick_inter5 = Parameter(1.0); 
    40     thick_inter6 = Parameter(1.0); 
    41     thick_inter7 = Parameter(1.0); 
    42     thick_inter8 = Parameter(1.0); 
    43     thick_inter9 = Parameter(1.0); 
    44     thick_inter10 = Parameter(1.0); 
    45  
    46  
    47     thick_flat1 = Parameter(100.0); 
    48     thick_flat2 = Parameter(100.0); 
    49     thick_flat3 = Parameter(100.0); 
    50     thick_flat4 = Parameter(100.0); 
    51     thick_flat5 = Parameter(100.0); 
    52     thick_flat6 = Parameter(100.0); 
    53     thick_flat7 = Parameter(100.0); 
    54     thick_flat8 = Parameter(100.0); 
    55     thick_flat9 = Parameter(100.0); 
    56     thick_flat10 = Parameter(100.0); 
    57  
    58  
    59     func_inter1 = Parameter(0); 
    60     func_inter2 = Parameter(0); 
    61     func_inter3 = Parameter(0); 
    62     func_inter4 = Parameter(0); 
    63     func_inter5 = Parameter(0); 
    64     func_inter6 = Parameter(0); 
    65     func_inter7 = Parameter(0); 
    66     func_inter8 = Parameter(0); 
    67     func_inter9 = Parameter(0); 
    68     func_inter10 = Parameter(0); 
    69  
    70     nu_inter1 = Parameter(2.5); 
    71     nu_inter2 = Parameter(2.5); 
    72     nu_inter3 = Parameter(2.5); 
    73     nu_inter4 = Parameter(2.5); 
    74     nu_inter5 = Parameter(2.5); 
    75     nu_inter6 = Parameter(2.5); 
    76     nu_inter7 = Parameter(2.5); 
    77     nu_inter8 = Parameter(2.5); 
    78     nu_inter9 = Parameter(2.5); 
    79     nu_inter10 = Parameter(2.5); 
    80  
    81     npts_inter = Parameter(35.0); 
    82     nu_inter0 = Parameter(2.5); 
    83     rad_core0 = Parameter(60.0, true); 
    84     rad_core0.set_min(0.0); 
     308  n_shells = Parameter(1.0); 
     309  scale = Parameter(1.0); 
     310  thick_inter0 = Parameter(1.0, true); 
     311  thick_inter0.set_min(0.0); 
     312  func_inter0 = Parameter(0); 
     313  sld_core0 = Parameter(2.07e-06); 
     314  sld_solv = Parameter(1.0e-06); 
     315  background = Parameter(0.0); 
     316 
     317 
     318  sld_flat1 = Parameter(2.7e-06); 
     319  sld_flat2 = Parameter(3.5e-06); 
     320  sld_flat3 = Parameter(4.0e-06); 
     321  sld_flat4 = Parameter(3.5e-06); 
     322  sld_flat5 = Parameter(4.0e-06); 
     323  sld_flat6 = Parameter(3.5e-06); 
     324  sld_flat7 = Parameter(4.0e-06); 
     325  sld_flat8 = Parameter(3.5e-06); 
     326  sld_flat9 = Parameter(4.0e-06); 
     327  sld_flat10 = Parameter(3.5e-06); 
     328 
     329 
     330  thick_inter1 = Parameter(1.0); 
     331  thick_inter2 = Parameter(1.0); 
     332  thick_inter3 = Parameter(1.0); 
     333  thick_inter4 = Parameter(1.0); 
     334  thick_inter5 = Parameter(1.0); 
     335  thick_inter6 = Parameter(1.0); 
     336  thick_inter7 = Parameter(1.0); 
     337  thick_inter8 = Parameter(1.0); 
     338  thick_inter9 = Parameter(1.0); 
     339  thick_inter10 = Parameter(1.0); 
     340 
     341 
     342  thick_flat1 = Parameter(100.0); 
     343  thick_flat2 = Parameter(100.0); 
     344  thick_flat3 = Parameter(100.0); 
     345  thick_flat4 = Parameter(100.0); 
     346  thick_flat5 = Parameter(100.0); 
     347  thick_flat6 = Parameter(100.0); 
     348  thick_flat7 = Parameter(100.0); 
     349  thick_flat8 = Parameter(100.0); 
     350  thick_flat9 = Parameter(100.0); 
     351  thick_flat10 = Parameter(100.0); 
     352 
     353 
     354  func_inter1 = Parameter(0); 
     355  func_inter2 = Parameter(0); 
     356  func_inter3 = Parameter(0); 
     357  func_inter4 = Parameter(0); 
     358  func_inter5 = Parameter(0); 
     359  func_inter6 = Parameter(0); 
     360  func_inter7 = Parameter(0); 
     361  func_inter8 = Parameter(0); 
     362  func_inter9 = Parameter(0); 
     363  func_inter10 = Parameter(0); 
     364 
     365  nu_inter1 = Parameter(2.5); 
     366  nu_inter2 = Parameter(2.5); 
     367  nu_inter3 = Parameter(2.5); 
     368  nu_inter4 = Parameter(2.5); 
     369  nu_inter5 = Parameter(2.5); 
     370  nu_inter6 = Parameter(2.5); 
     371  nu_inter7 = Parameter(2.5); 
     372  nu_inter8 = Parameter(2.5); 
     373  nu_inter9 = Parameter(2.5); 
     374  nu_inter10 = Parameter(2.5); 
     375 
     376  npts_inter = Parameter(35.0); 
     377  nu_inter0 = Parameter(2.5); 
     378  rad_core0 = Parameter(60.0, true); 
     379  rad_core0.set_min(0.0); 
    85380} 
    86381 
    … …  
    91386 */ 
    92387double SphereSLDModel :: operator()(double q) { 
    93         double dp[60]; 
    94         // Fill parameter array for IGOR library 
    95         // Add the background after averaging 
    96         dp[0] = n_shells(); 
    97         dp[1] = scale(); 
    98         dp[2] = thick_inter0(); 
    99         dp[3] = func_inter0(); 
    100         dp[4] = sld_core0(); 
    101         dp[5] = sld_solv(); 
    102         dp[6] = 0.0; 
    103  
    104         dp[7] = sld_flat1(); 
    105         dp[8] = sld_flat2(); 
    106         dp[9] = sld_flat3(); 
    107         dp[10] = sld_flat4(); 
    108         dp[11] = sld_flat5(); 
    109         dp[12] = sld_flat6(); 
    110         dp[13] = sld_flat7(); 
    111         dp[14] = sld_flat8(); 
    112         dp[15] = sld_flat9(); 
    113         dp[16] = sld_flat10(); 
    114  
    115         dp[17] = thick_inter1(); 
    116         dp[18] = thick_inter2(); 
    117         dp[19] = thick_inter3(); 
    118         dp[20] = thick_inter4(); 
    119         dp[21] = thick_inter5(); 
    120         dp[22] = thick_inter6(); 
    121         dp[23] = thick_inter7(); 
    122         dp[24] = thick_inter8(); 
    123         dp[25] = thick_inter9(); 
    124         dp[26] = thick_inter10(); 
    125  
    126         dp[27] = thick_flat1(); 
    127         dp[28] = thick_flat2(); 
    128         dp[29] = thick_flat3(); 
    129         dp[30] = thick_flat4(); 
    130         dp[31] = thick_flat5(); 
    131         dp[32] = thick_flat6(); 
    132         dp[33] = thick_flat7(); 
    133         dp[34] = thick_flat8(); 
    134         dp[35] = thick_flat9(); 
    135         dp[36] = thick_flat10(); 
    136  
    137         dp[37] = func_inter1(); 
    138         dp[38] = func_inter2(); 
    139         dp[39] = func_inter3(); 
    140         dp[40] = func_inter4(); 
    141         dp[41] = func_inter5(); 
    142         dp[42] = func_inter6(); 
    143         dp[43] = func_inter7(); 
    144         dp[44] = func_inter8(); 
    145         dp[45] = func_inter9(); 
    146         dp[46] = func_inter10(); 
    147  
    148         dp[47] = nu_inter1(); 
    149         dp[48] = nu_inter2(); 
    150         dp[49] = nu_inter3(); 
    151         dp[50] = nu_inter4(); 
    152         dp[51] = nu_inter5(); 
    153         dp[52] = nu_inter6(); 
    154         dp[53] = nu_inter7(); 
    155         dp[54] = nu_inter8(); 
    156         dp[55] = nu_inter9(); 
    157         dp[56] = nu_inter10(); 
    158  
    159  
    160         dp[57] = npts_inter(); 
    161         dp[58] = nu_inter0(); 
    162         dp[59] = rad_core0(); 
    163  
    164         // No polydispersion supported in this model. 
    165         // Get the dispersion points for the radius 
    166         vector<WeightPoint> weights_rad_core0; 
    167         rad_core0.get_weights(weights_rad_core0); 
    168         vector<WeightPoint> weights_thick_inter0; 
    169         thick_inter0.get_weights(weights_thick_inter0); 
    170         // Perform the computation, with all weight points 
    171         double sum = 0.0; 
    172         double norm = 0.0; 
    173         double vol = 0.0; 
    174  
    175         // Loop over core weight points 
    176         for(size_t i=0; i<weights_rad_core0.size(); i++) { 
    177                 dp[59] = weights_rad_core0[i].value; 
    178                 // Loop over thick_inter0 weight points 
    179                 for(size_t j=0; j<weights_thick_inter0.size(); j++) { 
    180                         dp[2] = weights_thick_inter0[j].value; 
    181  
    182                         //Un-normalize Sphere by volume 
    183                         sum += weights_rad_core0[i].weight * weights_thick_inter0[j].weight 
    184                                 * sphere_sld_kernel(dp,q) * pow((weights_rad_core0[i].value + 
    185                                                 weights_thick_inter0[j].value),3.0); 
    186                         //Find average volume 
    187                         vol += weights_rad_core0[i].weight * weights_thick_inter0[j].weight 
    188                                 * pow((weights_rad_core0[i].value+weights_thick_inter0[j].value),3.0); 
    189  
    190                         norm += weights_rad_core0[i].weight * weights_thick_inter0[j].weight; 
    191                 } 
    192         } 
    193  
    194         if (vol != 0.0 && norm != 0.0) { 
    195                 //Re-normalize by avg volume 
    196                 sum = sum/(vol/norm);} 
    197  
    198         return sum/norm + background(); 
     388  double dp[60]; 
     389  // Fill parameter array for IGOR library 
     390  // Add the background after averaging 
     391  dp[0] = n_shells(); 
     392  dp[1] = scale(); 
     393  dp[2] = thick_inter0(); 
     394  dp[3] = func_inter0(); 
     395  dp[4] = sld_core0(); 
     396  dp[5] = sld_solv(); 
     397  dp[6] = 0.0; 
     398 
     399  dp[7] = sld_flat1(); 
     400  dp[8] = sld_flat2(); 
     401  dp[9] = sld_flat3(); 
     402  dp[10] = sld_flat4(); 
     403  dp[11] = sld_flat5(); 
     404  dp[12] = sld_flat6(); 
     405  dp[13] = sld_flat7(); 
     406  dp[14] = sld_flat8(); 
     407  dp[15] = sld_flat9(); 
     408  dp[16] = sld_flat10(); 
     409 
     410  dp[17] = thick_inter1(); 
     411  dp[18] = thick_inter2(); 
     412  dp[19] = thick_inter3(); 
     413  dp[20] = thick_inter4(); 
     414  dp[21] = thick_inter5(); 
     415  dp[22] = thick_inter6(); 
     416  dp[23] = thick_inter7(); 
     417  dp[24] = thick_inter8(); 
     418  dp[25] = thick_inter9(); 
     419  dp[26] = thick_inter10(); 
     420 
     421  dp[27] = thick_flat1(); 
     422  dp[28] = thick_flat2(); 
     423  dp[29] = thick_flat3(); 
     424  dp[30] = thick_flat4(); 
     425  dp[31] = thick_flat5(); 
     426  dp[32] = thick_flat6(); 
     427  dp[33] = thick_flat7(); 
     428  dp[34] = thick_flat8(); 
     429  dp[35] = thick_flat9(); 
     430  dp[36] = thick_flat10(); 
     431 
     432  dp[37] = func_inter1(); 
     433  dp[38] = func_inter2(); 
     434  dp[39] = func_inter3(); 
     435  dp[40] = func_inter4(); 
     436  dp[41] = func_inter5(); 
     437  dp[42] = func_inter6(); 
     438  dp[43] = func_inter7(); 
     439  dp[44] = func_inter8(); 
     440  dp[45] = func_inter9(); 
     441  dp[46] = func_inter10(); 
     442 
     443  dp[47] = nu_inter1(); 
     444  dp[48] = nu_inter2(); 
     445  dp[49] = nu_inter3(); 
     446  dp[50] = nu_inter4(); 
     447  dp[51] = nu_inter5(); 
     448  dp[52] = nu_inter6(); 
     449  dp[53] = nu_inter7(); 
     450  dp[54] = nu_inter8(); 
     451  dp[55] = nu_inter9(); 
     452  dp[56] = nu_inter10(); 
     453 
     454 
     455  dp[57] = npts_inter(); 
     456  dp[58] = nu_inter0(); 
     457  dp[59] = rad_core0(); 
     458 
     459  // No polydispersion supported in this model. 
     460  // Get the dispersion points for the radius 
     461  vector<WeightPoint> weights_rad_core0; 
     462  rad_core0.get_weights(weights_rad_core0); 
     463  vector<WeightPoint> weights_thick_inter0; 
     464  thick_inter0.get_weights(weights_thick_inter0); 
     465  // Perform the computation, with all weight points 
     466  double sum = 0.0; 
     467  double norm = 0.0; 
     468  double vol = 0.0; 
     469 
     470  // Loop over core weight points 
     471  for(size_t i=0; i<weights_rad_core0.size(); i++) { 
     472    dp[59] = weights_rad_core0[i].value; 
     473    // Loop over thick_inter0 weight points 
     474    for(size_t j=0; j<weights_thick_inter0.size(); j++) { 
     475      dp[2] = weights_thick_inter0[j].value; 
     476 
     477      //Un-normalize Sphere by volume 
     478      sum += weights_rad_core0[i].weight * weights_thick_inter0[j].weight 
     479          * sphere_sld_kernel(dp,q) * pow((weights_rad_core0[i].value + 
     480              weights_thick_inter0[j].value),3.0); 
     481      //Find average volume 
     482      vol += weights_rad_core0[i].weight * weights_thick_inter0[j].weight 
     483          * pow((weights_rad_core0[i].value+weights_thick_inter0[j].value),3.0); 
     484 
     485      norm += weights_rad_core0[i].weight * weights_thick_inter0[j].weight; 
     486    } 
     487  } 
     488 
     489  if (vol != 0.0 && norm != 0.0) { 
     490    //Re-normalize by avg volume 
     491    sum = sum/(vol/norm);} 
     492 
     493  return sum/norm + background(); 
    199494} 
    200495 
    … …  
    206501 */ 
    207502double SphereSLDModel :: operator()(double qx, double qy) { 
    208         double q = sqrt(qx*qx + qy*qy); 
    209         return (*this).operator()(q); 
     503  double q = sqrt(qx*qx + qy*qy); 
     504  return (*this).operator()(q); 
    210505} 
    211506 
    … …  
    218513 */ 
    219514double SphereSLDModel :: evaluate_rphi(double q, double phi) { 
    220         return (*this).operator()(q); 
     515  return (*this).operator()(q); 
    221516} 
    222517 
    … …  
    231526// sld of solvent. 
    232527double SphereSLDModel :: calculate_ER() { 
    233         SphereSLDParameters dp; 
    234  
    235         dp.n_shells = n_shells(); 
    236  
    237         dp.rad_core0 = rad_core0(); 
    238         dp.thick_flat1 = thick_flat1(); 
    239         dp.thick_flat2 = thick_flat2(); 
    240         dp.thick_flat3 = thick_flat3(); 
    241         dp.thick_flat4 = thick_flat4(); 
    242         dp.thick_flat5 = thick_flat5(); 
    243         dp.thick_flat6 = thick_flat6(); 
    244         dp.thick_flat7 = thick_flat7(); 
    245         dp.thick_flat8 = thick_flat8(); 
    246         dp.thick_flat9 = thick_flat9(); 
    247         dp.thick_flat10 = thick_flat10(); 
    248  
    249         dp.thick_inter0 = thick_inter0(); 
    250         dp.thick_inter1 = thick_inter1(); 
    251         dp.thick_inter2 = thick_inter2(); 
    252         dp.thick_inter3 = thick_inter3(); 
    253         dp.thick_inter4 = thick_inter4(); 
    254         dp.thick_inter5 = thick_inter5(); 
    255         dp.thick_inter6 = thick_inter6(); 
    256         dp.thick_inter7 = thick_inter7(); 
    257         dp.thick_inter8 = thick_inter8(); 
    258         dp.thick_inter9 = thick_inter9(); 
    259         dp.thick_inter10 = thick_inter10(); 
    260  
    261         double rad_out = 0.0; 
    262         double out = 0.0; 
    263         // Perform the computation, with all weight points 
    264         double sum = 0.0; 
    265         double norm = 0.0; 
    266  
    267         // Get the dispersion points for the radius 
    268         vector<WeightPoint> weights_rad_core0; 
    269         rad_core0.get_weights(weights_rad_core0); 
    270  
    271         // Get the dispersion points for the thick 1 
    272         vector<WeightPoint> weights_thick_inter0; 
    273         thick_inter0.get_weights(weights_thick_inter0); 
    274         // Loop over radius weight points 
    275         for(size_t i=0; i<weights_rad_core0.size(); i++) { 
    276                 dp.rad_core0 = weights_rad_core0[i].value; 
    277                 // Loop over radius weight points 
    278                 for(size_t j=0; j<weights_thick_inter0.size(); j++) { 
    279                         dp.thick_inter0 = weights_thick_inter0[j].value; 
    280                         rad_out = dp.rad_core0 + dp.thick_inter0; 
    281                         if (dp.n_shells > 0) 
    282                                 rad_out += dp.thick_flat1 + dp.thick_inter1; 
    283                         if (dp.n_shells > 1) 
    284                                 rad_out += dp.thick_flat2 + dp.thick_inter2; 
    285                         if (dp.n_shells > 2) 
    286                                 rad_out += dp.thick_flat3 + dp.thick_inter3; 
    287                         if (dp.n_shells > 3) 
    288                                 rad_out += dp.thick_flat4 + dp.thick_inter4; 
    289                         if (dp.n_shells > 4) 
    290                                 rad_out += dp.thick_flat5 + dp.thick_inter5; 
    291                         if (dp.n_shells > 5) 
    292                                 rad_out += dp.thick_flat6 + dp.thick_inter6; 
    293                         if (dp.n_shells > 6) 
    294                                 rad_out += dp.thick_flat7 + dp.thick_inter7; 
    295                         if (dp.n_shells > 7) 
    296                                 rad_out += dp.thick_flat8 + dp.thick_inter8; 
    297                         if (dp.n_shells > 8) 
    298                                 rad_out += dp.thick_flat9 + dp.thick_inter9; 
    299                         if (dp.n_shells > 9) 
    300                                 rad_out += dp.thick_flat10 + dp.thick_inter10; 
    301                         sum += weights_rad_core0[i].weight*weights_thick_inter0[j].weight 
    302                                 * (rad_out); 
    303                         norm += weights_rad_core0[i].weight*weights_thick_inter0[j].weight; 
    304                 } 
    305         } 
    306         if (norm != 0){ 
    307                 //return the averaged value 
    308                 out =  sum/norm;} 
    309         else{ 
    310                 //return normal value 
    311                 out = dp.rad_core0 + dp.thick_inter0; 
    312                 if (dp.n_shells > 0) 
    313                         out += dp.thick_flat1 + dp.thick_inter1; 
    314                 if (dp.n_shells > 1) 
    315                         out += dp.thick_flat2 + dp.thick_inter2; 
    316                 if (dp.n_shells > 2) 
    317                         out += dp.thick_flat3 + dp.thick_inter3; 
    318                 if (dp.n_shells > 3) 
    319                         out += dp.thick_flat4 + dp.thick_inter4; 
    320                 if (dp.n_shells > 4) 
    321                         out += dp.thick_flat5 + dp.thick_inter5; 
    322                 if (dp.n_shells > 5) 
    323                         out += dp.thick_flat6 + dp.thick_inter6; 
    324                 if (dp.n_shells > 6) 
    325                         out += dp.thick_flat7 + dp.thick_inter7; 
    326                 if (dp.n_shells > 7) 
    327                         out += dp.thick_flat8 + dp.thick_inter8; 
    328                 if (dp.n_shells > 8) 
    329                         out += dp.thick_flat9 + dp.thick_inter9; 
    330                 if (dp.n_shells > 9) 
    331                         out += dp.thick_flat10 + dp.thick_inter10; 
    332                 } 
    333  
    334         return out; 
     528  SphereSLDParameters dp; 
     529 
     530  dp.n_shells = n_shells(); 
     531 
     532  dp.rad_core0 = rad_core0(); 
     533  dp.thick_flat1 = thick_flat1(); 
     534  dp.thick_flat2 = thick_flat2(); 
     535  dp.thick_flat3 = thick_flat3(); 
     536  dp.thick_flat4 = thick_flat4(); 
     537  dp.thick_flat5 = thick_flat5(); 
     538  dp.thick_flat6 = thick_flat6(); 
     539  dp.thick_flat7 = thick_flat7(); 
     540  dp.thick_flat8 = thick_flat8(); 
     541  dp.thick_flat9 = thick_flat9(); 
     542  dp.thick_flat10 = thick_flat10(); 
     543 
     544  dp.thick_inter0 = thick_inter0(); 
     545  dp.thick_inter1 = thick_inter1(); 
     546  dp.thick_inter2 = thick_inter2(); 
     547  dp.thick_inter3 = thick_inter3(); 
     548  dp.thick_inter4 = thick_inter4(); 
     549  dp.thick_inter5 = thick_inter5(); 
     550  dp.thick_inter6 = thick_inter6(); 
     551  dp.thick_inter7 = thick_inter7(); 
     552  dp.thick_inter8 = thick_inter8(); 
     553  dp.thick_inter9 = thick_inter9(); 
     554  dp.thick_inter10 = thick_inter10(); 
     555 
     556  double rad_out = 0.0; 
     557  double out = 0.0; 
     558  // Perform the computation, with all weight points 
     559  double sum = 0.0; 
     560  double norm = 0.0; 
     561 
     562  // Get the dispersion points for the radius 
     563  vector<WeightPoint> weights_rad_core0; 
     564  rad_core0.get_weights(weights_rad_core0); 
     565 
     566  // Get the dispersion points for the thick 1 
     567  vector<WeightPoint> weights_thick_inter0; 
     568  thick_inter0.get_weights(weights_thick_inter0); 
     569  // Loop over radius weight points 
     570  for(size_t i=0; i<weights_rad_core0.size(); i++) { 
     571    dp.rad_core0 = weights_rad_core0[i].value; 
     572    // Loop over radius weight points 
     573    for(size_t j=0; j<weights_thick_inter0.size(); j++) { 
     574      dp.thick_inter0 = weights_thick_inter0[j].value; 
     575      rad_out = dp.rad_core0 + dp.thick_inter0; 
     576      if (dp.n_shells > 0) 
     577        rad_out += dp.thick_flat1 + dp.thick_inter1; 
     578      if (dp.n_shells > 1) 
     579        rad_out += dp.thick_flat2 + dp.thick_inter2; 
     580      if (dp.n_shells > 2) 
     581        rad_out += dp.thick_flat3 + dp.thick_inter3; 
     582      if (dp.n_shells > 3) 
     583        rad_out += dp.thick_flat4 + dp.thick_inter4; 
     584      if (dp.n_shells > 4) 
     585        rad_out += dp.thick_flat5 + dp.thick_inter5; 
     586      if (dp.n_shells > 5) 
     587        rad_out += dp.thick_flat6 + dp.thick_inter6; 
     588      if (dp.n_shells > 6) 
     589        rad_out += dp.thick_flat7 + dp.thick_inter7; 
     590      if (dp.n_shells > 7) 
     591        rad_out += dp.thick_flat8 + dp.thick_inter8; 
     592      if (dp.n_shells > 8) 
     593        rad_out += dp.thick_flat9 + dp.thick_inter9; 
     594      if (dp.n_shells > 9) 
     595        rad_out += dp.thick_flat10 + dp.thick_inter10; 
     596      sum += weights_rad_core0[i].weight*weights_thick_inter0[j].weight 
     597          * (rad_out); 
     598      norm += weights_rad_core0[i].weight*weights_thick_inter0[j].weight; 
     599    } 
     600  } 
     601  if (norm != 0){ 
     602    //return the averaged value 
     603    out =  sum/norm;} 
     604  else{ 
     605    //return normal value 
     606    out = dp.rad_core0 + dp.thick_inter0; 
     607    if (dp.n_shells > 0) 
     608      out += dp.thick_flat1 + dp.thick_inter1; 
     609    if (dp.n_shells > 1) 
     610      out += dp.thick_flat2 + dp.thick_inter2; 
     611    if (dp.n_shells > 2) 
     612      out += dp.thick_flat3 + dp.thick_inter3; 
     613    if (dp.n_shells > 3) 
     614      out += dp.thick_flat4 + dp.thick_inter4; 
     615    if (dp.n_shells > 4) 
     616      out += dp.thick_flat5 + dp.thick_inter5; 
     617    if (dp.n_shells > 5) 
     618      out += dp.thick_flat6 + dp.thick_inter6; 
     619    if (dp.n_shells > 6) 
     620      out += dp.thick_flat7 + dp.thick_inter7; 
     621    if (dp.n_shells > 7) 
     622      out += dp.thick_flat8 + dp.thick_inter8; 
     623    if (dp.n_shells > 8) 
     624      out += dp.thick_flat9 + dp.thick_inter9; 
     625    if (dp.n_shells > 9) 
     626      out += dp.thick_flat10 + dp.thick_inter10; 
     627  } 
     628 
     629  return out; 
    335630 
    336631} 
Note: See TracChangeset for help on using the changeset viewer.