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/triaxialellipsoid.cpp

    r67424cd r82c11d3  
    1818 *   sansmodels/src/libigor 
    1919 * 
    20  *      TODO: refactor so that we pull in the old sansmodels.c_extensions 
    2120 */ 
    2221 
    2322#include <math.h> 
    24 #include "models.hh" 
    2523#include "parameters.hh" 
    2624#include <stdio.h> 
     25#include <stdlib.h> 
    2726using namespace std; 
     27#include "triaxial_ellipsoid.h" 
    2828 
    2929extern "C" { 
    30         #include "libCylinder.h" 
    31         #include "libStructureFactor.h" 
    32         #include "triaxial_ellipsoid.h" 
    33 } 
     30#include "libCylinder.h" 
     31#include "libStructureFactor.h" 
     32} 
     33 
     34typedef struct { 
     35  double scale; 
     36  double semi_axisA; 
     37  double semi_axisB; 
     38  double semi_axisC; 
     39  double sldEll; 
     40  double sldSolv; 
     41  double background; 
     42  double axis_theta; 
     43  double axis_phi; 
     44  double axis_psi; 
     45 
     46} TriaxialEllipsoidParameters; 
     47 
     48static double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) { 
     49  double t,a,b,c; 
     50  double kernel; 
     51 
     52  a = pars->semi_axisA ; 
     53  b = pars->semi_axisB ; 
     54  c = pars->semi_axisC ; 
     55 
     56  t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha)); 
     57  if (t==0.0){ 
     58    kernel  = 1.0; 
     59  }else{ 
     60    kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t); 
     61  } 
     62  return kernel*kernel; 
     63} 
     64 
     65 
     66/** 
     67 * Function to evaluate 2D scattering function 
     68 * @param pars: parameters of the triaxial ellipsoid 
     69 * @param q: q-value 
     70 * @param q_x: q_x / q 
     71 * @param q_y: q_y / q 
     72 * @return: function value 
     73 */ 
     74static double triaxial_ellipsoid_analytical_2D_scaled(TriaxialEllipsoidParameters *pars, double q, double q_x, double q_y) { 
     75  double cyl_x, cyl_y, cyl_z, ell_x, ell_y; 
     76  double q_z; 
     77  double cos_nu,nu; 
     78  double alpha, vol, cos_val; 
     79  double answer; 
     80  double pi = 4.0*atan(1.0); 
     81 
     82  //convert angle degree to radian 
     83  double theta = pars->axis_theta * pi/180.0; 
     84  double phi = pars->axis_phi * pi/180.0; 
     85  double psi = pars->axis_psi * pi/180.0; 
     86 
     87  // Cylinder orientation 
     88  cyl_x = sin(theta) * cos(phi); 
     89  cyl_y = sin(theta) * sin(phi); 
     90  cyl_z = cos(theta); 
     91 
     92  // q vector 
     93  q_z = 0.0; 
     94 
     95  //dx = 1.0; 
     96  //dy = 1.0; 
     97  // Compute the angle btw vector q and the 
     98  // axis of the cylinder 
     99  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
     100 
     101  // The following test should always pass 
     102  if (fabs(cos_val)>1.0) { 
     103    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); 
     104    return 0; 
     105  } 
     106 
     107  // Note: cos(alpha) = 0 and 1 will get an 
     108  // undefined value from CylKernel 
     109  alpha = acos( cos_val ); 
     110 
     111  //ellipse orientation: 
     112  // the elliptical corss section was transformed and projected 
     113  // into the detector plane already through sin(alpha)and furthermore psi remains as same 
     114  // on the detector plane. 
     115  // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt 
     116  // the wave vector q. 
     117 
     118  //x- y- component on the detector plane. 
     119  ell_x =  cos(psi); 
     120  ell_y =  sin(psi); 
     121 
     122  // calculate the axis of the ellipse wrt q-coord. 
     123  cos_nu = ell_x*q_x + ell_y*q_y; 
     124  nu = acos(cos_nu); 
     125 
     126  // Call the IGOR library function to get the kernel 
     127  answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu); 
     128 
     129  // Multiply by contrast^2 
     130  answer *= (pars->sldEll- pars->sldSolv)*(pars->sldEll- pars->sldSolv); 
     131 
     132  //normalize by cylinder volume 
     133  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
     134  vol = 4.0* pi/3.0  * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC; 
     135  answer *= vol; 
     136  //convert to [cm-1] 
     137  answer *= 1.0e8; 
     138  //Scale 
     139  answer *= pars->scale; 
     140 
     141  // add in the background 
     142  answer += pars->background; 
     143 
     144  return answer; 
     145} 
     146 
     147/** 
     148 * Function to evaluate 2D scattering function 
     149 * @param pars: parameters of the triaxial ellipsoid 
     150 * @param q: q-value 
     151 * @return: function value 
     152 */ 
     153static double triaxial_ellipsoid_analytical_2DXY(TriaxialEllipsoidParameters *pars, double qx, double qy) { 
     154  double q; 
     155  q = sqrt(qx*qx+qy*qy); 
     156  return triaxial_ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q); 
     157} 
     158 
     159 
    34160 
    35161TriaxialEllipsoidModel :: TriaxialEllipsoidModel() { 
    36         scale      = Parameter(1.0); 
    37         semi_axisA     = Parameter(35.0, true); 
    38         semi_axisA.set_min(0.0); 
    39         semi_axisB     = Parameter(100.0, true); 
    40         semi_axisB.set_min(0.0); 
    41         semi_axisC  = Parameter(400.0, true); 
    42         semi_axisC.set_min(0.0); 
    43         sldEll   = Parameter(1.0e-6); 
    44         sldSolv   = Parameter(6.3e-6); 
    45         background = Parameter(0.0); 
    46         axis_theta  = Parameter(57.325, true); 
    47         axis_phi    = Parameter(57.325, true); 
    48         axis_psi    = Parameter(0.0, true); 
     162  scale      = Parameter(1.0); 
     163  semi_axisA     = Parameter(35.0, true); 
     164  semi_axisA.set_min(0.0); 
     165  semi_axisB     = Parameter(100.0, true); 
     166  semi_axisB.set_min(0.0); 
     167  semi_axisC  = Parameter(400.0, true); 
     168  semi_axisC.set_min(0.0); 
     169  sldEll   = Parameter(1.0e-6); 
     170  sldSolv   = Parameter(6.3e-6); 
     171  background = Parameter(0.0); 
     172  axis_theta  = Parameter(57.325, true); 
     173  axis_phi    = Parameter(57.325, true); 
     174  axis_psi    = Parameter(0.0, true); 
    49175} 
    50176 
    … …  
    56182 */ 
    57183double TriaxialEllipsoidModel :: operator()(double q) { 
    58         double dp[7]; 
    59  
    60         // Fill parameter array for IGOR library 
    61         // Add the background after averaging 
    62         dp[0] = scale(); 
    63         dp[1] = semi_axisA(); 
    64         dp[2] = semi_axisB(); 
    65         dp[3] = semi_axisC(); 
    66         dp[4] = sldEll(); 
    67         dp[5] = sldSolv(); 
    68         dp[6] = 0.0; 
    69  
    70         // Get the dispersion points for the semi axis A 
    71         vector<WeightPoint> weights_semi_axisA; 
    72         semi_axisA.get_weights(weights_semi_axisA); 
    73  
    74         // Get the dispersion points for the semi axis B 
    75         vector<WeightPoint> weights_semi_axisB; 
    76         semi_axisB.get_weights(weights_semi_axisB); 
    77  
    78         // Get the dispersion points for the semi axis C 
    79         vector<WeightPoint> weights_semi_axisC; 
    80         semi_axisC.get_weights(weights_semi_axisC); 
    81  
    82         // Perform the computation, with all weight points 
    83         double sum = 0.0; 
    84         double norm = 0.0; 
    85         double vol = 0.0; 
    86  
    87         // Loop over semi axis A weight points 
    88         for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
    89                 dp[1] = weights_semi_axisA[i].value; 
    90  
    91                 // Loop over semi axis B weight points 
    92                 for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
    93                         dp[2] = weights_semi_axisB[j].value; 
    94  
    95                         // Loop over semi axis C weight points 
    96                         for(int k=0; k< (int)weights_semi_axisC.size(); k++) { 
    97                                 dp[3] = weights_semi_axisC[k].value; 
    98                                 //Un-normalize  by volume 
    99                                 sum += weights_semi_axisA[i].weight 
    100                                         * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q) 
    101                                         * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
    102                                 //Find average volume 
    103                                 vol += weights_semi_axisA[i].weight 
    104                                         * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight 
    105                                         * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
    106  
    107                                 norm += weights_semi_axisA[i].weight 
    108                                         * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight; 
    109                         } 
    110                 } 
    111         } 
    112         if (vol != 0.0 && norm != 0.0) { 
    113                 //Re-normalize by avg volume 
    114                 sum = sum/(vol/norm);} 
    115  
    116         return sum/norm + background(); 
     184  double dp[7]; 
     185 
     186  // Fill parameter array for IGOR library 
     187  // Add the background after averaging 
     188  dp[0] = scale(); 
     189  dp[1] = semi_axisA(); 
     190  dp[2] = semi_axisB(); 
     191  dp[3] = semi_axisC(); 
     192  dp[4] = sldEll(); 
     193  dp[5] = sldSolv(); 
     194  dp[6] = 0.0; 
     195 
     196  // Get the dispersion points for the semi axis A 
     197  vector<WeightPoint> weights_semi_axisA; 
     198  semi_axisA.get_weights(weights_semi_axisA); 
     199 
     200  // Get the dispersion points for the semi axis B 
     201  vector<WeightPoint> weights_semi_axisB; 
     202  semi_axisB.get_weights(weights_semi_axisB); 
     203 
     204  // Get the dispersion points for the semi axis C 
     205  vector<WeightPoint> weights_semi_axisC; 
     206  semi_axisC.get_weights(weights_semi_axisC); 
     207 
     208  // Perform the computation, with all weight points 
     209  double sum = 0.0; 
     210  double norm = 0.0; 
     211  double vol = 0.0; 
     212 
     213  // Loop over semi axis A weight points 
     214  for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
     215    dp[1] = weights_semi_axisA[i].value; 
     216 
     217    // Loop over semi axis B weight points 
     218    for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
     219      dp[2] = weights_semi_axisB[j].value; 
     220 
     221      // Loop over semi axis C weight points 
     222      for(int k=0; k< (int)weights_semi_axisC.size(); k++) { 
     223        dp[3] = weights_semi_axisC[k].value; 
     224        //Un-normalize  by volume 
     225        sum += weights_semi_axisA[i].weight 
     226            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q) 
     227        * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
     228        //Find average volume 
     229        vol += weights_semi_axisA[i].weight 
     230            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight 
     231            * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
     232 
     233        norm += weights_semi_axisA[i].weight 
     234            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight; 
     235      } 
     236    } 
     237  } 
     238  if (vol != 0.0 && norm != 0.0) { 
     239    //Re-normalize by avg volume 
     240    sum = sum/(vol/norm);} 
     241 
     242  return sum/norm + background(); 
    117243} 
    118244 
    … …  
    124250 */ 
    125251double TriaxialEllipsoidModel :: operator()(double qx, double qy) { 
    126         TriaxialEllipsoidParameters dp; 
    127         // Fill parameter array 
    128         dp.scale      = scale(); 
    129         dp.semi_axisA   = semi_axisA(); 
    130         dp.semi_axisB     = semi_axisB(); 
    131         dp.semi_axisC     = semi_axisC(); 
    132         dp.sldEll   = sldEll(); 
    133         dp.sldSolv   = sldSolv(); 
    134         dp.background = 0.0; 
    135         dp.axis_theta  = axis_theta(); 
    136         dp.axis_phi    = axis_phi(); 
    137         dp.axis_psi    = axis_psi(); 
    138  
    139         // Get the dispersion points for the semi_axis A 
    140         vector<WeightPoint> weights_semi_axisA; 
    141         semi_axisA.get_weights(weights_semi_axisA); 
    142  
    143         // Get the dispersion points for the semi_axis B 
    144         vector<WeightPoint> weights_semi_axisB; 
    145         semi_axisB.get_weights(weights_semi_axisB); 
    146  
    147         // Get the dispersion points for the semi_axis C 
    148         vector<WeightPoint> weights_semi_axisC; 
    149         semi_axisC.get_weights(weights_semi_axisC); 
    150  
    151         // Get angular averaging for theta 
    152         vector<WeightPoint> weights_theta; 
    153         axis_theta.get_weights(weights_theta); 
    154  
    155         // Get angular averaging for phi 
    156         vector<WeightPoint> weights_phi; 
    157         axis_phi.get_weights(weights_phi); 
    158  
    159         // Get angular averaging for psi 
    160         vector<WeightPoint> weights_psi; 
    161         axis_psi.get_weights(weights_psi); 
    162  
    163         // Perform the computation, with all weight points 
    164         double sum = 0.0; 
    165         double norm = 0.0; 
    166         double norm_vol = 0.0; 
    167         double vol = 0.0; 
    168         double pi = 4.0*atan(1.0); 
    169         // Loop over semi axis A weight points 
    170         for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
    171                 dp.semi_axisA = weights_semi_axisA[i].value; 
    172  
    173                 // Loop over semi axis B weight points 
    174                 for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
    175                         dp.semi_axisB = weights_semi_axisB[j].value; 
    176  
    177                         // Loop over semi axis C weight points 
    178                         for(int k=0; k < (int)weights_semi_axisC.size(); k++) { 
    179                         dp.semi_axisC = weights_semi_axisC[k].value; 
    180  
    181                                 // Average over theta distribution 
    182                                 for(int l=0; l< (int)weights_theta.size(); l++) { 
    183                                         dp.axis_theta = weights_theta[l].value; 
    184  
    185                                         // Average over phi distribution 
    186                                         for(int m=0; m <(int)weights_phi.size(); m++) { 
    187                                                 dp.axis_phi = weights_phi[m].value; 
    188                                                 // Average over psi distribution 
    189                                                 for(int n=0; n <(int)weights_psi.size(); n++) { 
    190                                                         dp.axis_psi = weights_psi[n].value; 
    191                                                         //Un-normalize  by volume 
    192                                                         double _ptvalue = weights_semi_axisA[i].weight 
    193                                                                 * weights_semi_axisB[j].weight 
    194                                                                 * weights_semi_axisC[k].weight 
    195                                                                 * weights_theta[l].weight 
    196                                                                 * weights_phi[m].weight 
    197                                                                 * weights_psi[n].weight 
    198                                                                 * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy) 
    199                                                                 * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
    200                                                         if (weights_theta.size()>1) { 
    201                                                                 _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0)); 
    202                                                         } 
    203                                                         sum += _ptvalue; 
    204                                                         //Find average volume 
    205                                                         vol += weights_semi_axisA[i].weight 
    206                                                                 * weights_semi_axisB[j].weight 
    207                                                                 * weights_semi_axisC[k].weight 
    208                                                                 * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
    209                                                         //Find norm for volume 
    210                                                         norm_vol += weights_semi_axisA[i].weight 
    211                                                                 * weights_semi_axisB[j].weight 
    212                                                                 * weights_semi_axisC[k].weight; 
    213  
    214                                                         norm += weights_semi_axisA[i].weight 
    215                                                                 * weights_semi_axisB[j].weight 
    216                                                                 * weights_semi_axisC[k].weight 
    217                                                                 * weights_theta[l].weight 
    218                                                                 * weights_phi[m].weight 
    219                                                                 * weights_psi[n].weight; 
    220                                                 } 
    221                                         } 
    222  
    223                                 } 
    224                         } 
    225                 } 
    226         } 
    227         // Averaging in theta needs an extra normalization 
    228         // factor to account for the sin(theta) term in the 
    229         // integration (see documentation). 
    230         if (weights_theta.size()>1) norm = norm / asin(1.0); 
    231  
    232         if (vol != 0.0 && norm_vol != 0.0) { 
    233                 //Re-normalize by avg volume 
    234                 sum = sum/(vol/norm_vol);} 
    235  
    236         return sum/norm + background(); 
     252  TriaxialEllipsoidParameters dp; 
     253  // Fill parameter array 
     254  dp.scale      = scale(); 
     255  dp.semi_axisA   = semi_axisA(); 
     256  dp.semi_axisB     = semi_axisB(); 
     257  dp.semi_axisC     = semi_axisC(); 
     258  dp.sldEll   = sldEll(); 
     259  dp.sldSolv   = sldSolv(); 
     260  dp.background = 0.0; 
     261  dp.axis_theta  = axis_theta(); 
     262  dp.axis_phi    = axis_phi(); 
     263  dp.axis_psi    = axis_psi(); 
     264 
     265  // Get the dispersion points for the semi_axis A 
     266  vector<WeightPoint> weights_semi_axisA; 
     267  semi_axisA.get_weights(weights_semi_axisA); 
     268 
     269  // Get the dispersion points for the semi_axis B 
     270  vector<WeightPoint> weights_semi_axisB; 
     271  semi_axisB.get_weights(weights_semi_axisB); 
     272 
     273  // Get the dispersion points for the semi_axis C 
     274  vector<WeightPoint> weights_semi_axisC; 
     275  semi_axisC.get_weights(weights_semi_axisC); 
     276 
     277  // Get angular averaging for theta 
     278  vector<WeightPoint> weights_theta; 
     279  axis_theta.get_weights(weights_theta); 
     280 
     281  // Get angular averaging for phi 
     282  vector<WeightPoint> weights_phi; 
     283  axis_phi.get_weights(weights_phi); 
     284 
     285  // Get angular averaging for psi 
     286  vector<WeightPoint> weights_psi; 
     287  axis_psi.get_weights(weights_psi); 
     288 
     289  // Perform the computation, with all weight points 
     290  double sum = 0.0; 
     291  double norm = 0.0; 
     292  double norm_vol = 0.0; 
     293  double vol = 0.0; 
     294  double pi = 4.0*atan(1.0); 
     295  // Loop over semi axis A weight points 
     296  for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
     297    dp.semi_axisA = weights_semi_axisA[i].value; 
     298 
     299    // Loop over semi axis B weight points 
     300    for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
     301      dp.semi_axisB = weights_semi_axisB[j].value; 
     302 
     303      // Loop over semi axis C weight points 
     304      for(int k=0; k < (int)weights_semi_axisC.size(); k++) { 
     305        dp.semi_axisC = weights_semi_axisC[k].value; 
     306 
     307        // Average over theta distribution 
     308        for(int l=0; l< (int)weights_theta.size(); l++) { 
     309          dp.axis_theta = weights_theta[l].value; 
     310 
     311          // Average over phi distribution 
     312          for(int m=0; m <(int)weights_phi.size(); m++) { 
     313            dp.axis_phi = weights_phi[m].value; 
     314            // Average over psi distribution 
     315            for(int n=0; n <(int)weights_psi.size(); n++) { 
     316              dp.axis_psi = weights_psi[n].value; 
     317              //Un-normalize  by volume 
     318              double _ptvalue = weights_semi_axisA[i].weight 
     319                  * weights_semi_axisB[j].weight 
     320                  * weights_semi_axisC[k].weight 
     321                  * weights_theta[l].weight 
     322                  * weights_phi[m].weight 
     323                  * weights_psi[n].weight 
     324                  * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy) 
     325              * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
     326              if (weights_theta.size()>1) { 
     327                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0)); 
     328              } 
     329              sum += _ptvalue; 
     330              //Find average volume 
     331              vol += weights_semi_axisA[i].weight 
     332                  * weights_semi_axisB[j].weight 
     333                  * weights_semi_axisC[k].weight 
     334                  * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value; 
     335              //Find norm for volume 
     336              norm_vol += weights_semi_axisA[i].weight 
     337                  * weights_semi_axisB[j].weight 
     338                  * weights_semi_axisC[k].weight; 
     339 
     340              norm += weights_semi_axisA[i].weight 
     341                  * weights_semi_axisB[j].weight 
     342                  * weights_semi_axisC[k].weight 
     343                  * weights_theta[l].weight 
     344                  * weights_phi[m].weight 
     345                  * weights_psi[n].weight; 
     346            } 
     347          } 
     348 
     349        } 
     350      } 
     351    } 
     352  } 
     353  // Averaging in theta needs an extra normalization 
     354  // factor to account for the sin(theta) term in the 
     355  // integration (see documentation). 
     356  if (weights_theta.size()>1) norm = norm / asin(1.0); 
     357 
     358  if (vol != 0.0 && norm_vol != 0.0) { 
     359    //Re-normalize by avg volume 
     360    sum = sum/(vol/norm_vol);} 
     361 
     362  return sum/norm + background(); 
    237363} 
    238364 
    … …  
    245371 */ 
    246372double TriaxialEllipsoidModel :: evaluate_rphi(double q, double phi) { 
    247         double qx = q*cos(phi); 
    248         double qy = q*sin(phi); 
    249         return (*this).operator()(qx, qy); 
     373  double qx = q*cos(phi); 
     374  double qy = q*sin(phi); 
     375  return (*this).operator()(qx, qy); 
    250376} 
    251377/** 
    … …  
    254380 */ 
    255381double TriaxialEllipsoidModel :: calculate_ER() { 
    256         TriaxialEllipsoidParameters dp; 
    257  
    258         dp.semi_axisA   = semi_axisA(); 
    259         dp.semi_axisB     = semi_axisB(); 
    260         //polar axis C 
    261         dp.semi_axisC     = semi_axisC(); 
    262  
    263         double rad_out = 0.0; 
    264         //Surface average radius at the equat. cross section. 
    265         double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB); 
    266  
    267         // Perform the computation, with all weight points 
    268         double sum = 0.0; 
    269         double norm = 0.0; 
    270  
    271         // Get the dispersion points for the semi_axis A 
    272         vector<WeightPoint> weights_semi_axisA; 
    273         semi_axisA.get_weights(weights_semi_axisA); 
    274  
    275         // Get the dispersion points for the semi_axis B 
    276         vector<WeightPoint> weights_semi_axisB; 
    277         semi_axisB.get_weights(weights_semi_axisB); 
    278  
    279         // Get the dispersion points for the semi_axis C 
    280         vector<WeightPoint> weights_semi_axisC; 
    281         semi_axisC.get_weights(weights_semi_axisC); 
    282  
    283         // Loop over semi axis A weight points 
    284         for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
    285                 dp.semi_axisA = weights_semi_axisA[i].value; 
    286  
    287                 // Loop over semi axis B weight points 
    288                 for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
    289                         dp.semi_axisB = weights_semi_axisB[j].value; 
    290  
    291                         // Loop over semi axis C weight points 
    292                         for(int k=0; k < (int)weights_semi_axisC.size(); k++) { 
    293                                 dp.semi_axisC = weights_semi_axisC[k].value; 
    294  
    295                                 //Calculate surface averaged radius 
    296                                 suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB); 
    297  
    298                                 //Sum 
    299                                 sum += weights_semi_axisA[i].weight 
    300                                         * weights_semi_axisB[j].weight 
    301                                         * weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0; 
    302                                 //Norm 
    303                                 norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight 
    304                                         * weights_semi_axisC[k].weight; 
    305                         } 
    306                 } 
    307         } 
    308         if (norm != 0){ 
    309                 //return the averaged value 
    310                 rad_out =  sum/norm;} 
    311         else{ 
    312                 //return normal value 
    313                 rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;} 
    314  
    315         return rad_out; 
    316 } 
     382  TriaxialEllipsoidParameters dp; 
     383 
     384  dp.semi_axisA   = semi_axisA(); 
     385  dp.semi_axisB     = semi_axisB(); 
     386  //polar axis C 
     387  dp.semi_axisC     = semi_axisC(); 
     388 
     389  double rad_out = 0.0; 
     390  //Surface average radius at the equat. cross section. 
     391  double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB); 
     392 
     393  // Perform the computation, with all weight points 
     394  double sum = 0.0; 
     395  double norm = 0.0; 
     396 
     397  // Get the dispersion points for the semi_axis A 
     398  vector<WeightPoint> weights_semi_axisA; 
     399  semi_axisA.get_weights(weights_semi_axisA); 
     400 
     401  // Get the dispersion points for the semi_axis B 
     402  vector<WeightPoint> weights_semi_axisB; 
     403  semi_axisB.get_weights(weights_semi_axisB); 
     404 
     405  // Get the dispersion points for the semi_axis C 
     406  vector<WeightPoint> weights_semi_axisC; 
     407  semi_axisC.get_weights(weights_semi_axisC); 
     408 
     409  // Loop over semi axis A weight points 
     410  for(int i=0; i< (int)weights_semi_axisA.size(); i++) { 
     411    dp.semi_axisA = weights_semi_axisA[i].value; 
     412 
     413    // Loop over semi axis B weight points 
     414    for(int j=0; j< (int)weights_semi_axisB.size(); j++) { 
     415      dp.semi_axisB = weights_semi_axisB[j].value; 
     416 
     417      // Loop over semi axis C weight points 
     418      for(int k=0; k < (int)weights_semi_axisC.size(); k++) { 
     419        dp.semi_axisC = weights_semi_axisC[k].value; 
     420 
     421        //Calculate surface averaged radius 
     422        suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB); 
     423 
     424        //Sum 
     425        sum += weights_semi_axisA[i].weight 
     426            * weights_semi_axisB[j].weight 
     427            * weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0; 
     428        //Norm 
     429        norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight 
     430            * weights_semi_axisC[k].weight; 
     431      } 
     432    } 
     433  } 
     434  if (norm != 0){ 
     435    //return the averaged value 
     436    rad_out =  sum/norm;} 
     437  else{ 
     438    //return normal value 
     439    rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;} 
     440 
     441  return rad_out; 
     442} 
Note: See TracChangeset for help on using the changeset viewer.