Changeset ea07075 in sasview for sansmodels/src
- Timestamp:
- Aug 3, 2009 5:41:22 PM (15 years ago)
- Branches:
- master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- e6fa43e
- Parents:
- 2a1be2d9
- Location:
- sansmodels/src/sans/models
- Files:
-
- 4 edited
-
c_extensions/flexible_cylinder.c (modified) (4 diffs)
-
c_extensions/flexible_cylinder.h (modified) (3 diffs)
-
c_models/CFlexibleCylinderModel.cpp (modified) (7 diffs)
-
c_models/flexiblecylinder.cpp (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/c_extensions/flexible_cylinder.c
r5068697 rea07075 19 19 double flexible_cylinder_analytical_1D(FlexibleCylinderParameters *pars, double q) { 20 20 double dp[6]; 21 21 22 22 // Fill paramater array 23 23 dp[0] = pars->scale; … … 29 29 30 30 // Call library function to evaluate model 31 return FlexExclVolCyl(dp, q); 31 return FlexExclVolCyl(dp, q); 32 32 } 33 33 /** … … 40 40 double q; 41 41 q = sqrt(qx*qx+qy*qy); 42 return flexible_cylinder_analytical_ 2D_scaled(pars, q, qx/q, qy/q);43 } 42 return flexible_cylinder_analytical_1D(pars, q); 43 } 44 44 45 45 … … 52 52 */ 53 53 double flexible_cylinder_analytical_2D(FlexibleCylinderParameters *pars, double q, double phi) { 54 return flexible_cylinder_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 55 } 56 57 /** 58 * Function to evaluate 2D scattering function 59 * @param pars: parameters of the flexible cylinder 60 * @param q: q-value 61 * @param q_x: q_x / q 62 * @param q_y: q_y / q 63 * @return: function value 64 */ 65 double flexible_cylinder_analytical_2D_scaled(FlexibleCylinderParameters *pars, double q, double q_x, double q_y) { 66 double cyl_x, cyl_y, cyl_z; 67 double q_z; 68 double alpha, vol, cos_val; 69 double answer; 70 71 54 return flexible_cylinder_analytical_1D(pars, q); 55 } 72 56 73 // parallelepiped orientation74 cyl_x = sin(pars->axis_theta) * cos(pars->axis_phi);75 cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi);76 cyl_z = cos(pars->axis_theta);77 78 // q vector79 q_z = 0;80 81 // Compute the angle btw vector q and the82 // axis of the parallelepiped83 cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;84 85 // The following test should always pass86 if (fabs(cos_val)>1.0) {87 printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");88 return 0;89 }90 91 // Note: cos(alpha) = 0 and 1 will get an92 // undefined value from PPKernel93 alpha = acos( cos_val );94 57 95 // Call the IGOR library function to get the kernel96 answer = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);97 98 // Multiply by contrast^299 answer *= pars->contrast*pars->contrast;100 101 //normalize by cylinder volume102 //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel103 vol = acos(-1.0) * pars->radius * pars->radius * pars->length;104 answer *= vol;105 106 //convert to [cm-1]107 answer *= 1.0e8;108 109 //Scale110 answer *= pars->scale;111 112 // add in the background113 answer += pars->background;114 115 return answer;116 }117 118 -
sansmodels/src/sans/models/c_extensions/flexible_cylinder.h
r5068697 rea07075 7 7 /** Structure definition for Flexible cylinder parameters 8 8 * [PYTHONCLASS] = FlexibleCylinderModel 9 * [DISP_PARAMS] = length, radius, axis_theta, axis_phi 10 [DESCRIPTION] = <text> Note : scale and contrast are both multiplicative factors in the model and are perfectly 11 correlated. One or both of these parameters must be held fixed during model fitting. 9 * [DISP_PARAMS] = length, kuhn_length, radius 10 [DESCRIPTION] = <text> Note : 'scale' and 'contrast' are both multiplicative factors in the 11 model and are perfectly correlated. One or 12 both of these parameters must be held fixed 13 during model fitting. 12 14 </text> 13 [FIXED]= <text>length.width; radius.width; axis_theta.width; axis_phi.width</text>14 [ORIENTATION_PARAMS]= <text> axis_phi; axis_theta; axis_phi.width; axis_theta.width</text>15 [FIXED]= <text>length.width; radius.width; , kuhn_length.width</text> 16 [ORIENTATION_PARAMS]= <text></text> 15 17 16 18 … … 32 34 // [DEFAULT]=contrast=5.3e-6 [1/A²] 33 35 double contrast; 34 /// Incoherent Background [1/cm] 36 /// Incoherent Background [1/cm] 35 37 // [DEFAULT]=background=0.0001 [1/cm] 36 38 double background; 37 /// Orientation of the flexible cylinder axis w/respect incoming beam [rad]38 // [DEFAULT]=axis_theta=1.0 [rad]39 double axis_theta;40 /// Orientation of the flexible cylinder in the plane of the detector [rad]41 // [DEFAULT]=axis_phi=1.0 [rad]42 double axis_phi;43 39 44 40 … … 53 49 double flexible_cylinder_analytical_2D(FlexibleCylinderParameters *pars, double q, double phi); 54 50 double flexible_cylinder_analytical_2DXY(FlexibleCylinderParameters *pars, double qx, double qy); 55 double flexible_cylinder_analytical_2D_scaled(FlexibleCylinderParameters *pars, double q, double q_x, double q_y);56 51 57 52 #endif -
sansmodels/src/sans/models/c_models/CFlexibleCylinderModel.cpp
r5068697 rea07075 84 84 // Initialize parameter dictionary 85 85 PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 86 PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000));87 86 PyDict_SetItemString(self->params,"length",Py_BuildValue("d",1000.000000)); 88 PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000));89 87 PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",20.000000)); 90 88 PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000100)); … … 98 96 PyDict_SetItemString(self->dispersion, "length", disp_dict); 99 97 disp_dict = PyDict_New(); 98 self->model->kuhn_length.dispersion->accept_as_source(visitor, self->model->kuhn_length.dispersion, disp_dict); 99 PyDict_SetItemString(self->dispersion, "kuhn_length", disp_dict); 100 disp_dict = PyDict_New(); 100 101 self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict); 101 102 PyDict_SetItemString(self->dispersion, "radius", disp_dict); 102 disp_dict = PyDict_New();103 self->model->axis_theta.dispersion->accept_as_source(visitor, self->model->axis_theta.dispersion, disp_dict);104 PyDict_SetItemString(self->dispersion, "axis_theta", disp_dict);105 disp_dict = PyDict_New();106 self->model->axis_phi.dispersion->accept_as_source(visitor, self->model->axis_phi.dispersion, disp_dict);107 PyDict_SetItemString(self->dispersion, "axis_phi", disp_dict);108 103 109 104 … … 159 154 // Reader parameter dictionary 160 155 self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 161 self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") );162 156 self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 163 self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") );164 157 self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") ); 165 158 self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); … … 171 164 disp_dict = PyDict_GetItemString(self->dispersion, "length"); 172 165 self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 166 disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length"); 167 self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict); 173 168 disp_dict = PyDict_GetItemString(self->dispersion, "radius"); 174 169 self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict); 175 disp_dict = PyDict_GetItemString(self->dispersion, "axis_theta");176 self->model->axis_theta.dispersion->accept_as_destination(visitor, self->model->axis_theta.dispersion, disp_dict);177 disp_dict = PyDict_GetItemString(self->dispersion, "axis_phi");178 self->model->axis_phi.dispersion->accept_as_destination(visitor, self->model->axis_phi.dispersion, disp_dict);179 170 180 171 … … 229 220 // Reader parameter dictionary 230 221 self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 231 self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") );232 222 self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 233 self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") );234 223 self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") ); 235 224 self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); … … 241 230 disp_dict = PyDict_GetItemString(self->dispersion, "length"); 242 231 self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 232 disp_dict = PyDict_GetItemString(self->dispersion, "kuhn_length"); 233 self->model->kuhn_length.dispersion->accept_as_destination(visitor, self->model->kuhn_length.dispersion, disp_dict); 243 234 disp_dict = PyDict_GetItemString(self->dispersion, "radius"); 244 235 self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict); 245 disp_dict = PyDict_GetItemString(self->dispersion, "axis_theta");246 self->model->axis_theta.dispersion->accept_as_destination(visitor, self->model->axis_theta.dispersion, disp_dict);247 disp_dict = PyDict_GetItemString(self->dispersion, "axis_phi");248 self->model->axis_phi.dispersion->accept_as_destination(visitor, self->model->axis_phi.dispersion, disp_dict);249 236 250 237 … … 304 291 if (!strcmp(par_name, "length")) { 305 292 self->model->length.dispersion = dispersion; 293 } else if (!strcmp(par_name, "kuhn_length")) { 294 self->model->kuhn_length.dispersion = dispersion; 306 295 } else if (!strcmp(par_name, "radius")) { 307 296 self->model->radius.dispersion = dispersion; 308 } else if (!strcmp(par_name, "axis_theta")) {309 self->model->axis_theta.dispersion = dispersion;310 } else if (!strcmp(par_name, "axis_phi")) {311 self->model->axis_phi.dispersion = dispersion;312 297 } else { 313 298 PyErr_SetString(CFlexibleCylinderModelError, -
sansmodels/src/sans/models/c_models/flexiblecylinder.cpp
r5068697 rea07075 19 19 * 20 20 * TODO: refactor so that we pull in the old sansmodels.c_extensions 21 * TODO: add 2d 21 * TODO: add 2d 22 22 */ 23 23 … … 43 43 contrast = Parameter(5.3e-6); 44 44 background = Parameter(0.0001); 45 axis_theta = Parameter(0.0, true);46 axis_phi = Parameter(0.0, true);47 45 } 48 46 … … 101 99 */ 102 100 double FlexibleCylinderModel :: operator()(double qx, double qy) { 103 FlexibleCylinderParameters dp; 104 // Fill parameter array 105 dp.scale = scale(); 106 dp.length = length(); 107 dp.kuhn_length= kuhn_length(); 108 dp.radius = radius(); 109 dp.contrast = contrast(); 110 dp.background = background(); 111 dp.axis_theta = axis_theta(); 112 dp.axis_phi = axis_phi(); 113 114 // Get the dispersion points for the length 115 vector<WeightPoint> weights_length; 116 length.get_weights(weights_length); 117 118 // Get the dispersion points for the radius 119 vector<WeightPoint> weights_radius; 120 radius.get_weights(weights_radius); 121 122 // Get angular averaging for theta 123 vector<WeightPoint> weights_theta; 124 axis_theta.get_weights(weights_theta); 125 126 // Get angular averaging for phi 127 vector<WeightPoint> weights_phi; 128 axis_phi.get_weights(weights_phi); 129 130 // Perform the computation, with all weight points 131 double sum = 0.0; 132 double norm = 0.0; 133 134 // Loop over length weight points 135 for(int i=0; i< (int)weights_length.size(); i++) { 136 dp.length = weights_length[i].value; 137 138 // Loop over radius weight points 139 for(int j=0; j< (int)weights_radius.size(); j++) { 140 dp.radius = weights_radius[j].value; 141 142 // Average over theta distribution 143 for(int k=0; k< (int)weights_theta.size(); k++) { 144 dp.axis_theta = weights_theta[k].value; 145 146 // Average over phi distribution 147 for(int l=0; l <(int)weights_phi.size(); l++) { 148 dp.axis_phi = weights_phi[l].value; 149 150 double _ptvalue = weights_length[i].weight 151 * weights_radius[j].weight 152 * weights_theta[k].weight 153 * weights_phi[l].weight 154 * flexible_cylinder_analytical_2DXY(&dp, qx, qy); 155 if (weights_theta.size()>1) { 156 _ptvalue *= sin(weights_theta[k].value); 157 } 158 sum += _ptvalue; 159 160 norm += weights_length[i].weight 161 * weights_radius[j].weight 162 * weights_theta[k].weight 163 * weights_phi[l].weight; 164 } 165 } 166 } 167 } 168 // Averaging in theta needs an extra normalization 169 // factor to account for the sin(theta) term in the 170 // integration (see documentation). 171 if (weights_theta.size()>1) norm = norm / asin(1.0); 172 return sum/norm + background(); 101 double q = sqrt(qx*qx + qy*qy); 102 return (*this).operator()(q); 173 103 } 174 104 … … 181 111 */ 182 112 double FlexibleCylinderModel :: evaluate_rphi(double q, double phi) { 183 double qx = q*cos(phi);184 double qy = q*sin(phi);185 return (*this).operator()(q x, qy);113 //double qx = q*cos(phi); 114 //double qy = q*sin(phi); 115 return (*this).operator()(q); 186 116 }
Note: See TracChangeset
for help on using the changeset viewer.
