source: sasview/src/sas/sascalc/calculator/c_extensions/sld2i.cpp @ 8548d739

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 8548d739 was b523c0e, checked in by Paul Kienzle <pkienzle@…>, 9 years ago

remove compiler warnings
  • Property mode set to 100644
File size: 5.9 KB
Line 
1/**
2Computes the (magnetic) scattering form sld (n and m) profile
3 */
4#include "sld2i.hh"
5#include <stdio.h>
6#include <math.h>
7using namespace std;
8extern "C" {
9        #include "libfunc.h"
10        #include "librefl.h"
11}
12/**
13 * Constructor for GenI
14 *
15 * binning
16 * //@param qx: array of Qx values
17 * //@param qy: array of Qy values
18 * //@param qz: array of Qz values
19 * @param x: array of x values
20 * @param y: array of y values
21 * @param z: array of z values
22 * @param sldn: array of sld n
23 * @param mx: array of sld mx
24 * @param my: array of sld my
25 * @param mz: array of sld mz
26 * @param in_spin: ratio of up spin in Iin
27 * @param out_spin: ratio of up spin in Iout
28 * @param s_theta: angle (from x-axis) of the up spin in degree
29 */
30GenI :: GenI(int npix, double* x, double* y, double* z, double* sldn,
31                        double* mx, double* my, double* mz, double* voli,
32                        double in_spin, double out_spin,
33                        double s_theta) {
34        //this->qx_val = qx;
35        //this->qy_val = qy;
36        //this->qz_val = qz;
37        this->n_pix = npix;
38        this->x_val = x;
39        this->y_val = y;
40        this->z_val = z;
41        this->sldn_val = sldn;
42        this->mx_val = mx;
43        this->my_val = my;
44        this->mz_val = mz;
45        this->vol_pix = voli;
46        this->inspin = in_spin;
47        this->outspin = out_spin;
48        this->stheta = s_theta;
49};
50
51/**
52 * Compute 2D anisotropic
53 */
54void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){
55        //npoints is given negative for angular averaging
56        // Assumes that q doesn't have qz component and sld_n is all real
57        //double q = 0.0;
58        //double Pi = 4.0*atan(1.0);
59        polar_sld b_sld;
60        double qr = 0.0;
61        complex iqr = cassign(0.0, 0.0);
62        complex ephase = cassign(0.0, 0.0);
63        complex comp_sld = cassign(0.0, 0.0);
64
65        complex sumj_uu;
66        complex sumj_ud;
67        complex sumj_du;
68        complex sumj_dd;
69        complex temp_fi;
70
71        double count = 0.0;
72        //check if this computation is for averaging
73
74        //Assume that pixel volumes are given in vol_pix in A^3 unit
75        //int x_size = 0; //in Ang
76        //int y_size = 0; //in Ang
77        //int z_size = 0; //in Ang
78       
79        // Loop over q-values and multiply apply matrix
80       
81        for(int i=0; i<npoints; i++){
82                //I_out[i] = 0.0;
83                sumj_uu = cassign(0.0, 0.0);
84                sumj_ud = cassign(0.0, 0.0);
85                sumj_du = cassign(0.0, 0.0);
86                sumj_dd = cassign(0.0, 0.0);           
87                //printf ("%d ", i);
88                //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]);
89
90                for(int j=0; j<n_pix; j++){
91                        if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0)
92                        {       
93                                //anisotropic
94                                temp_fi = cassign(0.0, 0.0);
95                                b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j],
96                                                         mx_val[j], my_val[j], mz_val[j],
97                                                         inspin, outspin, stheta);
98                                qr = (qx[i]*x_val[j] + qy[i]*y_val[j]);
99                                iqr = cassign(0.0, qr);
100                                ephase = cplx_exp(iqr);
101                               
102                                //Let's multiply pixel(atomic) volume here
103                                ephase = rcmult(vol_pix[j], ephase);
104                                //up_up
105                                if (inspin > 0.0 && outspin > 0.0){
106                                        comp_sld = cassign(b_sld.uu, 0.0);
107                                        temp_fi = cplx_mult(comp_sld, ephase);
108                                        sumj_uu = cplx_add(sumj_uu, temp_fi);
109                                }
110                                //down_down
111                                if (inspin < 1.0 && outspin < 1.0){
112                                        comp_sld = cassign(b_sld.dd, 0.0);
113                                        temp_fi = cplx_mult(comp_sld, ephase);
114                                        sumj_dd = cplx_add(sumj_dd, temp_fi);
115                                }
116                                //up_down
117                                if (inspin > 0.0 && outspin < 1.0){
118                                        comp_sld = cassign(b_sld.re_ud, b_sld.im_ud);
119                                        temp_fi = cplx_mult(comp_sld, ephase);
120                                        sumj_ud = cplx_add(sumj_ud, temp_fi);
121                                }
122                                //down_up
123                                if (inspin < 1.0 && outspin > 0.0){
124                                        comp_sld = cassign(b_sld.re_du, b_sld.im_du);
125                                        temp_fi = cplx_mult(comp_sld, ephase);
126                                        sumj_du = cplx_add(sumj_du, temp_fi);
127                                }
128
129
130                                if (i == 0){
131                                        count += vol_pix[j];
132                                }
133                        }
134                }
135                //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);
136
137                I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im);
138                I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im);
139                I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im);
140                I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im);
141
142                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
143        }
144        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
145}
146/**
147 * Compute 1D isotropic
148 * Isotropic: Assumes all slds are real (no magnetic)
149 * Also assumes there is no polarization: No dependency on spin
150 */
151void GenI :: genicom(int npoints, double *q, double *I_out){
152        //npoints is given negative for angular averaging
153        // Assumes that q doesn't have qz component and sld_n is all real
154        //double Pi = 4.0*atan(1.0);
155        int is_sym = 0;
156        double qr = 0.0;
157        double sumj;
158        double sld_j = 0.0;
159        double count = 0.0;
160        if (n_pix < 0 ){
161                is_sym = 1;
162                n_pix = n_pix * -1;
163        }
164        //Assume that pixel volumes are given in vol_pix in A^3 unit
165        // Loop over q-values and multiply apply matrix
166        for(int i=0; i<npoints; i++){
167                sumj =0.0;             
168                for(int j=0; j<n_pix; j++){
169                        //Isotropic: Assumes all slds are real (no magnetic)
170                        //Also assumes there is no polarization: No dependency on spin
171                        if (is_sym == 1){
172                                // approximation for a spherical symmetric particle
173                                qr = sqrt(x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j])*q[i];
174                                if (qr > 0.0){
175                                        qr = sin(qr) / qr;
176                                        sumj += sldn_val[j] * vol_pix[j] * qr;
177                                }
178                                else{
179                                        sumj += sldn_val[j] * vol_pix[j];
180                                }
181                        }
182                        else{
183                                //full calculation
184                                //pragma omp parallel for
185                                for(int k=0; k<n_pix; k++){
186                                        sld_j =  sldn_val[j] * sldn_val[k] * vol_pix[j] * vol_pix[k];
187                                        qr = (x_val[j]-x_val[k])*(x_val[j]-x_val[k])+
188                                                      (y_val[j]-y_val[k])*(y_val[j]-y_val[k])+ 
189                                                      (z_val[j]-z_val[k])*(z_val[j]-z_val[k]);
190                                        qr = sqrt(qr) * q[i];
191                                        if (qr > 0.0){
192                                                sumj += sld_j*sin(qr)/qr;
193                                        }
194                                        else{
195                                                sumj += sld_j;
196                                        }
197                                }
198                        }
199                        if (i == 0){
200                                count += vol_pix[j];
201                        }
202                }
203                I_out[i] = sumj;
204                if (is_sym == 1){
205                        I_out[i] *= sumj;
206                }
207                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
208        }
209        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
210}
Note: See TracBrowser for help on using the repository browser.