| 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 | * TODO: refactor so that we pull in the old sansmodels.c_extensions |
|---|
| 21 | */ |
|---|
| 22 | |
|---|
| 23 | #include <math.h> |
|---|
| 24 | #include "models.hh" |
|---|
| 25 | #include "parameters.hh" |
|---|
| 26 | #include <stdio.h> |
|---|
| 27 | using namespace std; |
|---|
| 28 | |
|---|
| 29 | extern "C" { |
|---|
| 30 | #include "libCylinder.h" |
|---|
| 31 | #include "libStructureFactor.h" |
|---|
| 32 | #include "hollow_cylinder.h" |
|---|
| 33 | } |
|---|
| 34 | |
|---|
| 35 | HollowCylinderModel :: HollowCylinderModel() { |
|---|
| 36 | scale = Parameter(1.0); |
|---|
| 37 | core_radius = Parameter(20.0, true); |
|---|
| 38 | core_radius.set_min(0.0); |
|---|
| 39 | radius = Parameter(30.0, true); |
|---|
| 40 | radius.set_min(0.0); |
|---|
| 41 | length = Parameter(400.0, true); |
|---|
| 42 | length.set_min(0.0); |
|---|
| 43 | sldCyl = Parameter(6.3e-6); |
|---|
| 44 | sldSolv = Parameter(1.0e-6); |
|---|
| 45 | background = Parameter(0.0); |
|---|
| 46 | axis_theta = Parameter(0.0, true); |
|---|
| 47 | axis_phi = Parameter(0.0, true); |
|---|
| 48 | } |
|---|
| 49 | |
|---|
| 50 | /** |
|---|
| 51 | * Function to evaluate 1D scattering function |
|---|
| 52 | * The NIST IGOR library is used for the actual calculation. |
|---|
| 53 | * @param q: q-value |
|---|
| 54 | * @return: function value |
|---|
| 55 | */ |
|---|
| 56 | double HollowCylinderModel :: operator()(double q) { |
|---|
| 57 | double dp[7]; |
|---|
| 58 | |
|---|
| 59 | dp[0] = scale(); |
|---|
| 60 | dp[1] = core_radius(); |
|---|
| 61 | dp[2] = radius(); |
|---|
| 62 | dp[3] = length(); |
|---|
| 63 | dp[4] = sldCyl(); |
|---|
| 64 | dp[5] = sldSolv(); |
|---|
| 65 | dp[6] = 0.0; |
|---|
| 66 | |
|---|
| 67 | // Get the dispersion points for the core radius |
|---|
| 68 | vector<WeightPoint> weights_core_radius; |
|---|
| 69 | core_radius.get_weights(weights_core_radius); |
|---|
| 70 | |
|---|
| 71 | // Get the dispersion points for the shell radius |
|---|
| 72 | vector<WeightPoint> weights_radius; |
|---|
| 73 | radius.get_weights(weights_radius); |
|---|
| 74 | |
|---|
| 75 | // Get the dispersion points for the length |
|---|
| 76 | vector<WeightPoint> weights_length; |
|---|
| 77 | length.get_weights(weights_length); |
|---|
| 78 | |
|---|
| 79 | // Perform the computation, with all weight points |
|---|
| 80 | double sum = 0.0; |
|---|
| 81 | double norm = 0.0; |
|---|
| 82 | double vol = 0.0; |
|---|
| 83 | |
|---|
| 84 | // Loop over core radius weight points |
|---|
| 85 | for(int i=0; i< (int)weights_core_radius.size(); i++) { |
|---|
| 86 | dp[1] = weights_core_radius[i].value; |
|---|
| 87 | |
|---|
| 88 | // Loop over length weight points |
|---|
| 89 | for(int j=0; j< (int)weights_length.size(); j++) { |
|---|
| 90 | dp[3] = weights_length[j].value; |
|---|
| 91 | |
|---|
| 92 | // Loop over shell radius weight points |
|---|
| 93 | for(int k=0; k< (int)weights_radius.size(); k++) { |
|---|
| 94 | dp[2] = weights_radius[k].value; |
|---|
| 95 | //Un-normalize by volume |
|---|
| 96 | sum += weights_core_radius[i].weight |
|---|
| 97 | * weights_length[j].weight |
|---|
| 98 | * weights_radius[k].weight |
|---|
| 99 | * HollowCylinder(dp, q) |
|---|
| 100 | * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2)) |
|---|
| 101 | * weights_length[j].value; |
|---|
| 102 | //Find average volume |
|---|
| 103 | vol += weights_core_radius[i].weight |
|---|
| 104 | * weights_length[j].weight |
|---|
| 105 | * weights_radius[k].weight |
|---|
| 106 | * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2)) |
|---|
| 107 | * weights_length[j].value; |
|---|
| 108 | |
|---|
| 109 | norm += weights_core_radius[i].weight |
|---|
| 110 | * weights_length[j].weight |
|---|
| 111 | * weights_radius[k].weight; |
|---|
| 112 | } |
|---|
| 113 | } |
|---|
| 114 | } |
|---|
| 115 | if (vol != 0.0 && norm != 0.0) { |
|---|
| 116 | //Re-normalize by avg volume |
|---|
| 117 | sum = sum/(vol/norm);} |
|---|
| 118 | |
|---|
| 119 | return sum/norm + background(); |
|---|
| 120 | } |
|---|
| 121 | |
|---|
| 122 | /** |
|---|
| 123 | * Function to evaluate 2D scattering function |
|---|
| 124 | * @param q_x: value of Q along x |
|---|
| 125 | * @param q_y: value of Q along y |
|---|
| 126 | * @return: function value |
|---|
| 127 | */ |
|---|
| 128 | double HollowCylinderModel :: operator()(double qx, double qy) { |
|---|
| 129 | HollowCylinderParameters dp; |
|---|
| 130 | // Fill parameter array |
|---|
| 131 | dp.scale = scale(); |
|---|
| 132 | dp.core_radius = core_radius(); |
|---|
| 133 | dp.radius = radius(); |
|---|
| 134 | dp.length = length(); |
|---|
| 135 | dp.sldCyl = sldCyl(); |
|---|
| 136 | dp.sldSolv = sldSolv(); |
|---|
| 137 | dp.background = 0.0; |
|---|
| 138 | dp.axis_theta = axis_theta(); |
|---|
| 139 | dp.axis_phi = axis_phi(); |
|---|
| 140 | |
|---|
| 141 | // Get the dispersion points for the core radius |
|---|
| 142 | vector<WeightPoint> weights_core_radius; |
|---|
| 143 | core_radius.get_weights(weights_core_radius); |
|---|
| 144 | |
|---|
| 145 | // Get the dispersion points for the shell radius |
|---|
| 146 | vector<WeightPoint> weights_radius; |
|---|
| 147 | radius.get_weights(weights_radius); |
|---|
| 148 | |
|---|
| 149 | // Get the dispersion points for the length |
|---|
| 150 | vector<WeightPoint> weights_length; |
|---|
| 151 | length.get_weights(weights_length); |
|---|
| 152 | |
|---|
| 153 | // Get angular averaging for theta |
|---|
| 154 | vector<WeightPoint> weights_theta; |
|---|
| 155 | axis_theta.get_weights(weights_theta); |
|---|
| 156 | |
|---|
| 157 | // Get angular averaging for phi |
|---|
| 158 | vector<WeightPoint> weights_phi; |
|---|
| 159 | axis_phi.get_weights(weights_phi); |
|---|
| 160 | |
|---|
| 161 | // Perform the computation, with all weight points |
|---|
| 162 | double sum = 0.0; |
|---|
| 163 | double norm = 0.0; |
|---|
| 164 | double norm_vol = 0.0; |
|---|
| 165 | double vol = 0.0; |
|---|
| 166 | double pi = 4.0*atan(1.0); |
|---|
| 167 | // Loop over core radius weight points |
|---|
| 168 | for(int i=0; i<(int)weights_core_radius.size(); i++) { |
|---|
| 169 | dp.core_radius = weights_core_radius[i].value; |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | // Loop over length weight points |
|---|
| 173 | for(int j=0; j<(int)weights_length.size(); j++) { |
|---|
| 174 | dp.length = weights_length[j].value; |
|---|
| 175 | |
|---|
| 176 | // Loop over shell radius weight points |
|---|
| 177 | for(int m=0; m< (int)weights_radius.size(); m++) { |
|---|
| 178 | dp.radius = weights_radius[m].value; |
|---|
| 179 | |
|---|
| 180 | // Average over theta distribution |
|---|
| 181 | for(int k=0; k< (int)weights_theta.size(); k++) { |
|---|
| 182 | dp.axis_theta = weights_theta[k].value; |
|---|
| 183 | |
|---|
| 184 | // Average over phi distribution |
|---|
| 185 | for(int l=0; l< (int)weights_phi.size(); l++) { |
|---|
| 186 | dp.axis_phi = weights_phi[l].value; |
|---|
| 187 | //Un-normalize by volume |
|---|
| 188 | double _ptvalue = weights_core_radius[i].weight |
|---|
| 189 | * weights_length[j].weight |
|---|
| 190 | * weights_radius[m].weight |
|---|
| 191 | * weights_theta[k].weight |
|---|
| 192 | * weights_phi[l].weight |
|---|
| 193 | * hollow_cylinder_analytical_2DXY(&dp, qx, qy) |
|---|
| 194 | / ((pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2)) |
|---|
| 195 | * weights_length[j].value); |
|---|
| 196 | if (weights_theta.size()>1) { |
|---|
| 197 | _ptvalue *= fabs(sin(weights_theta[k].value * pi/180.0)); |
|---|
| 198 | } |
|---|
| 199 | sum += _ptvalue; |
|---|
| 200 | //Find average volume |
|---|
| 201 | vol += weights_core_radius[i].weight |
|---|
| 202 | * weights_length[j].weight |
|---|
| 203 | * weights_radius[m].weight |
|---|
| 204 | * (pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2)) |
|---|
| 205 | * weights_length[j].value; |
|---|
| 206 | //Find norm for volume |
|---|
| 207 | norm_vol += weights_core_radius[i].weight |
|---|
| 208 | * weights_length[j].weight |
|---|
| 209 | * weights_radius[m].weight; |
|---|
| 210 | |
|---|
| 211 | norm += weights_core_radius[i].weight |
|---|
| 212 | * weights_length[j].weight |
|---|
| 213 | * weights_radius[m].weight |
|---|
| 214 | * weights_theta[k].weight |
|---|
| 215 | * weights_phi[l].weight; |
|---|
| 216 | |
|---|
| 217 | } |
|---|
| 218 | } |
|---|
| 219 | } |
|---|
| 220 | } |
|---|
| 221 | } |
|---|
| 222 | // Averaging in theta needs an extra normalization |
|---|
| 223 | // factor to account for the sin(theta) term in the |
|---|
| 224 | // integration (see documentation). |
|---|
| 225 | if (weights_theta.size()>1) norm = norm/asin(1.0); |
|---|
| 226 | if (vol != 0.0 || norm_vol != 0.0) { |
|---|
| 227 | //Re-normalize by avg volume |
|---|
| 228 | sum = sum*(vol/norm_vol);} |
|---|
| 229 | return sum/norm + background(); |
|---|
| 230 | } |
|---|
| 231 | |
|---|
| 232 | /** |
|---|
| 233 | * Function to evaluate 2D scattering function |
|---|
| 234 | * @param pars: parameters of the cylinder |
|---|
| 235 | * @param q: q-value |
|---|
| 236 | * @param phi: angle phi |
|---|
| 237 | * @return: function value |
|---|
| 238 | */ |
|---|
| 239 | double HollowCylinderModel :: evaluate_rphi(double q, double phi) { |
|---|
| 240 | double qx = q*cos(phi); |
|---|
| 241 | double qy = q*sin(phi); |
|---|
| 242 | return (*this).operator()(qx, qy); |
|---|
| 243 | } |
|---|
| 244 | /** |
|---|
| 245 | * Function to calculate effective radius |
|---|
| 246 | * @return: effective radius value |
|---|
| 247 | */ |
|---|
| 248 | double HollowCylinderModel :: calculate_ER() { |
|---|
| 249 | HollowCylinderParameters dp; |
|---|
| 250 | |
|---|
| 251 | dp.radius = radius(); |
|---|
| 252 | dp.length = length(); |
|---|
| 253 | |
|---|
| 254 | double rad_out = 0.0; |
|---|
| 255 | |
|---|
| 256 | // Perform the computation, with all weight points |
|---|
| 257 | double sum = 0.0; |
|---|
| 258 | double norm = 0.0; |
|---|
| 259 | |
|---|
| 260 | // Get the dispersion points for the major shell |
|---|
| 261 | vector<WeightPoint> weights_length; |
|---|
| 262 | length.get_weights(weights_length); |
|---|
| 263 | |
|---|
| 264 | // Get the dispersion points for the minor shell |
|---|
| 265 | vector<WeightPoint> weights_radius ; |
|---|
| 266 | radius.get_weights(weights_radius); |
|---|
| 267 | |
|---|
| 268 | // Loop over major shell weight points |
|---|
| 269 | for(int i=0; i< (int)weights_length.size(); i++) { |
|---|
| 270 | dp.length = weights_length[i].value; |
|---|
| 271 | for(int k=0; k< (int)weights_radius.size(); k++) { |
|---|
| 272 | dp.radius = weights_radius[k].value; |
|---|
| 273 | //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. |
|---|
| 274 | sum +=weights_length[i].weight |
|---|
| 275 | * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0; |
|---|
| 276 | norm += weights_length[i].weight* weights_radius[k].weight; |
|---|
| 277 | } |
|---|
| 278 | } |
|---|
| 279 | if (norm != 0){ |
|---|
| 280 | //return the averaged value |
|---|
| 281 | rad_out = sum/norm;} |
|---|
| 282 | else{ |
|---|
| 283 | //return normal value |
|---|
| 284 | //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER. |
|---|
| 285 | rad_out = DiamCyl(dp.length,dp.radius)/2.0;} |
|---|
| 286 | |
|---|
| 287 | return rad_out; |
|---|
| 288 | } |
|---|
