Changeset 9624cda in sasview for sansmodels/src/c_gen/sld2i.cpp

Timestamp:

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

File:

• sansmodels/src/c_gen/sld2i.cpp (modified) (9 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]);
+}