Changeset 9624cda in sasview for sansmodels

Timestamp:

Mar 13, 2013 11:10:51 AM (12 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:

49a8843

Parents:

5f248f6

Message:

fixed debye averging in gensans

Location:

sansmodels/src/c_gen

Files:

• sld2i.cpp (modified) (9 diffs)
• sld2i.hh (modified) (1 diff)
• sld2i_module.cpp (modified) (3 diffs)

sansmodels/src/c_gen/sld2i.cpp

@@ -50 +50 @@
 
 /**
- * Compute
- */
-void GenI :: genicom(int npoints, double *qx, double *qy, double *I_out){
+ * Compute 2D anisotropic
+ */
+void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){
         //npoints is given negative for angular averaging 
         // Assumes that q doesn't have qz component and sld_n is all real
-        double q = 0.0;
-        int is_avg = 0;
+        //double q = 0.0;
         //double Pi = 4.0*atan(1.0);
         polar_sld b_sld;
@@ -63 +62 @@
         complex ephase = cassign(0.0, 0.0);
         complex comp_sld = cassign(0.0, 0.0);
-        complex sumj;
+
         complex sumj_uu;
         complex sumj_ud;
@@ -72 +71 @@
         double count = 0.0;
         //check if this computation is for averaging
-        if (n_pix < 0 ){
-                is_avg = 1;
-                n_pix = fabs(n_pix);
-        }
+
         //Assume that pixel volumes are given in vol_pix in A^3 unit
         //int x_size = 0; //in Ang
         //int y_size = 0; //in Ang
         //int z_size = 0; //in Ang
+        
         // Loop over q-values and multiply apply matrix
-
+        
         for(int i=0; i<npoints; i++){
                 //I_out[i] = 0.0;
-                sumj = cassign(0.0, 0.0);
                 sumj_uu = cassign(0.0, 0.0);
                 sumj_ud = cassign(0.0, 0.0);
@@ -90 +86 @@
                 sumj_dd = cassign(0.0, 0.0);            
                 //printf ("%d ", i);
-                q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
+                //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
 
                 for(int j=0; j<n_pix; j++){
-                        
                         if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0)
                         {       
+                                //anisotropic
                                 temp_fi = cassign(0.0, 0.0);
                                 b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j],
-                                                                 mx_val[j], my_val[j], mz_val[j],
-                                                                 inspin, outspin, stheta);
-                                if (is_avg < 1){
-                                        //anisotropic
-                                        qr = (qx[i]*x_val[j] + qy[i]*y_val[j]);
-                                        iqr = cassign(0.0, qr);
-                                        ephase = cplx_exp(iqr);
-                                        }
-                                else{
-                                        //isotropic
-                                        qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j]);
-                                        qr *= q;
-                                        if (qr == 0.0){
-                                                ephase = cassign(0.0, 1.0);
-                                                }
-                                        else{
-                                                qr = sin(qr) / qr;
-                                                ephase = cassign(qr, 0.0);
-                                                }
-                                        }
+                                                         mx_val[j], my_val[j], mz_val[j],
+                                                         inspin, outspin, stheta);
+                                qr = (qx[i]*x_val[j] + qy[i]*y_val[j]);
+                                iqr = cassign(0.0, qr);
+                                ephase = cplx_exp(iqr);
+                                
                                 //Let's multiply pixel(atomic) volume here
                                 ephase = rcmult(vol_pix[j], ephase);
@@ -125 +107 @@
                                         temp_fi = cplx_mult(comp_sld, ephase);
                                         sumj_uu = cplx_add(sumj_uu, temp_fi);
-                                        }
+                                }
                                 //down_down
                                 if (inspin < 1.0 && outspin < 1.0){
@@ -131 +113 @@
                                         temp_fi = cplx_mult(comp_sld, ephase);
                                         sumj_dd = cplx_add(sumj_dd, temp_fi);
-                                        }
+                                }
                                 //up_down
                                 if (inspin > 0.0 && outspin < 1.0){
@@ -137 +119 @@
                                         temp_fi = cplx_mult(comp_sld, ephase);
                                         sumj_ud = cplx_add(sumj_ud, temp_fi);
-                                        }
+                                }
                                 //down_up
                                 if (inspin < 1.0 && outspin > 0.0){
@@ -143 +125 @@
                                         temp_fi = cplx_mult(comp_sld, ephase);
                                         sumj_du = cplx_add(sumj_du, temp_fi);
-                                        }
+                                }
+
+
                                 if (i == 0){
                                         count += vol_pix[j];
@@ -150 +134 @@
                 }
                 //printf("aa%d=%g %g %d\n", i, (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im), (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im), count);
+
                 I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im);
                 I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im);
                 I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im);
                 I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im);
+
                 I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
         }
         //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
 }
