[339ce67] | 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; |
---|
| 26 | |
---|
| 27 | extern "C" { |
---|
| 28 | #include "libCylinder.h" |
---|
[85b856b] | 29 | #include "GaussWeights.h" |
---|
[339ce67] | 30 | #include "barbell.h" |
---|
| 31 | } |
---|
| 32 | |
---|
| 33 | BarBellModel :: BarBellModel() { |
---|
| 34 | scale = Parameter(1.0); |
---|
| 35 | rad_bar = Parameter(20.0); |
---|
| 36 | rad_bar.set_min(0.0); |
---|
| 37 | len_bar = Parameter(400.0, true); |
---|
| 38 | len_bar.set_min(0.0); |
---|
| 39 | rad_bell = Parameter(40.0); |
---|
| 40 | rad_bell.set_min(0.0); |
---|
| 41 | sld_barbell = Parameter(1.0e-6); |
---|
| 42 | sld_solv = Parameter(6.3e-6); |
---|
| 43 | background = Parameter(0.0); |
---|
| 44 | theta = Parameter(0.0, true); |
---|
| 45 | phi = Parameter(0.0, true); |
---|
| 46 | } |
---|
| 47 | |
---|
[85b856b] | 48 | double bar2d_kernel(double dp[], double q, double alpha) { |
---|
| 49 | int j; |
---|
| 50 | double Pi; |
---|
| 51 | double scale,contr,bkg,sldc,slds; |
---|
| 52 | double len,rad,hDist,endRad; |
---|
| 53 | int nordj=76; |
---|
| 54 | double zi=alpha,yyy,answer; //running tally of integration |
---|
| 55 | double summj,vaj,vbj,zij; //for the inner integration |
---|
| 56 | double arg1,arg2,inner,be; |
---|
| 57 | |
---|
| 58 | |
---|
| 59 | scale = dp[0]; |
---|
| 60 | rad = dp[1]; |
---|
| 61 | len = dp[2]; |
---|
| 62 | endRad = dp[3]; |
---|
| 63 | sldc = dp[4]; |
---|
| 64 | slds = dp[5]; |
---|
| 65 | bkg = dp[6]; |
---|
| 66 | |
---|
| 67 | hDist = sqrt(fabs(endRad*endRad-rad*rad)); //by definition for this model |
---|
| 68 | |
---|
| 69 | contr = sldc-slds; |
---|
| 70 | |
---|
| 71 | Pi = 4.0*atan(1.0); |
---|
| 72 | vaj = -1.0*hDist/endRad; |
---|
| 73 | vbj = 1.0; //endpoints of inner integral |
---|
| 74 | |
---|
| 75 | summj=0.0; |
---|
| 76 | |
---|
| 77 | for(j=0;j<nordj;j++) { |
---|
| 78 | //20 gauss points for the inner integral |
---|
| 79 | zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; //the "t" dummy |
---|
| 80 | yyy = Gauss76Wt[j] * Dumb_kernel(dp,q,zij,zi); //uses the same Kernel as the Dumbbell, here L>0 |
---|
| 81 | summj += yyy; |
---|
| 82 | } |
---|
| 83 | //now calculate the value of the inner integral |
---|
| 84 | inner = (vbj-vaj)/2.0*summj; |
---|
| 85 | inner *= 4.0*Pi*endRad*endRad*endRad; |
---|
| 86 | |
---|
| 87 | //now calculate outer integrand |
---|
| 88 | arg1 = q*len/2.0*cos(zi); |
---|
| 89 | arg2 = q*rad*sin(zi); |
---|
| 90 | yyy = inner; |
---|
| 91 | |
---|
| 92 | if(arg2 == 0) { |
---|
| 93 | be = 0.5; |
---|
| 94 | } else { |
---|
| 95 | be = NR_BessJ1(arg2)/arg2; |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | if(arg1 == 0.0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h |
---|
| 99 | yyy += Pi*rad*rad*len*2.0*be; |
---|
| 100 | } else { |
---|
| 101 | yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; |
---|
| 102 | } |
---|
| 103 | yyy *= yyy; //sin(zi); |
---|
| 104 | answer = yyy; |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0); //divide by volume |
---|
| 108 | answer *= 1.0e8; //convert to cm^-1 |
---|
| 109 | answer *= contr*contr; |
---|
| 110 | answer *= scale; |
---|
| 111 | answer += bkg; |
---|
| 112 | |
---|
| 113 | return answer; |
---|
| 114 | } |
---|
[339ce67] | 115 | /** |
---|
| 116 | * Function to evaluate 1D scattering function |
---|
| 117 | * The NIST IGOR library is used for the actual calculation. |
---|
| 118 | * @param q: q-value |
---|
| 119 | * @return: function value |
---|
| 120 | */ |
---|
| 121 | double BarBellModel :: operator()(double q) { |
---|
| 122 | double dp[7]; |
---|
| 123 | |
---|
| 124 | // Fill parameter array for IGOR library |
---|
| 125 | // Add the background after averaging |
---|
| 126 | dp[0] = scale(); |
---|
| 127 | dp[1] = rad_bar(); |
---|
| 128 | dp[2] = len_bar(); |
---|
| 129 | dp[3] = rad_bell(); |
---|
| 130 | dp[4] = sld_barbell(); |
---|
| 131 | dp[5] = sld_solv(); |
---|
| 132 | dp[6] = 0.0; |
---|
| 133 | |
---|
| 134 | // Get the dispersion points for the rad_bar |
---|
| 135 | vector<WeightPoint> weights_rad_bar; |
---|
| 136 | rad_bar.get_weights(weights_rad_bar); |
---|
| 137 | // Get the dispersion points for the len_bar |
---|
| 138 | vector<WeightPoint> weights_len_bar; |
---|
| 139 | len_bar.get_weights(weights_len_bar); |
---|
| 140 | // Get the dispersion points for the rad_bell |
---|
| 141 | vector<WeightPoint> weights_rad_bell; |
---|
| 142 | rad_bell.get_weights(weights_rad_bell); |
---|
| 143 | |
---|
| 144 | // Perform the computation, with all weight points |
---|
| 145 | double sum = 0.0; |
---|
| 146 | double norm = 0.0; |
---|
| 147 | double vol = 0.0; |
---|
| 148 | double pi,hDist,result; |
---|
| 149 | double vol_i = 0.0; |
---|
| 150 | pi = 4.0*atan(1.0); |
---|
| 151 | // Loop over radius weight points |
---|
[34c2649] | 152 | for(size_t i=0; i<weights_rad_bar.size(); i++) { |
---|
[339ce67] | 153 | dp[1] = weights_rad_bar[i].value; |
---|
[34c2649] | 154 | for(size_t j=0; j<weights_len_bar.size(); j++) { |
---|
[339ce67] | 155 | dp[2] = weights_len_bar[j].value; |
---|
[34c2649] | 156 | for(size_t k=0; k<weights_rad_bell.size(); k++) { |
---|
[339ce67] | 157 | dp[3] = weights_rad_bell[k].value; |
---|
| 158 | |
---|
| 159 | //Un-normalize SphereForm by volume |
---|
| 160 | hDist = sqrt(fabs(dp[3]*dp[3]-dp[1]*dp[1])); |
---|
| 161 | vol_i = pi*dp[1]*dp[1]*dp[2]+2.0*pi*(2.0*dp[3]*dp[3]*dp[3]/3.0 |
---|
| 162 | +dp[3]*dp[3]*hDist-hDist*hDist*hDist/3.0); |
---|
| 163 | result = Barbell(dp, q) * vol_i; |
---|
| 164 | // This FIXES a singualrity the kernel in libigor. |
---|
| 165 | if ( result == INFINITY || result == NAN){ |
---|
| 166 | result = 0.0; |
---|
| 167 | } |
---|
| 168 | sum += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight |
---|
| 169 | * result; |
---|
| 170 | //Find average volume |
---|
| 171 | vol += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight |
---|
| 172 | * vol_i; |
---|
| 173 | |
---|
| 174 | norm += weights_rad_bar[i].weight*weights_len_bar[j].weight*weights_rad_bell[k].weight; |
---|
| 175 | } |
---|
| 176 | } |
---|
| 177 | } |
---|
| 178 | |
---|
| 179 | if (vol != 0.0 && norm != 0.0) { |
---|
| 180 | //Re-normalize by avg volume |
---|
| 181 | sum = sum/(vol/norm);} |
---|
| 182 | return sum/norm + background(); |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | /** |
---|
| 186 | * Function to evaluate 2D scattering function |
---|
| 187 | * @param q_x: value of Q along x |
---|
| 188 | * @param q_y: value of Q along y |
---|
| 189 | * @return: function value |
---|
| 190 | */ |
---|
| 191 | double BarBellModel :: operator()(double qx, double qy) { |
---|
[85b856b] | 192 | double dp[7]; |
---|
[339ce67] | 193 | |
---|
[85b856b] | 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; |
---|
[339ce67] | 203 | |
---|
[ef93acb] | 204 | double _theta = theta(); |
---|
| 205 | double _phi = phi(); |
---|
[339ce67] | 206 | |
---|
| 207 | // Get the dispersion points for the rad_bar |
---|
| 208 | vector<WeightPoint> weights_rad_bar; |
---|
| 209 | rad_bar.get_weights(weights_rad_bar); |
---|
| 210 | |
---|
| 211 | // Get the dispersion points for the len_bar |
---|
| 212 | vector<WeightPoint> weights_len_bar; |
---|
| 213 | len_bar.get_weights(weights_len_bar); |
---|
| 214 | |
---|
| 215 | // Get the dispersion points for the rad_bell |
---|
| 216 | vector<WeightPoint> weights_rad_bell; |
---|
| 217 | rad_bell.get_weights(weights_rad_bell); |
---|
| 218 | |
---|
| 219 | // Get angular averaging for theta |
---|
| 220 | vector<WeightPoint> weights_theta; |
---|
| 221 | theta.get_weights(weights_theta); |
---|
| 222 | |
---|
| 223 | // Get angular averaging for phi |
---|
| 224 | vector<WeightPoint> weights_phi; |
---|
| 225 | phi.get_weights(weights_phi); |
---|
| 226 | |
---|
| 227 | |
---|
| 228 | // Perform the computation, with all weight points |
---|
| 229 | double sum = 0.0; |
---|
| 230 | double norm = 0.0; |
---|
| 231 | double norm_vol = 0.0; |
---|
| 232 | double vol = 0.0; |
---|
[34c2649] | 233 | double pi,hDist; |
---|
[339ce67] | 234 | double vol_i = 0.0; |
---|
| 235 | pi = 4.0*atan(1.0); |
---|
| 236 | |
---|
| 237 | // Loop over radius weight points |
---|
[34c2649] | 238 | for(size_t i=0; i<weights_rad_bar.size(); i++) { |
---|
[85b856b] | 239 | dp[1] = weights_rad_bar[i].value; |
---|
[34c2649] | 240 | for(size_t j=0; j<weights_len_bar.size(); j++) { |
---|
[85b856b] | 241 | dp[2] = weights_len_bar[j].value; |
---|
[34c2649] | 242 | for(size_t k=0; k<weights_rad_bell.size(); k++) { |
---|
[85b856b] | 243 | dp[3] = weights_rad_bell[k].value; |
---|
[339ce67] | 244 | // Average over theta distribution |
---|
[34c2649] | 245 | for(size_t l=0; l< weights_theta.size(); l++) { |
---|
[85b856b] | 246 | _theta = weights_theta[l].value; |
---|
[339ce67] | 247 | // Average over phi distribution |
---|
[34c2649] | 248 | for(size_t m=0; m< weights_phi.size(); m++) { |
---|
[85b856b] | 249 | _phi = weights_phi[m].value; |
---|
[339ce67] | 250 | //Un-normalize Form by volume |
---|
[85b856b] | 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 |
---|
[318b5bbb] | 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); |
---|
[85b856b] | 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); |
---|
[339ce67] | 278 | |
---|
| 279 | double _ptvalue = weights_rad_bar[i].weight |
---|
| 280 | * weights_len_bar[j].weight |
---|
| 281 | * weights_rad_bell[k].weight |
---|
| 282 | * weights_theta[l].weight |
---|
| 283 | * weights_phi[m].weight |
---|
| 284 | * vol_i |
---|
[85b856b] | 285 | * output; |
---|
[339ce67] | 286 | //* pow(weights_rad[i].value,3.0); |
---|
[85b856b] | 287 | |
---|
[339ce67] | 288 | // Consider when there is infinte or nan. |
---|
| 289 | if ( _ptvalue == INFINITY || _ptvalue == NAN){ |
---|
| 290 | _ptvalue = 0.0; |
---|
| 291 | } |
---|
| 292 | if (weights_theta.size()>1) { |
---|
[318b5bbb] | 293 | _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0)); |
---|
[339ce67] | 294 | } |
---|
| 295 | sum += _ptvalue; |
---|
| 296 | // This model dose not need the volume of spheres correction!!! |
---|
| 297 | //Find average volume |
---|
| 298 | vol += weights_rad_bar[i].weight |
---|
| 299 | * weights_len_bar[j].weight |
---|
| 300 | * weights_rad_bell[k].weight |
---|
| 301 | * vol_i; |
---|
| 302 | //Find norm for volume |
---|
| 303 | norm_vol += weights_rad_bar[i].weight |
---|
| 304 | * weights_len_bar[j].weight |
---|
| 305 | * weights_rad_bell[k].weight; |
---|
| 306 | |
---|
| 307 | norm += weights_rad_bar[i].weight |
---|
| 308 | * weights_len_bar[j].weight |
---|
| 309 | * weights_rad_bell[k].weight |
---|
| 310 | * weights_theta[l].weight |
---|
| 311 | * weights_phi[m].weight; |
---|
| 312 | } |
---|
| 313 | } |
---|
| 314 | } |
---|
| 315 | } |
---|
| 316 | } |
---|
| 317 | // Averaging in theta needs an extra normalization |
---|
| 318 | // factor to account for the sin(theta) term in the |
---|
| 319 | // integration (see documentation). |
---|
| 320 | if (weights_theta.size()>1) norm = norm / asin(1.0); |
---|
| 321 | |
---|
| 322 | if (vol != 0.0 && norm_vol != 0.0) { |
---|
| 323 | //Re-normalize by avg volume |
---|
| 324 | sum = sum/(vol/norm_vol);} |
---|
| 325 | |
---|
| 326 | return sum/norm + background(); |
---|
| 327 | } |
---|
| 328 | |
---|
| 329 | /** |
---|
| 330 | * Function to evaluate 2D scattering function |
---|
| 331 | * @param pars: parameters of the SCCrystal |
---|
| 332 | * @param q: q-value |
---|
| 333 | * @param phi: angle phi |
---|
| 334 | * @return: function value |
---|
| 335 | */ |
---|
| 336 | double BarBellModel :: evaluate_rphi(double q, double phi) { |
---|
| 337 | return (*this).operator()(q); |
---|
| 338 | } |
---|
| 339 | |
---|
| 340 | /** |
---|
| 341 | * Function to calculate effective radius |
---|
| 342 | * @return: effective radius value |
---|
| 343 | */ |
---|
| 344 | double BarBellModel :: calculate_ER() { |
---|
| 345 | //NOT implemented yet!!! |
---|
[34c2649] | 346 | return 0.0; |
---|
[339ce67] | 347 | } |
---|
[e08bd5b] | 348 | /** |
---|
| 349 | * Function to calculate volf_ratio for shell/tot |
---|
| 350 | * @return: volf_ratio value |
---|
| 351 | */ |
---|
| 352 | double BarBellModel :: calculate_VR() { |
---|
| 353 | return 1.0; |
---|
| 354 | } |
---|