Changeset 3c102d4 in sasview for sansmodels/src/sans/models

Timestamp:

Aug 18, 2009 11:56:49 AM (15 years ago)

Author:

Jae Cho <jhjcho@…>

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

Message:

fixed problems in 2d

Location:

sansmodels/src/sans/models

Files:

• TriaxialEllipsoidModel.py (modified) (2 diffs)
• c_extensions/parallelepiped.c (modified) (5 diffs)
• c_extensions/stacked_disks.c (modified) (7 diffs)
• c_extensions/triaxial_ellipsoid.c (modified) (3 diffs)
• c_extensions/triaxial_ellipsoid.h (modified) (2 diffs)
• c_models/CTriaxialEllipsoidModel.cpp (modified) (1 diff)
• c_models/parallelepiped.cpp (modified) (8 diffs)
• c_models/triaxialellipsoid.cpp (modified) (3 diffs)

sansmodels/src/sans/models/TriaxialEllipsoidModel.py

@@ -33 +33 @@
         List of default parameters:
          scale           = 1.0 
-         semi_axisB      = 35.0 [A]
-         semi_axisA      = 100.0 [A]
+         semi_axisA      = 35.0 [A]
+         semi_axisB      = 100.0 [A]
          semi_axisC      = 400.0 [A]
          contrast        = 5.3e-006 [1/A²]
          background      = 0.0 [1/cm]
-         axis_theta      = 0.0 [rad]
-         axis_phi        = 0.0 [rad]
+         axis_theta      = 1.0 [rad]
+         axis_phi        = 1.0 [rad]
          axis_psi        = 0.0 [rad]
 
@@ -61 +61 @@
         self.details = {}
         self.details['scale'] = ['', None, None]
+        self.details['semi_axisA'] = ['[A]', None, None]
         self.details['semi_axisB'] = ['[A]', None, None]
-        self.details['semi_axisA'] = ['[A]', None, None]
         self.details['semi_axisC'] = ['[A]', None, None]
         self.details['contrast'] = ['[1/A²]', None, None]

sansmodels/src/sans/models/c_extensions/parallelepiped.c

@@ -98 +98 @@
  */
 double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) {
-        double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y, parallel_z;
+        double parallel_x, parallel_y, parallel_z;
         double q_z;
         double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
@@ -111 +111 @@
 
     // parallelepiped c axis orientation
-    cparallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi);
-    cparallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi);
-    cparallel_z = cos(pars->parallel_theta);
+    parallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi);
+    parallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi);
+    parallel_z = cos(pars->parallel_theta);
 
     // q vector
@@ -120 +120 @@
     // Compute the angle btw vector q and the
     // axis of the parallelepiped
-    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
+    cos_val_c = parallel_x*q_x + parallel_y*q_y + parallel_z*q_z;
     alpha = acos(cos_val_c);
 
@@ -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);
     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);
+
+    //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);
+    //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);
     cos_val_a = (parallel_x*q_x + parallel_y*q_y);
 
@@ -131 +134 @@
 
     // parallelepiped b axis orientation
-    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);
-    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);
+    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);
+    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);
     // axis of the parallelepiped
-    cos_val_b = (bparallel_x*q_x + bparallel_y*q_y) ;
+    cos_val_b = (parallel_x*q_x + parallel_y*q_y) ;
 
 

sansmodels/src/sans/models/c_extensions/stacked_disks.c

