[c724ccd] | 1 | /** |
---|
| 2 | This software was developed by the University of Tennessee as part of the |
---|
| 3 | Distributed Data Analysis of Neutron Scattering Experiments (DANSE) |
---|
| 4 | project funded by the US National Science Foundation. |
---|
| 5 | |
---|
| 6 | If you use DANSE applications to do scientific research that leads to |
---|
| 7 | publication, we ask that you acknowledge the use of the software with the |
---|
| 8 | following sentence: |
---|
| 9 | |
---|
| 10 | "This work benefited from DANSE software developed under NSF award DMR-0520547." |
---|
| 11 | |
---|
| 12 | copyright 2008, University of Tennessee |
---|
| 13 | */ |
---|
| 14 | |
---|
| 15 | /** |
---|
| 16 | * Scattering model classes |
---|
| 17 | * The classes use the IGOR library found in |
---|
| 18 | * sansmodels/src/libigor |
---|
| 19 | * |
---|
| 20 | */ |
---|
| 21 | |
---|
| 22 | #include <math.h> |
---|
| 23 | #include "parameters.hh" |
---|
| 24 | #include <stdio.h> |
---|
| 25 | using namespace std; |
---|
[011e0e4] | 26 | #include "corefourshell.h" |
---|
[c724ccd] | 27 | |
---|
| 28 | extern "C" { |
---|
| 29 | #include "libSphere.h" |
---|
[318b5bbb] | 30 | #include "libmultifunc/libfunc.h" |
---|
[c724ccd] | 31 | } |
---|
| 32 | |
---|
[011e0e4] | 33 | typedef struct { |
---|
| 34 | double scale; |
---|
| 35 | double rad_core0; |
---|
| 36 | double sld_core0; |
---|
| 37 | double thick_shell1; |
---|
| 38 | double sld_shell1; |
---|
| 39 | double thick_shell2; |
---|
| 40 | double sld_shell2; |
---|
| 41 | double thick_shell3; |
---|
| 42 | double sld_shell3; |
---|
| 43 | double thick_shell4; |
---|
| 44 | double sld_shell4; |
---|
| 45 | double sld_solv; |
---|
| 46 | double background; |
---|
[318b5bbb] | 47 | double M0_sld_shell1; |
---|
| 48 | double M_theta_shell1; |
---|
| 49 | double M_phi_shell1; |
---|
| 50 | double M0_sld_shell2; |
---|
| 51 | double M_theta_shell2; |
---|
| 52 | double M_phi_shell2; |
---|
| 53 | double M0_sld_shell3; |
---|
| 54 | double M_theta_shell3; |
---|
| 55 | double M_phi_shell3; |
---|
| 56 | double M0_sld_shell4; |
---|
| 57 | double M_theta_shell4; |
---|
| 58 | double M_phi_shell4; |
---|
| 59 | double M0_sld_core0; |
---|
| 60 | double M_theta_core0; |
---|
| 61 | double M_phi_core0; |
---|
| 62 | double M0_sld_solv; |
---|
| 63 | double M_theta_solv; |
---|
| 64 | double M_phi_solv; |
---|
| 65 | double Up_frac_i; |
---|
| 66 | double Up_frac_f; |
---|
| 67 | double Up_theta; |
---|
[011e0e4] | 68 | } CoreFourShellParameters; |
---|
| 69 | |
---|
[c724ccd] | 70 | CoreFourShellModel :: CoreFourShellModel() { |
---|
| 71 | scale = Parameter(1.0); |
---|
[339ce67] | 72 | rad_core0 = Parameter(60.0, true); |
---|
| 73 | rad_core0.set_min(0.0); |
---|
| 74 | sld_core0 = Parameter(6.4e-6); |
---|
[c724ccd] | 75 | thick_shell1 = Parameter(10.0, true); |
---|
| 76 | thick_shell1.set_min(0.0); |
---|
| 77 | sld_shell1 = Parameter(1.0e-6); |
---|
| 78 | thick_shell2 = Parameter(10.0, true); |
---|
| 79 | thick_shell2.set_min(0.0); |
---|
| 80 | sld_shell2 = Parameter(2.0e-6); |
---|
| 81 | thick_shell3 = Parameter(10.0, true); |
---|
| 82 | thick_shell3.set_min(0.0); |
---|
| 83 | sld_shell3 = Parameter(3.0e-6); |
---|
| 84 | thick_shell4 = Parameter(10.0, true); |
---|
| 85 | thick_shell4.set_min(0.0); |
---|
| 86 | sld_shell4 = Parameter(4.0e-6); |
---|
| 87 | sld_solv = Parameter(6.4e-6); |
---|
| 88 | background = Parameter(0.001); |
---|
[318b5bbb] | 89 | M0_sld_shell1 = Parameter(0.0e-6); |
---|
| 90 | M_theta_shell1 = Parameter(0.0); |
---|
| 91 | M_phi_shell1 = Parameter(0.0); |
---|
| 92 | M0_sld_shell2 = Parameter(0.0e-6); |
---|
| 93 | M_theta_shell2 = Parameter(0.0); |
---|
| 94 | M_phi_shell2 = Parameter(0.0); |
---|
| 95 | M0_sld_shell3 = Parameter(0.0e-6); |
---|
| 96 | M_theta_shell3 = Parameter(0.0); |
---|
| 97 | M_phi_shell3 = Parameter(0.0); |
---|
| 98 | M0_sld_shell4 = Parameter(0.0e-6); |
---|
| 99 | M_theta_shell4 = Parameter(0.0); |
---|
| 100 | M_phi_shell4 = Parameter(0.0); |
---|
| 101 | M0_sld_core0 = Parameter(0.0e-6); |
---|
| 102 | M_theta_core0 = Parameter(0.0); |
---|
| 103 | M_phi_core0 = Parameter(0.0); |
---|
| 104 | M0_sld_solv = Parameter(0.0e-6); |
---|
| 105 | M_theta_solv = Parameter(0.0); |
---|
| 106 | M_phi_solv = Parameter(0.0); |
---|
| 107 | Up_frac_i = Parameter(0.5); |
---|
| 108 | Up_frac_f = Parameter(0.5); |
---|
| 109 | Up_theta = Parameter(0.0); |
---|
| 110 | |
---|
[c724ccd] | 111 | } |
---|
| 112 | |
---|
| 113 | /** |
---|
| 114 | * Function to evaluate 1D scattering function |
---|
| 115 | * The NIST IGOR library is used for the actual calculation. |
---|
| 116 | * @param q: q-value |
---|
| 117 | * @return: function value |
---|
| 118 | */ |
---|
| 119 | double CoreFourShellModel :: operator()(double q) { |
---|
| 120 | double dp[13]; |
---|
| 121 | |
---|
| 122 | // Fill parameter array for IGOR library |
---|
| 123 | // Add the background after averaging |
---|
[318b5bbb] | 124 | |
---|
[c724ccd] | 125 | dp[0] = scale(); |
---|
[339ce67] | 126 | dp[1] = rad_core0(); |
---|
| 127 | dp[2] = sld_core0(); |
---|
[c724ccd] | 128 | dp[3] = thick_shell1(); |
---|
| 129 | dp[4] = sld_shell1(); |
---|
| 130 | dp[5] = thick_shell2(); |
---|
| 131 | dp[6] = sld_shell2(); |
---|
| 132 | dp[7] = thick_shell3(); |
---|
| 133 | dp[8] = sld_shell3(); |
---|
| 134 | dp[9] = thick_shell4(); |
---|
| 135 | dp[10] = sld_shell4(); |
---|
| 136 | dp[11] = sld_solv(); |
---|
| 137 | dp[12] = 0.0; |
---|
| 138 | |
---|
| 139 | // Get the dispersion points for the radius |
---|
| 140 | vector<WeightPoint> weights_rad; |
---|
[339ce67] | 141 | rad_core0.get_weights(weights_rad); |
---|
[c724ccd] | 142 | |
---|
| 143 | // Get the dispersion points for the thick 1 |
---|
| 144 | vector<WeightPoint> weights_s1; |
---|
| 145 | thick_shell1.get_weights(weights_s1); |
---|
| 146 | |
---|
| 147 | // Get the dispersion points for the thick 2 |
---|
| 148 | vector<WeightPoint> weights_s2; |
---|
| 149 | thick_shell2.get_weights(weights_s2); |
---|
| 150 | |
---|
| 151 | // Get the dispersion points for the thick 3 |
---|
| 152 | vector<WeightPoint> weights_s3; |
---|
| 153 | thick_shell3.get_weights(weights_s3); |
---|
| 154 | |
---|
| 155 | // Get the dispersion points for the thick 4 |
---|
| 156 | vector<WeightPoint> weights_s4; |
---|
| 157 | thick_shell4.get_weights(weights_s4); |
---|
| 158 | |
---|
| 159 | // Perform the computation, with all weight points |
---|
| 160 | double sum = 0.0; |
---|
| 161 | double norm = 0.0; |
---|
| 162 | double vol = 0.0; |
---|
| 163 | |
---|
| 164 | // Loop over radius weight points |
---|
[34c2649] | 165 | for(size_t i=0; i<weights_rad.size(); i++) { |
---|
[c724ccd] | 166 | dp[1] = weights_rad[i].value; |
---|
| 167 | // Loop over radius weight points |
---|
[34c2649] | 168 | for(size_t j=0; j<weights_s1.size(); j++) { |
---|
[c724ccd] | 169 | dp[3] = weights_s1[j].value; |
---|
| 170 | // Loop over radius weight points |
---|
[34c2649] | 171 | for(size_t k=0; k<weights_s2.size(); k++) { |
---|
[c724ccd] | 172 | dp[5] = weights_s2[k].value; |
---|
| 173 | // Loop over radius weight points |
---|
[34c2649] | 174 | for(size_t l=0; l<weights_s3.size(); l++) { |
---|
[c724ccd] | 175 | dp[7] = weights_s3[l].value; |
---|
| 176 | // Loop over radius weight points |
---|
[34c2649] | 177 | for(size_t m=0; m<weights_s4.size(); m++) { |
---|
[c724ccd] | 178 | dp[9] = weights_s4[m].value; |
---|
| 179 | //Un-normalize FourShell by volume |
---|
| 180 | sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight |
---|
| 181 | * FourShell(dp, q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3); |
---|
| 182 | //Find average volume |
---|
| 183 | vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight |
---|
| 184 | * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3); |
---|
| 185 | |
---|
| 186 | norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight; |
---|
| 187 | } |
---|
| 188 | } |
---|
| 189 | } |
---|
| 190 | } |
---|
| 191 | } |
---|
| 192 | |
---|
| 193 | if (vol != 0.0 && norm != 0.0) { |
---|
| 194 | //Re-normalize by avg volume |
---|
| 195 | sum = sum/(vol/norm);} |
---|
| 196 | return sum/norm + background(); |
---|
| 197 | } |
---|
| 198 | |
---|
| 199 | /** |
---|
| 200 | * Function to evaluate 2D scattering function |
---|
| 201 | * @param q_x: value of Q along x |
---|
| 202 | * @param q_y: value of Q along y |
---|
| 203 | * @return: function value |
---|
| 204 | */ |
---|
[318b5bbb] | 205 | |
---|
| 206 | static double corefourshell_analytical_2D_scaled(CoreFourShellParameters *pars, double q, double q_x, double q_y) { |
---|
| 207 | double dp[13]; |
---|
| 208 | |
---|
| 209 | // Fill parameter array for IGOR library |
---|
| 210 | // Add the background after averaging |
---|
| 211 | |
---|
| 212 | dp[0] = pars->scale; |
---|
| 213 | dp[1] = pars->rad_core0; |
---|
| 214 | dp[2] = 0.0; //sld_core0; |
---|
| 215 | dp[3] = pars->thick_shell1; |
---|
| 216 | dp[4] = 0.0; //sld_shell1; |
---|
| 217 | dp[5] = pars->thick_shell2; |
---|
| 218 | dp[6] = 0.0; //sld_shell2; |
---|
| 219 | dp[7] = pars->thick_shell3; |
---|
| 220 | dp[8] = 0.0; //sld_shell3; |
---|
| 221 | dp[9] = pars->thick_shell4; |
---|
| 222 | dp[10] = 0.0; //sld_shell4; |
---|
| 223 | dp[11] = 0.0; //sld_solv; |
---|
| 224 | dp[12] = 0.0; |
---|
| 225 | |
---|
| 226 | double sld_core0 = pars->sld_core0; |
---|
| 227 | double sld_shell1 = pars->sld_shell1; |
---|
| 228 | double sld_shell2 = pars->sld_shell2; |
---|
| 229 | double sld_shell3 = pars->sld_shell3; |
---|
| 230 | double sld_shell4 = pars->sld_shell4; |
---|
| 231 | double sld_solv = pars->sld_solv; |
---|
| 232 | double answer = 0.0; |
---|
| 233 | double m_max0 = pars->M0_sld_core0; |
---|
| 234 | double m_max_shell1 = pars->M0_sld_shell1; |
---|
| 235 | double m_max_shell2 = pars->M0_sld_shell2; |
---|
| 236 | double m_max_shell3 = pars->M0_sld_shell3; |
---|
| 237 | double m_max_shell4 = pars->M0_sld_shell4; |
---|
| 238 | double m_max_solv = pars->M0_sld_solv; |
---|
| 239 | |
---|
| 240 | if (m_max0 < 1.0e-32 && m_max_solv < 1.0e-32 && m_max_shell1 < 1.0e-32 && m_max_shell2 < |
---|
| 241 | 1.0e-32 && m_max_shell3 < 1.0e-32 && m_max_shell4 < 1.0e-32){ |
---|
| 242 | dp[2] = sld_core0; |
---|
| 243 | dp[4] = sld_shell1; |
---|
| 244 | dp[6] = sld_shell2; |
---|
| 245 | dp[8] = sld_shell3; |
---|
| 246 | dp[10] = sld_shell4; |
---|
| 247 | dp[11] = sld_solv; |
---|
| 248 | answer = FourShell(dp, q); |
---|
| 249 | } |
---|
| 250 | else{ |
---|
| 251 | double qx = q_x; |
---|
| 252 | double qy = q_y; |
---|
| 253 | double s_theta = pars->Up_theta; |
---|
| 254 | double m_phi0 = pars->M_phi_core0; |
---|
| 255 | double m_theta0 = pars->M_theta_core0; |
---|
| 256 | double m_phi_shell1 = pars->M_phi_shell1; |
---|
| 257 | double m_theta_shell1 = pars->M_theta_shell1; |
---|
| 258 | double m_phi_shell2 = pars->M_phi_shell2; |
---|
| 259 | double m_theta_shell2 = pars->M_theta_shell2; |
---|
| 260 | double m_phi_shell3 = pars->M_phi_shell3; |
---|
| 261 | double m_theta_shell3 = pars->M_theta_shell3; |
---|
| 262 | double m_phi_shell4 = pars->M_phi_shell4; |
---|
| 263 | double m_theta_shell4 = pars->M_theta_shell4; |
---|
| 264 | double m_phi_solv = pars->M_phi_solv; |
---|
| 265 | double m_theta_solv = pars->M_theta_solv; |
---|
| 266 | double in_spin = pars->Up_frac_i; |
---|
| 267 | double out_spin = pars->Up_frac_f; |
---|
| 268 | polar_sld p_sld_core0; |
---|
| 269 | polar_sld p_sld_shell1; |
---|
| 270 | polar_sld p_sld_shell2; |
---|
| 271 | polar_sld p_sld_shell3; |
---|
| 272 | polar_sld p_sld_shell4; |
---|
| 273 | polar_sld p_sld_solv; |
---|
| 274 | //Find (b+m) slds |
---|
| 275 | p_sld_core0 = cal_msld(1, qx, qy, sld_core0, m_max0, m_theta0, m_phi0, |
---|
| 276 | in_spin, out_spin, s_theta); |
---|
| 277 | p_sld_shell1 = cal_msld(1, qx, qy, sld_shell1, m_max_shell1, m_theta_shell1, m_phi_shell1, |
---|
| 278 | in_spin, out_spin, s_theta); |
---|
| 279 | p_sld_shell2 = cal_msld(1, qx, qy, sld_shell2, m_max_shell2, m_theta_shell2, m_phi_shell2, |
---|
| 280 | in_spin, out_spin, s_theta); |
---|
| 281 | p_sld_shell3 = cal_msld(1, qx, qy, sld_shell3, m_max_shell3, m_theta_shell3, m_phi_shell3, |
---|
| 282 | in_spin, out_spin, s_theta); |
---|
| 283 | p_sld_shell4 = cal_msld(1, qx, qy, sld_shell4, m_max_shell4, m_theta_shell4, m_phi_shell4, |
---|
| 284 | in_spin, out_spin, s_theta); |
---|
| 285 | p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv, |
---|
| 286 | in_spin, out_spin, s_theta); |
---|
| 287 | //up_up |
---|
| 288 | if (in_spin > 0.0 && out_spin > 0.0){ |
---|
| 289 | dp[2] = p_sld_core0.uu; |
---|
| 290 | dp[4] = p_sld_shell1.uu; |
---|
| 291 | dp[6] = p_sld_shell2.uu; |
---|
| 292 | dp[8] = p_sld_shell3.uu; |
---|
| 293 | dp[10] = p_sld_shell4.uu; |
---|
| 294 | dp[11] = p_sld_solv.uu; |
---|
| 295 | answer += FourShell(dp, q); |
---|
| 296 | } |
---|
| 297 | //down_down |
---|
| 298 | if (in_spin < 1.0 && out_spin < 1.0){ |
---|
| 299 | dp[2] = p_sld_core0.dd; |
---|
| 300 | dp[4] = p_sld_shell1.dd; |
---|
| 301 | dp[6] = p_sld_shell2.dd; |
---|
| 302 | dp[8] = p_sld_shell3.dd; |
---|
| 303 | dp[10] = p_sld_shell4.dd; |
---|
| 304 | dp[11] = p_sld_solv.dd; |
---|
| 305 | answer += FourShell(dp, q); |
---|
| 306 | } |
---|
| 307 | //up_down |
---|
| 308 | if (in_spin > 0.0 && out_spin < 1.0){ |
---|
| 309 | dp[2] = p_sld_core0.re_ud; |
---|
| 310 | dp[4] = p_sld_shell1.re_ud; |
---|
| 311 | dp[6] = p_sld_shell2.re_ud; |
---|
| 312 | dp[8] = p_sld_shell3.re_ud; |
---|
| 313 | dp[10] = p_sld_shell4.re_ud; |
---|
| 314 | dp[11] = p_sld_solv.re_ud; |
---|
| 315 | answer += FourShell(dp, q); |
---|
| 316 | dp[2] = p_sld_core0.im_ud; |
---|
| 317 | dp[4] = p_sld_shell1.im_ud; |
---|
| 318 | dp[6] = p_sld_shell2.im_ud; |
---|
| 319 | dp[8] = p_sld_shell3.im_ud; |
---|
| 320 | dp[10] = p_sld_shell4.im_ud; |
---|
| 321 | dp[11] = p_sld_solv.im_ud; |
---|
| 322 | answer += FourShell(dp, q); |
---|
| 323 | } |
---|
| 324 | //down_up |
---|
| 325 | if (in_spin < 1.0 && out_spin > 0.0){ |
---|
| 326 | dp[2] = p_sld_core0.re_du; |
---|
| 327 | dp[4] = p_sld_shell1.re_du; |
---|
| 328 | dp[6] = p_sld_shell2.re_du; |
---|
| 329 | dp[8] = p_sld_shell3.re_du; |
---|
| 330 | dp[10] = p_sld_shell4.re_du; |
---|
| 331 | dp[11] = p_sld_solv.re_du; |
---|
| 332 | answer += FourShell(dp, q); |
---|
| 333 | dp[2] = p_sld_core0.im_du; |
---|
| 334 | dp[4] = p_sld_shell1.im_du; |
---|
| 335 | dp[6] = p_sld_shell2.im_du; |
---|
| 336 | dp[8] = p_sld_shell3.im_du; |
---|
| 337 | dp[10] = p_sld_shell4.im_du; |
---|
| 338 | dp[11] = p_sld_solv.im_du; |
---|
| 339 | answer += FourShell(dp, q); |
---|
| 340 | } |
---|
| 341 | } |
---|
| 342 | // Already normalized |
---|
| 343 | // add in the background |
---|
| 344 | answer += pars->background; |
---|
| 345 | return answer; |
---|
| 346 | } |
---|
| 347 | |
---|
| 348 | /** |
---|
| 349 | * Function to evaluate 2D scattering function |
---|
| 350 | * @param pars: parameters |
---|
| 351 | * @param q: q-value |
---|
| 352 | * @return: function value |
---|
| 353 | */ |
---|
| 354 | static double corefourshell_analytical_2DXY(CoreFourShellParameters *pars, double qx, double qy) { |
---|
| 355 | double q; |
---|
| 356 | q = sqrt(qx*qx+qy*qy); |
---|
| 357 | return corefourshell_analytical_2D_scaled(pars, q, qx/q, qy/q); |
---|
| 358 | } |
---|
| 359 | |
---|
| 360 | |
---|
[c724ccd] | 361 | double CoreFourShellModel :: operator()(double qx, double qy) { |
---|
[318b5bbb] | 362 | CoreFourShellParameters dp; |
---|
| 363 | dp.scale = scale(); |
---|
| 364 | dp.rad_core0 = rad_core0(); |
---|
| 365 | dp.sld_core0 = sld_core0(); |
---|
| 366 | dp.thick_shell1 = thick_shell1(); |
---|
| 367 | dp.sld_shell1 = sld_shell1(); |
---|
| 368 | dp.thick_shell2 = thick_shell2(); |
---|
| 369 | dp.sld_shell2 = sld_shell2(); |
---|
| 370 | dp.thick_shell3 = thick_shell3(); |
---|
| 371 | dp.sld_shell3 = sld_shell3(); |
---|
| 372 | dp.thick_shell4 = thick_shell4(); |
---|
| 373 | dp.sld_shell4 = sld_shell4(); |
---|
| 374 | dp.sld_solv = sld_solv(); |
---|
| 375 | dp.background = 0.0; |
---|
| 376 | dp.M0_sld_shell1 = M0_sld_shell1(); |
---|
| 377 | dp.M_theta_shell1 = M_theta_shell1(); |
---|
| 378 | dp.M_phi_shell1 = M_phi_shell1(); |
---|
| 379 | dp.M0_sld_shell2 = M0_sld_shell2(); |
---|
| 380 | dp.M_theta_shell2 = M_theta_shell2(); |
---|
| 381 | dp.M_phi_shell2 = M_phi_shell2(); |
---|
| 382 | dp.M0_sld_shell3 = M0_sld_shell3(); |
---|
| 383 | dp.M_theta_shell3 = M_theta_shell3(); |
---|
| 384 | dp.M_phi_shell3 = M_phi_shell3(); |
---|
| 385 | dp.M0_sld_shell4 = M0_sld_shell4(); |
---|
| 386 | dp.M_theta_shell4 = M_theta_shell4(); |
---|
| 387 | dp.M_phi_shell4 = M_phi_shell4(); |
---|
| 388 | dp.M0_sld_core0 = M0_sld_core0(); |
---|
| 389 | dp.M_theta_core0 = M_theta_core0(); |
---|
| 390 | dp.M_phi_core0 = M_phi_core0(); |
---|
| 391 | dp.M0_sld_solv = M0_sld_solv(); |
---|
| 392 | dp.M_theta_solv = M_theta_solv(); |
---|
| 393 | dp.M_phi_solv = M_phi_solv(); |
---|
| 394 | dp.Up_frac_i = Up_frac_i(); |
---|
| 395 | dp.Up_frac_f = Up_frac_f(); |
---|
| 396 | dp.Up_theta = Up_theta(); |
---|
| 397 | |
---|
| 398 | // Get the dispersion points for the radius |
---|
| 399 | vector<WeightPoint> weights_rad; |
---|
| 400 | rad_core0.get_weights(weights_rad); |
---|
| 401 | |
---|
| 402 | // Get the dispersion points for the thick 1 |
---|
| 403 | vector<WeightPoint> weights_s1; |
---|
| 404 | thick_shell1.get_weights(weights_s1); |
---|
| 405 | |
---|
| 406 | // Get the dispersion points for the thick 2 |
---|
| 407 | vector<WeightPoint> weights_s2; |
---|
| 408 | thick_shell2.get_weights(weights_s2); |
---|
| 409 | |
---|
| 410 | // Get the dispersion points for the thick 3 |
---|
| 411 | vector<WeightPoint> weights_s3; |
---|
| 412 | thick_shell3.get_weights(weights_s3); |
---|
| 413 | |
---|
| 414 | // Get the dispersion points for the thick 4 |
---|
| 415 | vector<WeightPoint> weights_s4; |
---|
| 416 | thick_shell4.get_weights(weights_s4); |
---|
| 417 | |
---|
| 418 | // Perform the computation, with all weight points |
---|
| 419 | double sum = 0.0; |
---|
| 420 | double norm = 0.0; |
---|
| 421 | double vol = 0.0; |
---|
| 422 | |
---|
| 423 | // Loop over radius weight points |
---|
| 424 | for(size_t i=0; i<weights_rad.size(); i++) { |
---|
| 425 | dp.rad_core0 = weights_rad[i].value; |
---|
| 426 | // Loop over radius weight points |
---|
| 427 | for(size_t j=0; j<weights_s1.size(); j++) { |
---|
| 428 | dp.thick_shell1 = weights_s1[j].value; |
---|
| 429 | // Loop over radius weight points |
---|
| 430 | for(size_t k=0; k<weights_s2.size(); k++) { |
---|
| 431 | dp.thick_shell2 = weights_s2[k].value; |
---|
| 432 | // Loop over radius weight points |
---|
| 433 | for(size_t l=0; l<weights_s3.size(); l++) { |
---|
| 434 | dp.thick_shell3 = weights_s3[l].value; |
---|
| 435 | // Loop over radius weight points |
---|
| 436 | for(size_t m=0; m<weights_s4.size(); m++) { |
---|
| 437 | dp.thick_shell4 = weights_s4[m].value; |
---|
| 438 | //Un-normalize FourShell by volume |
---|
| 439 | sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight |
---|
| 440 | * corefourshell_analytical_2DXY(&dp, qx, qy) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3); |
---|
| 441 | //Find average volume |
---|
| 442 | vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight |
---|
| 443 | * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3); |
---|
| 444 | |
---|
| 445 | norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight; |
---|
| 446 | } |
---|
| 447 | } |
---|
| 448 | } |
---|
| 449 | } |
---|
| 450 | } |
---|
| 451 | |
---|
| 452 | if (vol != 0.0 && norm != 0.0) { |
---|
| 453 | //Re-normalize by avg volume |
---|
| 454 | sum = sum/(vol/norm);} |
---|
| 455 | return sum/norm + background(); |
---|
[c724ccd] | 456 | } |
---|
| 457 | |
---|
| 458 | /** |
---|
| 459 | * Function to evaluate 2D scattering function |
---|
| 460 | * @param pars: parameters of the sphere |
---|
| 461 | * @param q: q-value |
---|
| 462 | * @param phi: angle phi |
---|
| 463 | * @return: function value |
---|
| 464 | */ |
---|
| 465 | double CoreFourShellModel :: evaluate_rphi(double q, double phi) { |
---|
[318b5bbb] | 466 | double qx = q*cos(phi); |
---|
| 467 | double qy = q*sin(phi); |
---|
| 468 | return (*this).operator()(qx, qy); |
---|
[c724ccd] | 469 | } |
---|
| 470 | |
---|
| 471 | /** |
---|
| 472 | * Function to calculate effective radius |
---|
| 473 | * @return: effective radius value |
---|
| 474 | */ |
---|
| 475 | double CoreFourShellModel :: calculate_ER() { |
---|
| 476 | CoreFourShellParameters dp; |
---|
| 477 | |
---|
| 478 | dp.scale = scale(); |
---|
[339ce67] | 479 | dp.rad_core0 = rad_core0(); |
---|
| 480 | dp.sld_core0 = sld_core0(); |
---|
[c724ccd] | 481 | dp.thick_shell1 = thick_shell1(); |
---|
| 482 | dp.sld_shell1 = sld_shell1(); |
---|
| 483 | dp.thick_shell2 = thick_shell2(); |
---|
| 484 | dp.sld_shell2 = sld_shell2(); |
---|
| 485 | dp.thick_shell3 = thick_shell3(); |
---|
| 486 | dp.sld_shell3 = sld_shell3(); |
---|
| 487 | dp.thick_shell4 = thick_shell4(); |
---|
| 488 | dp.sld_shell4 = sld_shell4(); |
---|
| 489 | dp.sld_solv = sld_solv(); |
---|
| 490 | dp.background = 0.0; |
---|
| 491 | |
---|
| 492 | // Get the dispersion points for the radius |
---|
| 493 | vector<WeightPoint> weights_rad; |
---|
[339ce67] | 494 | rad_core0.get_weights(weights_rad); |
---|
[c724ccd] | 495 | |
---|
| 496 | // Get the dispersion points for the thick 1 |
---|
| 497 | vector<WeightPoint> weights_s1; |
---|
| 498 | thick_shell1.get_weights(weights_s1); |
---|
| 499 | |
---|
| 500 | // Get the dispersion points for the thick 2 |
---|
| 501 | vector<WeightPoint> weights_s2; |
---|
| 502 | thick_shell2.get_weights(weights_s2); |
---|
| 503 | |
---|
| 504 | // Get the dispersion points for the thick 3 |
---|
| 505 | vector<WeightPoint> weights_s3; |
---|
| 506 | thick_shell3.get_weights(weights_s3); |
---|
| 507 | |
---|
| 508 | // Get the dispersion points for the thick 4 |
---|
| 509 | vector<WeightPoint> weights_s4; |
---|
| 510 | thick_shell4.get_weights(weights_s4); |
---|
| 511 | |
---|
| 512 | double rad_out = 0.0; |
---|
| 513 | // Perform the computation, with all weight points |
---|
| 514 | double sum = 0.0; |
---|
| 515 | double norm = 0.0; |
---|
| 516 | |
---|
| 517 | // Loop over radius weight points |
---|
[34c2649] | 518 | for(size_t i=0; i<weights_rad.size(); i++) { |
---|
[339ce67] | 519 | dp.rad_core0 = weights_rad[i].value; |
---|
[c724ccd] | 520 | // Loop over radius weight points |
---|
[34c2649] | 521 | for(size_t j=0; j<weights_s1.size(); j++) { |
---|
[c724ccd] | 522 | dp.thick_shell1 = weights_s1[j].value; |
---|
| 523 | // Loop over radius weight points |
---|
[34c2649] | 524 | for(size_t k=0; k<weights_s2.size(); k++) { |
---|
[c724ccd] | 525 | dp.thick_shell2 = weights_s2[k].value; |
---|
| 526 | // Loop over radius weight points |
---|
[34c2649] | 527 | for(size_t l=0; l<weights_s3.size(); l++) { |
---|
[c724ccd] | 528 | dp.thick_shell3 = weights_s3[l].value; |
---|
| 529 | // Loop over radius weight points |
---|
[34c2649] | 530 | for(size_t m=0; m<weights_s4.size(); m++) { |
---|
[c724ccd] | 531 | dp.thick_shell4 = weights_s4[m].value; |
---|
| 532 | //Un-normalize FourShell by volume |
---|
| 533 | sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight |
---|
[339ce67] | 534 | * (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4); |
---|
[c724ccd] | 535 | norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight; |
---|
| 536 | } |
---|
| 537 | } |
---|
| 538 | } |
---|
| 539 | } |
---|
| 540 | } |
---|
| 541 | if (norm != 0){ |
---|
| 542 | //return the averaged value |
---|
| 543 | rad_out = sum/norm;} |
---|
| 544 | else{ |
---|
| 545 | //return normal value |
---|
[339ce67] | 546 | rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4;} |
---|
[c724ccd] | 547 | |
---|
| 548 | return rad_out; |
---|
| 549 | } |
---|
[e08bd5b] | 550 | double CoreFourShellModel :: calculate_VR() { |
---|
| 551 | return 1.0; |
---|
| 552 | } |
---|