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