Changeset 3c102d4 in sasview for sansmodels/src/sans
- Timestamp:
- Aug 18, 2009 11:56:49 AM (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:
- eddff027
- Parents:
- 7d11b81
- Location:
- sansmodels/src/sans/models
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/TriaxialEllipsoidModel.py
r975ec8e r3c102d4 33 33 List of default parameters: 34 34 scale = 1.0 35 semi_axis B= 35.0 [A]36 semi_axis A= 100.0 [A]35 semi_axisA = 35.0 [A] 36 semi_axisB = 100.0 [A] 37 37 semi_axisC = 400.0 [A] 38 38 contrast = 5.3e-006 [1/A²] 39 39 background = 0.0 [1/cm] 40 axis_theta = 0.0 [rad]41 axis_phi = 0.0 [rad]40 axis_theta = 1.0 [rad] 41 axis_phi = 1.0 [rad] 42 42 axis_psi = 0.0 [rad] 43 43 … … 61 61 self.details = {} 62 62 self.details['scale'] = ['', None, None] 63 self.details['semi_axisA'] = ['[A]', None, None] 63 64 self.details['semi_axisB'] = ['[A]', None, None] 64 self.details['semi_axisA'] = ['[A]', None, None]65 65 self.details['semi_axisC'] = ['[A]', None, None] 66 66 self.details['contrast'] = ['[1/A²]', None, None] -
sansmodels/src/sans/models/c_extensions/parallelepiped.c
r2cb89e7 r3c102d4 98 98 */ 99 99 double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) { 100 double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y,parallel_x, parallel_y, parallel_z;100 double parallel_x, parallel_y, parallel_z; 101 101 double q_z; 102 102 double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC; … … 111 111 112 112 // parallelepiped c axis orientation 113 cparallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi);114 cparallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi);115 cparallel_z = cos(pars->parallel_theta);113 parallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi); 114 parallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi); 115 parallel_z = cos(pars->parallel_theta); 116 116 117 117 // q vector … … 120 120 // Compute the angle btw vector q and the 121 121 // axis of the parallelepiped 122 cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;122 cos_val_c = parallel_x*q_x + parallel_y*q_y + parallel_z*q_z; 123 123 alpha = acos(cos_val_c); 124 124 … … 126 126 parallel_x = -(1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi); 127 127 parallel_y = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi); 128 129 //parallel_x = -(1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi); 130 //parallel_y = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi); 128 131 cos_val_a = (parallel_x*q_x + parallel_y*q_y); 129 132 … … 131 134 132 135 // parallelepiped b axis orientation 133 bparallel_x = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);134 bparallel_y = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);136 parallel_x = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi); 137 parallel_y = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi); 135 138 // axis of the parallelepiped 136 cos_val_b = ( bparallel_x*q_x + bparallel_y*q_y) ;139 cos_val_b = (parallel_x*q_x + parallel_y*q_y) ; 137 140 138 141 -
sansmodels/src/sans/models/c_extensions/stacked_disks.c
r5068697 r3c102d4 19 19 double stacked_disks_analytical_1D(StackedDisksParameters *pars, double q) { 20 20 double dp[10]; 21 21 22 22 // Fill paramater array 23 23 dp[0] = pars->scale; 24 24 dp[1] = pars->radius; 25 dp[2] = pars-> length;26 dp[3] = pars-> thickness;25 dp[2] = pars->core_thick; 26 dp[3] = pars->layer_thick; 27 27 dp[4] = pars->core_sld; 28 28 dp[5] = pars->layer_sld; 29 29 dp[6] = pars->solvent_sld; 30 dp[7] = pars->n layers;31 dp[8] = pars->s pacing;30 dp[7] = pars->n_stacking; 31 dp[8] = pars->sigma_d; 32 32 dp[9] = pars->background; 33 33 34 34 // Call library function to evaluate model 35 return StackedDiscs(dp, q); 35 return StackedDiscs(dp, q); 36 36 } 37 37 /** … … 45 45 q = sqrt(qx*qx+qy*qy); 46 46 return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q); 47 } 47 } 48 48 49 49 … … 57 57 double stacked_disks_analytical_2D(StackedDisksParameters *pars, double q, double phi) { 58 58 return stacked_disks_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 59 } 60 59 } 60 61 61 /** 62 62 * Function to evaluate 2D scattering function … … 71 71 double q_z; 72 72 double alpha, vol, cos_val; 73 double d, dum ;73 double d, dum, halfheight; 74 74 double answer; 75 76 75 76 77 77 78 78 // parallelepiped orientation … … 80 80 cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi); 81 81 cyl_z = cos(pars->axis_theta); 82 82 83 83 // q vector 84 84 q_z = 0; 85 85 86 86 // Compute the angle btw vector q and the 87 87 // axis of the parallelepiped 88 88 cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 89 89 90 90 // The following test should always pass 91 91 if (fabs(cos_val)>1.0) { … … 93 93 return 0; 94 94 } 95 95 96 96 // Note: cos(alpha) = 0 and 1 will get an 97 97 // undefined value from Stackdisc_kern … … 99 99 100 100 // Call the IGOR library function to get the kernel 101 dum =0.1; 102 d= 0.1; 101 d = 2 * pars->layer_thick + pars->core_thick; 102 halfheight = pars->core_thick/2.0; 103 dum =alpha ; 103 104 answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld, 104 pars->solvent_sld, pars->length,pars->thickness, dum, pars->spacing, d,pars->nlayers);105 105 pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha); 106 106 107 // Multiply by contrast^2 107 108 //answer *= pars->contrast*pars->contrast; 108 109 109 110 //normalize by staked disks volume 110 //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel 111 vol = acos(-1.0) * pars->radius * pars->radius * pars->length; 112 answer *= vol; 113 111 vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking; 112 answer /= vol; 113 114 114 //convert to [cm-1] 115 115 answer *= 1.0e8; 116 116 117 117 //Scale 118 118 answer *= pars->scale; 119 119 120 120 // add in the background 121 121 answer += pars->background; 122 122 123 123 return answer; 124 124 } 125 126 125 126 -
sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.c
r975ec8e r3c102d4 35 35 double t,a,b,c; 36 36 double kernel; 37 double pi = acos(-1.0);37 double pi = 4.0*atan(1.0); 38 38 39 39 a = pars->semi_axisA ; … … 42 42 43 43 t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha)); 44 if (t==0 ){44 if (t==0.0){ 45 45 kernel = 1.0; 46 46 }else{ 47 kernel = 3 *(sin(t)-t*cos(t))/(t*t*t);47 kernel = 3.0*(sin(t)-t*cos(t))/(t*t*t); 48 48 } 49 49 return kernel*kernel; … … 137 137 //normalize by cylinder volume 138 138 //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 139 vol = 4 /3 * pi* pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;139 vol = 4.0* pi/3.0 * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC; 140 140 answer *= vol; 141 142 141 //convert to [cm-1] 143 142 answer *= 1.0e8; 144 145 143 //Scale 146 144 answer *= pars->scale; -
sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.h
r975ec8e r3c102d4 17 17 // [DEFAULT]=scale=1.0 18 18 double scale; 19 /// semi -axis Bof the triaxial_ellipsoid [A]20 // [DEFAULT]=semi_axis B= 35.0 [A]19 /// semi -axis A of the triaxial_ellipsoid [A] 20 // [DEFAULT]=semi_axisA= 35.0 [A] 21 21 double semi_axisA; 22 22 /// semi -axis B of the triaxial_ellipsoid [A] 23 // [DEFAULT]=semi_axis A=100.0 [A]23 // [DEFAULT]=semi_axisB=100.0 [A] 24 24 double semi_axisB; 25 25 /// semi -axis C of the triaxial_ellipsoid [A] … … 33 33 double background; 34 34 /// Orientation of the triaxial_ellipsoid axis w/respect incoming beam [rad] 35 // [DEFAULT]=axis_theta= 0.0 [rad]35 // [DEFAULT]=axis_theta=1.0 [rad] 36 36 double axis_theta; 37 37 /// Orientation of the triaxial_ellipsoid in the plane of the detector [rad] 38 // [DEFAULT]=axis_phi= 0.0 [rad]38 // [DEFAULT]=axis_phi=1.0 [rad] 39 39 double axis_phi; 40 40 /// Orientation of the cross section of the triaxial_ellipsoid in the plane of the detector [rad] -
sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp
r975ec8e r3c102d4 88 88 PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 89 89 PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.000000)); 90 PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d", 0.000000));91 PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d", 100.000000));92 PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d", 35.000000));90 PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000)); 91 PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",35.000000)); 92 PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",100.000000)); 93 93 PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.000000)); 94 PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d", 0.000000));94 PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000)); 95 95 PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000)); 96 96 PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005)); -
sansmodels/src/sans/models/c_models/parallelepiped.cpp
r975ec8e r3c102d4 35 35 ParallelepipedModel :: ParallelepipedModel() { 36 36 scale = Parameter(1.0); 37 short_ edgeA= Parameter(35.0, true);38 short_ edgeA.set_max(1.0);39 long er_edgeB= Parameter(75.0, true);40 long er_edgeB.set_min(1.0);41 long uest_edgeC= Parameter(400.0, true);42 long uest_edgeC.set_min(1.0);37 short_a = Parameter(35.0, true); 38 short_a.set_max(1.0); 39 long_b = Parameter(75.0, true); 40 long_b.set_min(1.0); 41 longer_c = Parameter(400.0, true); 42 longer_c.set_min(1.0); 43 43 contrast = Parameter(53.e-7); 44 44 background = Parameter(0.0); … … 60 60 // Add the background after averaging 61 61 dp[0] = scale(); 62 dp[1] = short_ edgeA();63 dp[2] = long er_edgeB();64 dp[3] = long uest_edgeC();62 dp[1] = short_a(); 63 dp[2] = long_b(); 64 dp[3] = longer_c(); 65 65 dp[4] = contrast(); 66 66 dp[5] = 0.0; 67 67 68 68 // Get the dispersion points for the short_edgeA 69 vector<WeightPoint> weights_short_ edgeA;70 short_ edgeA.get_weights(weights_short_edgeA);69 vector<WeightPoint> weights_short_a; 70 short_a.get_weights(weights_short_a); 71 71 72 72 // Get the dispersion points for the longer_edgeB 73 vector<WeightPoint> weights_long er_edgeB;74 long er_edgeB.get_weights(weights_longer_edgeB);73 vector<WeightPoint> weights_long_b; 74 long_b.get_weights(weights_long_b); 75 75 76 76 // Get the dispersion points for the longuest_edgeC 77 vector<WeightPoint> weights_long uest_edgeC;78 long uest_edgeC.get_weights(weights_longuest_edgeC);77 vector<WeightPoint> weights_longer_c; 78 longer_c.get_weights(weights_longer_c); 79 79 80 80 … … 85 85 86 86 // Loop over short_edgeA weight points 87 for(int i=0; i< (int)weights_short_ edgeA.size(); i++) {88 dp[1] = weights_short_ edgeA[i].value;87 for(int i=0; i< (int)weights_short_a.size(); i++) { 88 dp[1] = weights_short_a[i].value; 89 89 90 90 // Loop over longer_edgeB weight points 91 for(int j=0; j< (int)weights_long er_edgeB.size(); j++) {92 dp[2] = weights_long er_edgeB[i].value;91 for(int j=0; j< (int)weights_long_b.size(); j++) { 92 dp[2] = weights_long_b[j].value; 93 93 94 94 // Loop over longuest_edgeC weight points 95 for(int k=0; k< (int)weights_long uest_edgeC.size(); k++) {96 dp[3] = weights_long uest_edgeC[j].value;97 98 sum += weights_short_ edgeA[i].weight * weights_longer_edgeB[j].weight99 * weights_long uest_edgeC[k].weight * Parallelepiped(dp, q);100 101 norm += weights_short_ edgeA[i].weight102 * weights_long er_edgeB[j].weight * weights_longuest_edgeC[k].weight;95 for(int k=0; k< (int)weights_longer_c.size(); k++) { 96 dp[3] = weights_longer_c[k].value; 97 98 sum += weights_short_a[i].weight * weights_long_b[j].weight 99 * weights_longer_c[k].weight * Parallelepiped(dp, q); 100 101 norm += weights_short_a[i].weight 102 * weights_long_b[j].weight * weights_longer_c[k].weight; 103 103 } 104 104 } … … 116 116 // Fill parameter array 117 117 dp.scale = scale(); 118 dp.short_ edgeA = short_edgeA();119 dp.long er_edgeB = longer_edgeB();120 dp.long uest_edgeC = longuest_edgeC();118 dp.short_a = short_a(); 119 dp.long_b = long_b(); 120 dp.longer_c = longer_c(); 121 121 dp.contrast = contrast(); 122 122 dp.background = 0.0; … … 128 128 129 129 // Get the dispersion points for the short_edgeA 130 vector<WeightPoint> weights_short_ edgeA;131 short_ edgeA.get_weights(weights_short_edgeA);130 vector<WeightPoint> weights_short_a; 131 short_a.get_weights(weights_short_a); 132 132 133 133 // Get the dispersion points for the longer_edgeB 134 vector<WeightPoint> weights_long er_edgeB;135 long er_edgeB.get_weights(weights_longer_edgeB);134 vector<WeightPoint> weights_long_b; 135 long_b.get_weights(weights_long_b); 136 136 137 137 // Get angular averaging for the longuest_edgeC 138 vector<WeightPoint> weights_long uest_edgeC;139 long uest_edgeC.get_weights(weights_longuest_edgeC);138 vector<WeightPoint> weights_longer_c; 139 longer_c.get_weights(weights_longer_c); 140 140 141 141 // Get angular averaging for theta … … 156 156 157 157 // Loop over radius weight points 158 for(int i=0; i< (int)weights_short_ edgeA.size(); i++) {159 dp.short_ edgeA = weights_short_edgeA[i].value;158 for(int i=0; i< (int)weights_short_a.size(); i++) { 159 dp.short_a = weights_short_a[i].value; 160 160 161 161 // Loop over longer_edgeB weight points 162 for(int j=0; j< (int)weights_long er_edgeB.size(); j++) {163 dp.long er_edgeB = weights_longer_edgeB[j].value;162 for(int j=0; j< (int)weights_long_b.size(); j++) { 163 dp.long_b = weights_long_b[j].value; 164 164 165 165 // Average over longuest_edgeC distribution 166 for(int k=0; k< (int)weights_long uest_edgeC.size(); k++) {167 dp.long uest_edgeC = weights_longuest_edgeC[k].value;166 for(int k=0; k< (int)weights_longer_c.size(); k++) { 167 dp.longer_c = weights_longer_c[k].value; 168 168 169 169 // Average over theta distribution … … 179 179 dp.parallel_psi = weights_parallel_psi[n].value; 180 180 181 double _ptvalue = weights_short_ edgeA[i].weight182 * weights_long er_edgeB[j].weight183 * weights_long uest_edgeC[k].weight181 double _ptvalue = weights_short_a[i].weight 182 * weights_long_b[j].weight 183 * weights_longer_c[k].weight 184 184 * weights_parallel_theta[l].weight 185 185 * weights_parallel_phi[m].weight … … 191 191 sum += _ptvalue; 192 192 193 norm += weights_short_ edgeA[i].weight194 * weights_long er_edgeB[j].weight195 * weights_long uest_edgeC[k].weight193 norm += weights_short_a[i].weight 194 * weights_long_b[j].weight 195 * weights_longer_c[k].weight 196 196 * weights_parallel_theta[l].weight 197 197 * weights_parallel_phi[m].weight -
sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp
r975ec8e r3c102d4 34 34 TriaxialEllipsoidModel :: TriaxialEllipsoidModel() { 35 35 scale = Parameter(1.0); 36 semi_axisA = Parameter( 20.0, true);36 semi_axisA = Parameter(35.0, true); 37 37 semi_axisA.set_min(0.0); 38 semi_axisB = Parameter( 20.0, true);38 semi_axisB = Parameter(100.0, true); 39 39 semi_axisB.set_min(0.0); 40 40 semi_axisC = Parameter(400.0, true); … … 42 42 contrast = Parameter(5.3e-6); 43 43 background = Parameter(0.0); 44 axis_theta = Parameter( 0.0, true);45 axis_phi = Parameter( 0.0, true);44 axis_theta = Parameter(1.0, true); 45 axis_phi = Parameter(1.0, true); 46 46 axis_psi = Parameter(0.0, true); 47 47 } … … 190 190 * weights_theta[l].weight 191 191 * weights_phi[m].weight 192 * weights_psi[ m].weight;192 * weights_psi[n].weight; 193 193 } 194 194 }
Note: See TracChangeset
for help on using the changeset viewer.