Changes in src/sas/models/c_extension/c_models/barbell.cpp [db46d13:79492222] in sasview
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/models/c_extension/c_models/barbell.cpp
rdb46d13 r79492222 30 30 #include "barbell.h" 31 31 } 32 33 // Convenience parameter structure34 typedef struct {35 double scale;36 double rad_bar;37 double len_bar;38 double rad_bell;39 double sld_barbell;40 double sld_solv;41 double background;42 double theta;43 double phi;44 } BarBellParameters;45 32 46 33 BarBellModel :: BarBellModel() { … … 59 46 } 60 47 61 double bar bell2d_kernel(double dp[], double q, double alpha) {48 double bar2d_kernel(double dp[], double q, double alpha) { 62 49 int j; 63 50 double Pi; … … 126 113 return answer; 127 114 } 128 129 /**130 * Function to evaluate 2D scattering function131 * @param pars: parameters of the BarBell132 * @param q: q-value133 * @param q_x: q_x / q134 * @param q_y: q_y / q135 * @return: function value136 */137 static double barbell_analytical_2D_scaled(BarBellParameters *pars, double q, double q_x, double q_y) {138 double cyl_x, cyl_y;//, cyl_z;139 //double q_z;140 double alpha, cos_val;141 double answer;142 double dp[7];143 //convert angle degree to radian144 double pi = 4.0*atan(1.0);145 double theta = pars->theta * pi/180.0;146 double phi = pars->phi * pi/180.0;147 148 dp[0] = pars->scale;149 dp[1] = pars->rad_bar;150 dp[2] = pars->len_bar;151 dp[3] = pars->rad_bell;152 dp[4] = pars->sld_barbell;153 dp[5] = pars->sld_solv;154 dp[6] = pars->background;155 156 157 //double Pi = 4.0*atan(1.0);158 // Cylinder orientation159 cyl_x = cos(theta) * cos(phi);160 cyl_y = sin(theta);161 //cyl_z = -cos(theta) * sin(phi);162 163 // q vector164 //q_z = 0;165 166 // Compute the angle btw vector q and the167 // axis of the cylinder168 cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;169 170 // The following test should always pass171 if (fabs(cos_val)>1.0) {172 printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");173 return 0;174 }175 176 // Note: cos(alpha) = 0 and 1 will get an177 // undefined value from CylKernel178 alpha = acos( cos_val );179 180 answer = barbell2d_kernel(dp, q, alpha);181 182 183 return answer;184 185 }186 187 /**188 * Function to evaluate 2D scattering function189 * @param pars: parameters of the BarBell190 * @param q: q-value191 * @return: function value192 */193 static double barbell_analytical_2DXY(BarBellParameters *pars, double qx, double qy){194 double q;195 q = sqrt(qx*qx+qy*qy);196 return barbell_analytical_2D_scaled(pars, q, qx/q, qy/q);197 }198 199 115 /** 200 116 * Function to evaluate 1D scattering function … … 274 190 */ 275 191 double BarBellModel :: operator()(double qx, double qy) { 276 BarBellParameters dp; 277 278 // Fill parameter array for IGOR library 279 // Add the background after averaging 280 dp.scale = scale(); 281 dp.rad_bar = rad_bar(); 282 dp.len_bar = len_bar(); 283 dp.rad_bell = rad_bell(); 284 dp.sld_barbell = sld_barbell(); 285 dp.sld_solv = sld_solv(); 286 dp.background = 0.0; 287 dp.theta = theta(); 288 dp.phi = phi(); 192 double dp[7]; 193 194 // Fill parameter array for IGOR library 195 // Add the background after averaging 196 dp[0] = scale(); 197 dp[1] = rad_bar(); 198 dp[2] = len_bar(); 199 dp[3] = rad_bell(); 200 dp[4] = sld_barbell(); 201 dp[5] = sld_solv(); 202 dp[6] = 0.0; 203 204 double _theta = theta(); 205 double _phi = phi(); 289 206 290 207 // Get the dispersion points for the rad_bar … … 320 237 // Loop over radius weight points 321 238 for(size_t i=0; i<weights_rad_bar.size(); i++) { 322 dp .rad_bar= weights_rad_bar[i].value;239 dp[1] = weights_rad_bar[i].value; 323 240 for(size_t j=0; j<weights_len_bar.size(); j++) { 324 dp .len_bar= weights_len_bar[j].value;241 dp[2] = weights_len_bar[j].value; 325 242 for(size_t k=0; k<weights_rad_bell.size(); k++) { 326 dp .rad_bell= weights_rad_bell[k].value;243 dp[3] = weights_rad_bell[k].value; 327 244 // Average over theta distribution 328 245 for(size_t l=0; l< weights_theta.size(); l++) { 329 dp.theta = weights_theta[l].value;246 _theta = weights_theta[l].value; 330 247 // Average over phi distribution 331 248 for(size_t m=0; m< weights_phi.size(); m++) { 332 dp.phi = weights_phi[m].value;249 _phi = weights_phi[m].value; 333 250 //Un-normalize Form by volume 334 hDist = -sqrt(fabs(dp.rad_bell*dp.rad_bell-dp.rad_bar*dp.rad_bar)); 335 vol_i = pi*dp.rad_bar*dp.rad_bar*dp.len_bar+2.0*pi/3.0*((dp.rad_bell-hDist)*(dp.rad_bell-hDist)* 336 (2*dp.rad_bell+hDist)); 251 hDist = sqrt(fabs(dp[3]*dp[3]-dp[1]*dp[1])); 252 vol_i = pi*dp[1]*dp[1]*dp[2]+2.0*pi*(2.0*dp[3]*dp[3]*dp[3]/3.0 253 +dp[3]*dp[3]*hDist-hDist*hDist*hDist/3.0); 254 255 const double q = sqrt(qx*qx+qy*qy); 256 //convert angle degree to radian 257 const double pi = 4.0*atan(1.0); 258 259 // Cylinder orientation 260 const double cyl_x = cos(_theta * pi/180.0) * cos(_phi * pi/180.0); 261 const double cyl_y = sin(_theta * pi/180.0); 262 263 // Compute the angle btw vector q and the 264 // axis of the cylinder (assume qz = 0) 265 const double cos_val = cyl_x*qx + cyl_y*qy; 266 267 // The following test should always pass 268 if (fabs(cos_val)>1.0) { 269 return 0; 270 } 271 272 // Note: cos(alpha) = 0 and 1 will get an 273 // undefined value from CylKernel 274 const double alpha = acos( cos_val ); 275 276 // Call the IGOR library function to get the kernel 277 const double output = bar2d_kernel(dp, q, alpha)/sin(alpha); 278 337 279 double _ptvalue = weights_rad_bar[i].weight 338 280 * weights_len_bar[j].weight … … 341 283 * weights_phi[m].weight 342 284 * vol_i 343 * barbell_analytical_2DXY(&dp, qx, qy); 344 //* output; 285 * output; 345 286 //* pow(weights_rad[i].value,3.0); 346 287
Note: See TracChangeset
for help on using the changeset viewer.