@@ -19 +19 @@
 double stacked_disks_analytical_1D(StackedDisksParameters *pars, double q) {
         double dp[10];
-        
+
         // Fill paramater array
         dp[0] = pars->scale;
         dp[1] = pars->radius;
-        dp[2] = pars->length;
-        dp[3] = pars->thickness;
+        dp[2] = pars->core_thick;
+        dp[3] = pars->layer_thick;
         dp[4] = pars->core_sld;
         dp[5] = pars->layer_sld;
         dp[6] = pars->solvent_sld;
-        dp[7] = pars->nlayers;
-        dp[8] = pars->spacing;
+        dp[7] = pars->n_stacking;
+        dp[8] = pars->sigma_d;
         dp[9] = pars->background;
 
         // Call library function to evaluate model
-        return StackedDiscs(dp, q);     
+        return StackedDiscs(dp, q);
 }
 /**
@@ -45 +45 @@
         q = sqrt(qx*qx+qy*qy);
     return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q);
-} 
+}
 
 
@@ -57 +57 @@
 double stacked_disks_analytical_2D(StackedDisksParameters *pars, double q, double phi) {
     return stacked_disks_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
-} 
-        
+}
+
 /**
  * Function to evaluate 2D scattering function
@@ -71 +71 @@
         double q_z;
         double alpha, vol, cos_val;
-        double d, dum;
+        double d, dum, halfheight;
         double answer;
-        
-        
+
+
 
     // parallelepiped orientation
@@ -80 +80 @@
     cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi);
     cyl_z = cos(pars->axis_theta);
-     
+
     // q vector
     q_z = 0;
-        
+
     // Compute the angle btw vector q and the
     // axis of the parallelepiped
     cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
-    
+
     // The following test should always pass
     if (fabs(cos_val)>1.0) {
@@ -93 +93 @@
         return 0;
     }
-    
+
     // Note: cos(alpha) = 0 and 1 will get an
     // undefined value from Stackdisc_kern
@@ -99 +99 @@
 
         // Call the IGOR library function to get the kernel
-        dum =0.1;
-        d= 0.1;
+        d = 2 * pars->layer_thick + pars->core_thick;
+        halfheight = pars->core_thick/2.0;
+        dum =alpha ;
         answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld,
-                pars->solvent_sld,pars->length,pars->thickness, dum, pars->spacing, d,pars->nlayers);
-        
+                pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha);
+
         // Multiply by contrast^2
         //answer *= pars->contrast*pars->contrast;
-        
+
         //normalize by staked disks volume
-        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
-    vol = acos(-1.0) * pars->radius * pars->radius * pars->length;
-        answer *= vol;
-        
+    vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking;
+        answer /= vol;
+
         //convert to [cm-1]
         answer *= 1.0e8;
-        
+
         //Scale
         answer *= pars->scale;
-        
+
         // add in the background
         answer += pars->background;
-        
+
         return answer;
 }
-    
 
+

sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.c

@@ -35 +35 @@
         double t,a,b,c;
         double kernel;
-        double pi = acos(-1.0);
+        double pi = 4.0*atan(1.0);
 
         a = pars->semi_axisA ;
@@ -42 +42 @@
 
         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));
-        if (t==0){
+        if (t==0.0){
                 kernel  = 1.0;
         }else{
-                kernel  = 3*(sin(t)-t*cos(t))/(t*t*t);
+                kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t);
         }
         return kernel*kernel;
@@ -137 +137 @@
         //normalize by cylinder volume
         //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    vol = 4/3 * pi * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;
+    vol = 4.0* pi/3.0  * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;
         answer *= vol;
-
         //convert to [cm-1]
         answer *= 1.0e8;
-
         //Scale
         answer *= pars->scale;

sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.h

@@ -17 +17 @@
     //  [DEFAULT]=scale=1.0
     double scale;
-    /// semi -axis B of the triaxial_ellipsoid [A]
-    //  [DEFAULT]=semi_axisB= 35.0 [A]
+    /// semi -axis A of the triaxial_ellipsoid [A]
+    //  [DEFAULT]=semi_axisA= 35.0 [A]
     double semi_axisA;
     /// semi -axis B of the triaxial_ellipsoid [A]
-    //  [DEFAULT]=semi_axisA=100.0 [A]
+    //  [DEFAULT]=semi_axisB=100.0 [A]
     double semi_axisB;
         /// semi -axis C of the triaxial_ellipsoid [A]
@@ -33 +33 @@
         double background;
     /// Orientation of the triaxial_ellipsoid axis w/respect incoming beam [rad]
-    //  [DEFAULT]=axis_theta=0.0 [rad]
+    //  [DEFAULT]=axis_theta=1.0 [rad]
     double axis_theta;
     /// Orientation of the triaxial_ellipsoid in the plane of the detector [rad]
-    //  [DEFAULT]=axis_phi=0.0 [rad]
+    //  [DEFAULT]=axis_phi=1.0 [rad]
     double axis_phi;
     /// Orientation of the cross section of the triaxial_ellipsoid in the plane of the detector [rad]

sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp

@@ -88 +88 @@
         PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000));
         PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000));
-        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",100.000000));
-        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",35.000000));
+        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000));
+        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",35.000000));
+        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",100.000000));
         PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.000000));
-        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000));
+        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000));
         PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000));
         PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005));

sansmodels/src/sans/models/c_models/parallelepiped.cpp

@@ -35 +35 @@
 ParallelepipedModel :: ParallelepipedModel() {
         scale      = Parameter(1.0);
-        short_edgeA     = Parameter(35.0, true);
-        short_edgeA.set_max(1.0);
-        longer_edgeB     = Parameter(75.0, true);
-        longer_edgeB.set_min(1.0);
-        longuest_edgeC     = Parameter(400.0, true);
-        longuest_edgeC.set_min(1.0);
+        short_a     = Parameter(35.0, true);
+        short_a.set_max(1.0);
+        long_b     = Parameter(75.0, true);
+        long_b.set_min(1.0);
+        longer_c     = Parameter(400.0, true);
+        longer_c.set_min(1.0);
         contrast   = Parameter(53.e-7);
         background = Parameter(0.0);
@@ -60 +60 @@
         // Add the background after averaging
         dp[0] = scale();
-        dp[1] = short_edgeA();
-        dp[2] = longer_edgeB();
-        dp[3] = longuest_edgeC();
+        dp[1] = short_a();
+        dp[2] = long_b();
+        dp[3] = longer_c();
         dp[4] = contrast();
         dp[5] = 0.0;
 
         // Get the dispersion points for the short_edgeA
-        vector<WeightPoint> weights_short_edgeA;
-        short_edgeA.get_weights(weights_short_edgeA);
+        vector<WeightPoint> weights_short_a;
+        short_a.get_weights(weights_short_a);
 
         // Get the dispersion points for the longer_edgeB
-        vector<WeightPoint> weights_longer_edgeB;
-        longer_edgeB.get_weights(weights_longer_edgeB);
+        vector<WeightPoint> weights_long_b;
+        long_b.get_weights(weights_long_b);
 
         // Get the dispersion points for the longuest_edgeC
-        vector<WeightPoint> weights_longuest_edgeC;
-        longuest_edgeC.get_weights(weights_longuest_edgeC);
+        vector<WeightPoint> weights_longer_c;
+        longer_c.get_weights(weights_longer_c);
 
 
@@ -85 +85 @@
 
         // Loop over short_edgeA weight points
-        for(int i=0; i< (int)weights_short_edgeA.size(); i++) {
-                dp[1] = weights_short_edgeA[i].value;
+        for(int i=0; i< (int)weights_short_a.size(); i++) {
+                dp[1] = weights_short_a[i].value;
 
                 // Loop over longer_edgeB weight points
-                for(int j=0; j< (int)weights_longer_edgeB.size(); j++) {
-                        dp[2] = weights_longer_edgeB[i].value;
+                for(int j=0; j< (int)weights_long_b.size(); j++) {
+                        dp[2] = weights_long_b[j].value;
 
                         // Loop over longuest_edgeC weight points
-                        for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) {
-                                dp[3] = weights_longuest_edgeC[j].value;
-
-                                sum += weights_short_edgeA[i].weight * weights_longer_edgeB[j].weight
-                                        * weights_longuest_edgeC[k].weight * Parallelepiped(dp, q);
-
-                                norm += weights_short_edgeA[i].weight
-                                         * weights_longer_edgeB[j].weight * weights_longuest_edgeC[k].weight;
+                        for(int k=0; k< (int)weights_longer_c.size(); k++) {
+                                dp[3] = weights_longer_c[k].value;
+
+                                sum += weights_short_a[i].weight * weights_long_b[j].weight
+                                        * weights_longer_c[k].weight * Parallelepiped(dp, q);
+
+                                norm += weights_short_a[i].weight
+                                         * weights_long_b[j].weight * weights_longer_c[k].weight;
                         }
                 }
@@ -116 +116 @@
         // Fill parameter array
         dp.scale      = scale();
-        dp.short_edgeA   = short_edgeA();
-        dp.longer_edgeB   = longer_edgeB();
-        dp.longuest_edgeC  = longuest_edgeC();
+        dp.short_a   = short_a();
+        dp.long_b   = long_b();
+        dp.longer_c  = longer_c();
         dp.contrast   = contrast();
         dp.background = 0.0;
@@ -128 +128 @@
 
         // Get the dispersion points for the short_edgeA
-        vector<WeightPoint> weights_short_edgeA;
-        short_edgeA.get_weights(weights_short_edgeA);
+        vector<WeightPoint> weights_short_a;
+        short_a.get_weights(weights_short_a);
 
         // Get the dispersion points for the longer_edgeB
-        vector<WeightPoint> weights_longer_edgeB;
-        longer_edgeB.get_weights(weights_longer_edgeB);
+        vector<WeightPoint> weights_long_b;
+        long_b.get_weights(weights_long_b);
 
         // Get angular averaging for the longuest_edgeC
-        vector<WeightPoint> weights_longuest_edgeC;
-        longuest_edgeC.get_weights(weights_longuest_edgeC);
+        vector<WeightPoint> weights_longer_c;
+        longer_c.get_weights(weights_longer_c);
 
         // Get angular averaging for theta
@@ -156 +156 @@
 
         // Loop over radius weight points
-        for(int i=0; i< (int)weights_short_edgeA.size(); i++) {
-                dp.short_edgeA = weights_short_edgeA[i].value;
+        for(int i=0; i< (int)weights_short_a.size(); i++) {
+                dp.short_a = weights_short_a[i].value;
 
                 // Loop over longer_edgeB weight points
-                for(int j=0; j< (int)weights_longer_edgeB.size(); j++) {
-                        dp.longer_edgeB = weights_longer_edgeB[j].value;
+                for(int j=0; j< (int)weights_long_b.size(); j++) {
+                        dp.long_b = weights_long_b[j].value;
 
                         // Average over longuest_edgeC distribution
-                        for(int k=0; k< (int)weights_longuest_edgeC.size(); k++) {
-                                dp.longuest_edgeC = weights_longuest_edgeC[k].value;
+                        for(int k=0; k< (int)weights_longer_c.size(); k++) {
+                                dp.longer_c = weights_longer_c[k].value;
 
                                 // Average over theta distribution
@@ -179 +179 @@
                                                         dp.parallel_psi = weights_parallel_psi[n].value;
 
-                                                        double _ptvalue = weights_short_edgeA[i].weight
-                                                                * weights_longer_edgeB[j].weight
-                                                                * weights_longuest_edgeC[k].weight
+                                                        double _ptvalue = weights_short_a[i].weight
+                                                                * weights_long_b[j].weight
+                                                                * weights_longer_c[k].weight
                                                                 * weights_parallel_theta[l].weight
                                                                 * weights_parallel_phi[m].weight
@@ -191 +191 @@
                                                         sum += _ptvalue;
 
-                                                        norm += weights_short_edgeA[i].weight
-                                                                * weights_longer_edgeB[j].weight
-                                                                * weights_longuest_edgeC[k].weight
+                                                        norm += weights_short_a[i].weight
+                                                                * weights_long_b[j].weight
+                                                                * weights_longer_c[k].weight
                                                                 * weights_parallel_theta[l].weight
                                                                 * weights_parallel_phi[m].weight

sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp

@@ -34 +34 @@
 TriaxialEllipsoidModel :: TriaxialEllipsoidModel() {
         scale      = Parameter(1.0);
-        semi_axisA     = Parameter(20.0, true);
+        semi_axisA     = Parameter(35.0, true);
         semi_axisA.set_min(0.0);
-        semi_axisB     = Parameter(20.0, true);
+        semi_axisB     = Parameter(100.0, true);
         semi_axisB.set_min(0.0);
         semi_axisC  = Parameter(400.0, true);
@@ -42 +42 @@
         contrast   = Parameter(5.3e-6);
         background = Parameter(0.0);
-        axis_theta  = Parameter(0.0, true);
-        axis_phi    = Parameter(0.0, true);
+        axis_theta  = Parameter(1.0, true);
+        axis_phi    = Parameter(1.0, true);
         axis_psi    = Parameter(0.0, true);
 }
@@ -190 +190 @@
                                                                 * weights_theta[l].weight
                                                                 * weights_phi[m].weight
-                                                                * weights_psi[m].weight;
+                                                                * weights_psi[n].weight;
                                                 }
                                         }