Ignore:
Timestamp:
Jan 5, 2012 6:24:51 PM (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:
20d91bd
Parents:
98fdccd
Message:

refactored bunch of models

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/c_models/stackeddisks.cpp

    r67424cd r82c11d3  
    2323 
    2424#include <math.h> 
    25 #include "models.hh" 
    2625#include "parameters.hh" 
    2726#include <stdio.h> 
    2827using namespace std; 
     28#include "stacked_disks.h" 
    2929 
    3030extern "C" { 
    31         #include "libCylinder.h" 
    32         #include "libStructureFactor.h" 
    33         #include "stacked_disks.h" 
     31#include "libCylinder.h" 
     32#include "libStructureFactor.h" 
     33} 
     34 
     35typedef struct { 
     36  double scale; 
     37  double radius; 
     38  double core_thick; 
     39  double layer_thick; 
     40  double core_sld; 
     41  double layer_sld; 
     42  double solvent_sld; 
     43  double n_stacking; 
     44  double sigma_d; 
     45  double background; 
     46  double axis_theta; 
     47  double axis_phi; 
     48} StackedDisksParameters; 
     49 
     50 
     51/** 
     52 * Function to evaluate 2D scattering function 
     53 * @param pars: parameters of the staked disks 
     54 * @param q: q-value 
     55 * @param q_x: q_x / q 
     56 * @param q_y: q_y / q 
     57 * @return: function value 
     58 */ 
     59static double stacked_disks_analytical_2D_scaled(StackedDisksParameters *pars, double q, double q_x, double q_y) { 
     60  double cyl_x, cyl_y, cyl_z; 
     61  double q_z; 
     62  double alpha, vol, cos_val; 
     63  double d, dum, halfheight; 
     64  double answer; 
     65 
     66 
     67 
     68  // parallelepiped orientation 
     69  cyl_x = sin(pars->axis_theta) * cos(pars->axis_phi); 
     70  cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi); 
     71  cyl_z = cos(pars->axis_theta); 
     72 
     73  // q vector 
     74  q_z = 0; 
     75 
     76  // Compute the angle btw vector q and the 
     77  // axis of the parallelepiped 
     78  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
     79 
     80  // The following test should always pass 
     81  if (fabs(cos_val)>1.0) { 
     82    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     83    return 0; 
     84  } 
     85 
     86  // Note: cos(alpha) = 0 and 1 will get an 
     87  // undefined value from Stackdisc_kern 
     88  alpha = acos( cos_val ); 
     89 
     90  // Call the IGOR library function to get the kernel 
     91  d = 2 * pars->layer_thick + pars->core_thick; 
     92  halfheight = pars->core_thick/2.0; 
     93  dum =alpha ; 
     94  answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld, 
     95      pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha); 
     96 
     97  // Multiply by contrast^2 
     98  //answer *= pars->contrast*pars->contrast; 
     99 
     100  //normalize by staked disks volume 
     101  vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking; 
     102  answer /= vol; 
     103 
     104  //convert to [cm-1] 
     105  answer *= 1.0e8; 
     106 
     107  //Scale 
     108  answer *= pars->scale; 
     109 
     110  // add in the background 
     111  answer += pars->background; 
     112 
     113  return answer; 
     114} 
     115 
     116/** 
     117 * Function to evaluate 2D scattering function 
     118 * @param pars: parameters of the staked disks 
     119 * @param q: q-value 
     120 * @return: function value 
     121 */ 
     122static double stacked_disks_analytical_2DXY(StackedDisksParameters *pars, double qx, double qy) { 
     123  double q; 
     124  q = sqrt(qx*qx+qy*qy); 
     125  return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    34126} 
    35127 
    36128StackedDisksModel :: StackedDisksModel() { 
    37         scale      = Parameter(1.0); 
    38         radius     = Parameter(3000.0, true); 
    39         radius.set_min(0.0); 
    40         core_thick  = Parameter(10.0, true); 
    41         core_thick.set_min(0.0); 
    42         layer_thick     = Parameter(15.0); 
    43         layer_thick.set_min(0.0); 
    44         core_sld = Parameter(4.0e-6); 
    45         layer_sld  = Parameter(-4.0e-7); 
    46         solvent_sld  = Parameter(5.0e-6); 
    47         n_stacking   = Parameter(1); 
    48         sigma_d   = Parameter(0); 
    49         background = Parameter(0.001); 
    50         axis_theta  = Parameter(0.0, true); 
    51         axis_phi    = Parameter(0.0, true); 
     129  scale      = Parameter(1.0); 
     130  radius     = Parameter(3000.0, true); 
     131  radius.set_min(0.0); 
     132  core_thick  = Parameter(10.0, true); 
     133  core_thick.set_min(0.0); 
     134  layer_thick     = Parameter(15.0); 
     135  layer_thick.set_min(0.0); 
     136  core_sld = Parameter(4.0e-6); 
     137  layer_sld  = Parameter(-4.0e-7); 
     138  solvent_sld  = Parameter(5.0e-6); 
     139  n_stacking   = Parameter(1); 
     140  sigma_d   = Parameter(0); 
     141  background = Parameter(0.001); 
     142  axis_theta  = Parameter(0.0, true); 
     143  axis_phi    = Parameter(0.0, true); 
    52144} 
    53145 
    … …  
    59151 */ 
    60152double StackedDisksModel :: operator()(double q) { 
    61         double dp[10]; 
    62  
    63         // Fill parameter array for IGOR library 
    64         // Add the background after averaging 
    65         dp[0] = scale(); 
    66         dp[1] = radius(); 
    67         dp[2] = core_thick(); 
    68         dp[3] = layer_thick(); 
    69         dp[4] = core_sld(); 
    70         dp[5] = layer_sld(); 
    71         dp[6] = solvent_sld(); 
    72         dp[7] = n_stacking(); 
    73         dp[8] = sigma_d(); 
    74         dp[9] = 0.0; 
    75  
    76         // Get the dispersion points for the radius 
    77         vector<WeightPoint> weights_radius; 
    78         radius.get_weights(weights_radius); 
    79  
    80         // Get the dispersion points for the core_thick 
    81         vector<WeightPoint> weights_core_thick; 
    82         core_thick.get_weights(weights_core_thick); 
    83  
    84         // Get the dispersion points for the layer_thick 
    85         vector<WeightPoint> weights_layer_thick; 
    86         layer_thick.get_weights(weights_layer_thick); 
    87  
    88         // Perform the computation, with all weight points 
    89         double sum = 0.0; 
    90         double norm = 0.0; 
    91         double vol = 0.0; 
    92  
    93         // Loop over length weight points 
    94         for(int i=0; i< (int)weights_radius.size(); i++) { 
    95                 dp[1] = weights_radius[i].value; 
    96  
    97                 // Loop over radius weight points 
    98                 for(int j=0; j< (int)weights_core_thick.size(); j++) { 
    99                         dp[2] = weights_core_thick[j].value; 
    100  
    101                         // Loop over thickness weight points 
    102                         for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
    103                                 dp[3] = weights_layer_thick[k].value; 
    104                                 //Un-normalize by volume 
    105                                 sum += weights_radius[i].weight 
    106                                         * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q) 
    107                                         *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value); 
    108                                 //Find average volume 
    109                                 vol += weights_radius[i].weight 
    110                                         * weights_core_thick[j].weight * weights_layer_thick[k].weight 
    111                                         *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value); 
    112                                 norm += weights_radius[i].weight 
    113                                         * weights_core_thick[j].weight* weights_layer_thick[k].weight; 
    114                         } 
    115                 } 
    116         } 
    117         if (vol != 0.0 && norm != 0.0) { 
    118                 //Re-normalize by avg volume 
    119                 sum = sum/(vol/norm);} 
    120  
    121         return sum/norm + background(); 
     153  double dp[10]; 
     154 
     155  // Fill parameter array for IGOR library 
     156  // Add the background after averaging 
     157  dp[0] = scale(); 
     158  dp[1] = radius(); 
     159  dp[2] = core_thick(); 
     160  dp[3] = layer_thick(); 
     161  dp[4] = core_sld(); 
     162  dp[5] = layer_sld(); 
     163  dp[6] = solvent_sld(); 
     164  dp[7] = n_stacking(); 
     165  dp[8] = sigma_d(); 
     166  dp[9] = 0.0; 
     167 
     168  // Get the dispersion points for the radius 
     169  vector<WeightPoint> weights_radius; 
     170  radius.get_weights(weights_radius); 
     171 
     172  // Get the dispersion points for the core_thick 
     173  vector<WeightPoint> weights_core_thick; 
     174  core_thick.get_weights(weights_core_thick); 
     175 
     176  // Get the dispersion points for the layer_thick 
     177  vector<WeightPoint> weights_layer_thick; 
     178  layer_thick.get_weights(weights_layer_thick); 
     179 
     180  // Perform the computation, with all weight points 
     181  double sum = 0.0; 
     182  double norm = 0.0; 
     183  double vol = 0.0; 
     184 
     185  // Loop over length weight points 
     186  for(int i=0; i< (int)weights_radius.size(); i++) { 
     187    dp[1] = weights_radius[i].value; 
     188 
     189    // Loop over radius weight points 
     190    for(int j=0; j< (int)weights_core_thick.size(); j++) { 
     191      dp[2] = weights_core_thick[j].value; 
     192 
     193      // Loop over thickness weight points 
     194      for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
     195        dp[3] = weights_layer_thick[k].value; 
     196        //Un-normalize by volume 
     197        sum += weights_radius[i].weight 
     198            * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q) 
     199        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value); 
     200        //Find average volume 
     201        vol += weights_radius[i].weight 
     202            * weights_core_thick[j].weight * weights_layer_thick[k].weight 
     203            *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value); 
     204        norm += weights_radius[i].weight 
     205            * weights_core_thick[j].weight* weights_layer_thick[k].weight; 
     206      } 
     207    } 
     208  } 
     209  if (vol != 0.0 && norm != 0.0) { 
     210    //Re-normalize by avg volume 
     211    sum = sum/(vol/norm);} 
     212 
     213  return sum/norm + background(); 
    122214} 
    123215 
    … …  
    129221 */ 
    130222double StackedDisksModel :: operator()(double qx, double qy) { 
    131         StackedDisksParameters dp; 
    132         // Fill parameter array 
    133         dp.scale      = scale(); 
    134         dp.core_thick    = core_thick(); 
    135         dp.radius         = radius(); 
    136         dp.layer_thick  = layer_thick(); 
    137         dp.core_sld   = core_sld(); 
    138         dp.layer_sld  = layer_sld(); 
    139         dp.solvent_sld= solvent_sld(); 
    140         dp.n_stacking     = n_stacking(); 
    141         dp.sigma_d   = sigma_d(); 
    142         dp.background = 0.0; 
    143         dp.axis_theta = axis_theta(); 
    144         dp.axis_phi   = axis_phi(); 
    145  
    146         // Get the dispersion points for the length 
    147         vector<WeightPoint> weights_core_thick; 
    148         core_thick.get_weights(weights_core_thick); 
    149  
    150         // Get the dispersion points for the radius 
    151         vector<WeightPoint> weights_radius; 
    152         radius.get_weights(weights_radius); 
    153  
    154         // Get the dispersion points for the thickness 
    155         vector<WeightPoint> weights_layer_thick; 
    156         layer_thick.get_weights(weights_layer_thick); 
    157  
    158         // Get angular averaging for theta 
    159         vector<WeightPoint> weights_theta; 
    160         axis_theta.get_weights(weights_theta); 
    161  
    162         // Get angular averaging for phi 
    163         vector<WeightPoint> weights_phi; 
    164         axis_phi.get_weights(weights_phi); 
    165  
    166         // Perform the computation, with all weight points 
    167         double sum = 0.0; 
    168         double norm = 0.0; 
    169         double norm_vol = 0.0; 
    170         double vol = 0.0; 
    171  
    172         // Loop over length weight points 
    173         for(int i=0; i< (int)weights_core_thick.size(); i++) { 
    174                 dp.core_thick = weights_core_thick[i].value; 
    175  
    176                 // Loop over radius weight points 
    177                 for(int j=0; j< (int)weights_radius.size(); j++) { 
    178                         dp.radius = weights_radius[j].value; 
    179  
    180                                 // Loop over thickness weight points 
    181                                 for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
    182                                 dp.layer_thick = weights_layer_thick[k].value; 
    183  
    184                                         for(int l=0; l< (int)weights_theta.size(); l++) { 
    185                                         dp.axis_theta = weights_theta[l].value; 
    186  
    187                                                 // Average over phi distribution 
    188                                                 for(int m=0; m <(int)weights_phi.size(); m++) { 
    189                                                         dp.axis_phi = weights_phi[m].value; 
    190  
    191                                                         //Un-normalize by volume 
    192                                                         double _ptvalue = weights_core_thick[i].weight 
    193                                                                 * weights_radius[j].weight 
    194                                                                 * weights_layer_thick[k].weight 
    195                                                                 * weights_theta[l].weight 
    196                                                                 * weights_phi[m].weight 
    197                                                                 * stacked_disks_analytical_2DXY(&dp, qx, qy) 
    198                                                                 *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value); 
    199                                                         if (weights_theta.size()>1) { 
    200                                                                 _ptvalue *= fabs(sin(weights_theta[l].value)); 
    201                                                         } 
    202                                                         sum += _ptvalue; 
    203                                                         //Find average volume 
    204                                                         vol += weights_radius[j].weight 
    205                                                                 * weights_core_thick[i].weight * weights_layer_thick[k].weight 
    206                                                                 *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value); 
    207                                                         //Find norm for volume 
    208                                                         norm_vol += weights_radius[j].weight 
    209                                                                 * weights_core_thick[i].weight * weights_layer_thick[k].weight; 
    210  
    211                                                         norm += weights_core_thick[i].weight 
    212                                                                 * weights_radius[j].weight 
    213                                                                 * weights_layer_thick[k].weight 
    214                                                                 * weights_theta[l].weight 
    215                                                                 * weights_phi[m].weight; 
    216                                                 } 
    217                                 } 
    218                         } 
    219                 } 
    220         } 
    221         // Averaging in theta needs an extra normalization 
    222         // factor to account for the sin(theta) term in the 
    223         // integration (see documentation). 
    224         if (weights_theta.size()>1) norm = norm / asin(1.0); 
    225         if (vol != 0.0 && norm_vol != 0.0) { 
    226                 //Re-normalize by avg volume 
    227                 sum = sum/(vol/norm_vol);} 
    228         return sum/norm + background(); 
     223  StackedDisksParameters dp; 
     224  // Fill parameter array 
     225  dp.scale      = scale(); 
     226  dp.core_thick    = core_thick(); 
     227  dp.radius       = radius(); 
     228  dp.layer_thick  = layer_thick(); 
     229  dp.core_sld   = core_sld(); 
     230  dp.layer_sld  = layer_sld(); 
     231  dp.solvent_sld= solvent_sld(); 
     232  dp.n_stacking   = n_stacking(); 
     233  dp.sigma_d   = sigma_d(); 
     234  dp.background = 0.0; 
     235  dp.axis_theta = axis_theta(); 
     236  dp.axis_phi   = axis_phi(); 
     237 
     238  // Get the dispersion points for the length 
     239  vector<WeightPoint> weights_core_thick; 
     240  core_thick.get_weights(weights_core_thick); 
     241 
     242  // Get the dispersion points for the radius 
     243  vector<WeightPoint> weights_radius; 
     244  radius.get_weights(weights_radius); 
     245 
     246  // Get the dispersion points for the thickness 
     247  vector<WeightPoint> weights_layer_thick; 
     248  layer_thick.get_weights(weights_layer_thick); 
     249 
     250  // Get angular averaging for theta 
     251  vector<WeightPoint> weights_theta; 
     252  axis_theta.get_weights(weights_theta); 
     253 
     254  // Get angular averaging for phi 
     255  vector<WeightPoint> weights_phi; 
     256  axis_phi.get_weights(weights_phi); 
     257 
     258  // Perform the computation, with all weight points 
     259  double sum = 0.0; 
     260  double norm = 0.0; 
     261  double norm_vol = 0.0; 
     262  double vol = 0.0; 
     263 
     264  // Loop over length weight points 
     265  for(int i=0; i< (int)weights_core_thick.size(); i++) { 
     266    dp.core_thick = weights_core_thick[i].value; 
     267 
     268    // Loop over radius weight points 
     269    for(int j=0; j< (int)weights_radius.size(); j++) { 
     270      dp.radius = weights_radius[j].value; 
     271 
     272      // Loop over thickness weight points 
     273      for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
     274        dp.layer_thick = weights_layer_thick[k].value; 
     275 
     276        for(int l=0; l< (int)weights_theta.size(); l++) { 
     277          dp.axis_theta = weights_theta[l].value; 
     278 
     279          // Average over phi distribution 
     280          for(int m=0; m <(int)weights_phi.size(); m++) { 
     281            dp.axis_phi = weights_phi[m].value; 
     282 
     283            //Un-normalize by volume 
     284            double _ptvalue = weights_core_thick[i].weight 
     285                * weights_radius[j].weight 
     286                * weights_layer_thick[k].weight 
     287                * weights_theta[l].weight 
     288                * weights_phi[m].weight 
     289                * stacked_disks_analytical_2DXY(&dp, qx, qy) 
     290            *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value); 
     291            if (weights_theta.size()>1) { 
     292              _ptvalue *= fabs(sin(weights_theta[l].value)); 
     293            } 
     294            sum += _ptvalue; 
     295            //Find average volume 
     296            vol += weights_radius[j].weight 
     297                * weights_core_thick[i].weight * weights_layer_thick[k].weight 
     298                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value); 
     299            //Find norm for volume 
     300            norm_vol += weights_radius[j].weight 
     301                * weights_core_thick[i].weight * weights_layer_thick[k].weight; 
     302 
     303            norm += weights_core_thick[i].weight 
     304                * weights_radius[j].weight 
     305                * weights_layer_thick[k].weight 
     306                * weights_theta[l].weight 
     307                * weights_phi[m].weight; 
     308          } 
     309        } 
     310      } 
     311    } 
     312  } 
     313  // Averaging in theta needs an extra normalization 
     314  // factor to account for the sin(theta) term in the 
     315  // integration (see documentation). 
     316  if (weights_theta.size()>1) norm = norm / asin(1.0); 
     317  if (vol != 0.0 && norm_vol != 0.0) { 
     318    //Re-normalize by avg volume 
     319    sum = sum/(vol/norm_vol);} 
     320  return sum/norm + background(); 
    229321} 
    230322 
    … …  
    237329 */ 
    238330double StackedDisksModel :: evaluate_rphi(double q, double phi) { 
    239         double qx = q*cos(phi); 
    240         double qy = q*sin(phi); 
    241         return (*this).operator()(qx, qy); 
     331  double qx = q*cos(phi); 
     332  double qy = q*sin(phi); 
     333  return (*this).operator()(qx, qy); 
    242334} 
    243335/** 
    … …  
    246338 */ 
    247339double StackedDisksModel :: calculate_ER() { 
    248         StackedDisksParameters dp; 
    249  
    250         dp.core_thick    = core_thick(); 
    251         dp.radius         = radius(); 
    252         dp.layer_thick  = layer_thick(); 
    253         dp.n_stacking     = n_stacking(); 
    254  
    255         double rad_out = 0.0; 
    256         if (dp.n_stacking <= 0.0){ 
    257                 return rad_out; 
    258         } 
    259  
    260         // Perform the computation, with all weight points 
    261         double sum = 0.0; 
    262         double norm = 0.0; 
    263  
    264         // Get the dispersion points for the length 
    265         vector<WeightPoint> weights_core_thick; 
    266         core_thick.get_weights(weights_core_thick); 
    267  
    268         // Get the dispersion points for the radius 
    269         vector<WeightPoint> weights_radius; 
    270         radius.get_weights(weights_radius); 
    271  
    272         // Get the dispersion points for the thickness 
    273         vector<WeightPoint> weights_layer_thick; 
    274         layer_thick.get_weights(weights_layer_thick); 
    275  
    276         // Loop over major shell weight points 
    277         for(int i=0; i< (int)weights_core_thick.size(); i++) { 
    278                 dp.core_thick = weights_core_thick[i].value; 
    279                 for(int j=0; j< (int)weights_layer_thick.size(); j++) { 
    280                         dp.layer_thick = weights_layer_thick[j].value; 
    281                         for(int k=0; k< (int)weights_radius.size(); k++) { 
    282                                 dp.radius = weights_radius[k].value; 
    283                                 //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. 
    284                                 sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight 
    285                                         * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0; 
    286                                 norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight; 
    287                         } 
    288                 } 
    289         } 
    290         if (norm != 0){ 
    291                 //return the averaged value 
    292                 rad_out =  sum/norm;} 
    293         else{ 
    294                 //return normal value 
    295                 //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. 
    296                 rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;} 
    297  
    298         return rad_out; 
    299 } 
     340  StackedDisksParameters dp; 
     341 
     342  dp.core_thick    = core_thick(); 
     343  dp.radius       = radius(); 
     344  dp.layer_thick  = layer_thick(); 
     345  dp.n_stacking   = n_stacking(); 
     346 
     347  double rad_out = 0.0; 
     348  if (dp.n_stacking <= 0.0){ 
     349    return rad_out; 
     350  } 
     351 
     352  // Perform the computation, with all weight points 
     353  double sum = 0.0; 
     354  double norm = 0.0; 
     355 
     356  // Get the dispersion points for the length 
     357  vector<WeightPoint> weights_core_thick; 
     358  core_thick.get_weights(weights_core_thick); 
     359 
     360  // Get the dispersion points for the radius 
     361  vector<WeightPoint> weights_radius; 
     362  radius.get_weights(weights_radius); 
     363 
     364  // Get the dispersion points for the thickness 
     365  vector<WeightPoint> weights_layer_thick; 
     366  layer_thick.get_weights(weights_layer_thick); 
     367 
     368  // Loop over major shell weight points 
     369  for(int i=0; i< (int)weights_core_thick.size(); i++) { 
     370    dp.core_thick = weights_core_thick[i].value; 
     371    for(int j=0; j< (int)weights_layer_thick.size(); j++) { 
     372      dp.layer_thick = weights_layer_thick[j].value; 
     373      for(int k=0; k< (int)weights_radius.size(); k++) { 
     374        dp.radius = weights_radius[k].value; 
     375        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. 
     376        sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight 
     377            * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0; 
     378        norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight; 
     379      } 
     380    } 
     381  } 
     382  if (norm != 0){ 
     383    //return the averaged value 
     384    rad_out =  sum/norm;} 
     385  else{ 
     386    //return normal value 
     387    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. 
     388    rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;} 
     389 
     390  return rad_out; 
     391} 
Note: See TracChangeset for help on using the changeset viewer.