+/**
+ * Compute 1D isotropic
+ * Isotropic: Assumes all slds are real (no magnetic)
+ * Also assumes there is no polarization: No dependency on spin
+ */
+void GenI :: genicom(int npoints, double *q, double *I_out){
+        //npoints is given negative for angular averaging 
+        // Assumes that q doesn't have qz component and sld_n is all real
+        //double Pi = 4.0*atan(1.0);
+        int is_sym = 0;
+        double qr = 0.0;
+        double sumj;
+        double sld_j = 0.0;
+        double count = 0.0;
+        if (n_pix < 0 ){
+                is_sym = 1;
+                n_pix = n_pix * -1;
+        }
+        //Assume that pixel volumes are given in vol_pix in A^3 unit
+        // Loop over q-values and multiply apply matrix
+        for(int i=0; i<npoints; i++){
+                sumj =0.0;              
+                for(int j=0; j<n_pix; j++){
+                        //Isotropic: Assumes all slds are real (no magnetic)
+                        //Also assumes there is no polarization: No dependency on spin
+                        if (is_sym == 1){
+                                // approximation for a spherical symmetric particle
+                                qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j])*q[i];
+                                if (qr > 0.0){
+                                        qr = sin(qr) / qr;
+                                        sumj += sldn_val[j] * vol_pix[j] * qr;
+                                }
+                                else{
+                                        sumj += sldn_val[j] * vol_pix[j];
+                                }
+                        }
+                        else{
+                                //full calculation
+                                #pragma omp parallel for
+                                for(int k=0; k<n_pix; k++){
+                                        sld_j =  sldn_val[j] * sldn_val[k] * vol_pix[j] * vol_pix[k];
+                                        qr = (x_val[j]-x_val[k])*(x_val[j]-x_val[k])+
+                                                      (y_val[j]-y_val[k])*(y_val[j]-y_val[k])+ 
+                                                      (z_val[j]-z_val[k])*(z_val[j]-z_val[k]);
+                                        qr = sqrt(qr) * q[i];
+                                        if (qr > 0.0){
+                                                sumj += sld_j*sin(qr)/qr;
+                                        }
+                                        else{
+                                                sumj += sld_j;
+                                        }
+                                }
+                        }
+                        if (i == 0){
+                                count += vol_pix[j];
+                        }
+                }
+                I_out[i] = sumj;
+                if (is_sym == 1){
+                        I_out[i] *= sumj;
+                }
+                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
+        }
+        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
+}

sansmodels/src/c_gen/sld2i.hh

@@ -39 +39 @@
                         double s_theta);
         // compute function
-        virtual void genicom(int npoints, double* qx, double* qy, double *I_out);
+        virtual void genicomXY(int npoints, double* qx, double* qy, double *I_out);
+        virtual void genicom(int npoints, double* q, double *I_out);
 };
 

sansmodels/src/c_gen/sld2i_module.cpp

@@ -75 +75 @@
 
 /**
- * GenI the given input according to a given object
+ * GenI the given input (2D) according to a given object
  */
-PyObject * genicom_input(PyObject *, PyObject *args) {
+PyObject * genicom_inputXY(PyObject *, PyObject *args) {
         int npoints;
         PyObject *qx_obj;
@@ -100 +100 @@
         GenI* s = static_cast<GenI *>(temp);
 
-        s->genicom(npoints, qx, qy, I_out);
+        s->genicomXY(npoints, qx, qy, I_out);
         //return PyCObject_FromVoidPtr(s, del_genicom);
         return Py_BuildValue("i",1);
 }
 
+/**
+ * GenI the given 1D input according to a given object
+ */
+PyObject * genicom_input(PyObject *, PyObject *args) {
+        int npoints;
+        PyObject *q_obj;
+        double *q;
+        PyObject *I_out_obj;
+        Py_ssize_t n_out;
+        double *I_out;
+        PyObject *gen_obj;
+
+        if (!PyArg_ParseTuple(args, "OiOO",  &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL;
+        OUTVECTOR(q_obj, q, n_out);
+        OUTVECTOR(I_out_obj, I_out, n_out);
+
+        // Sanity check
+        //if(n_in!=n_out) return Py_BuildValue("i",-1);
+
+        // Set the array pointers
+        void *temp = PyCObject_AsVoidPtr(gen_obj);
+        GenI* s = static_cast<GenI *>(temp);
+
+        s->genicom(npoints, q, I_out);
+        //return PyCObject_FromVoidPtr(s, del_genicom);
+        return Py_BuildValue("i",1);
+}
 
 /**
@@ -113 +140 @@
                   "Create a new GenI object"},
         {"genicom",(PyCFunction)genicom_input, METH_VARARGS,
-                  "genicom the given input arrays"},
+                  "genicom the given 1d input arrays"},
+        {"genicomXY",(PyCFunction)genicom_inputXY, METH_VARARGS,
+                  "genicomXY the given 2d input arrays"},
     {NULL}
 };