source: sasview/sansmodels/src/c_gen/sld2i.cpp @ 2d9c7266

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 2d9c7266 was 318b5bbb, checked in by Jae Cho <jhjcho@…>, 12 years ago

Added polarization and magnetic stuffs

  • Property mode set to 100644
File size: 3.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 "libmultifunc/libfunc.h"
10        #include "libmultifunc/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,
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->inspin = in_spin;
46        this->outspin = out_spin;
47        this->stheta = s_theta;
48};
49
50/**
51 * Compute
52 */
53void GenI :: genicom(int npoints, double *qx, double *qy, double *I_out){
54        double q = 0.0;
55        //double Pi = 4.0*atan(1.0);
56        polar_sld b_sld;
57        double qr = 0.0;
58        complex iqr = cassign(0.0, 0.0);
59        complex ephase = cassign(0.0, 0.0);
60        complex comp_sld = cassign(0.0, 0.0);
61        complex sumj;
62        complex sumj_uu;
63        complex sumj_ud;
64        complex sumj_du;
65        complex sumj_dd;
66        complex temp_fi;
67
68        int count = 0;
69        //Assume that all pixels are in same size and int * 1A
70        //int x_size = 0; //in Ang
71        //int y_size = 0; //in Ang
72        //int z_size = 0; //in Ang
73        // Loop over q-values and multiply apply matrix
74
75        for(int i=0; i<npoints; i++){
76                //I_out[i] = 0.0;
77                sumj = cassign(0.0, 0.0);
78                sumj_uu = cassign(0.0, 0.0);
79                sumj_ud = cassign(0.0, 0.0);
80                sumj_du = cassign(0.0, 0.0);
81                sumj_dd = cassign(0.0, 0.0);           
82                //printf ("%d ", i);
83                q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]);
84
85                for(int j=0; j<n_pix; j++){
86                       
87                        if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0)
88                        {       
89                                temp_fi = cassign(0.0, 0.0);
90                                b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j],
91                                                                 mx_val[j], my_val[j], mz_val[j],
92                                                                 inspin, outspin, stheta);
93                                qr = (qx[i]*x_val[j] + qy[i]*y_val[j]);
94                                iqr = cassign(0.0, qr);
95                                ephase = cplx_exp(iqr);
96                                //up_up
97                                if (inspin > 0.0 && outspin > 0.0){
98                                        comp_sld = cassign(b_sld.uu, 0.0);
99                                        temp_fi = cplx_mult(comp_sld, ephase);
100                                        sumj_uu = cplx_add(sumj_uu, temp_fi);
101                                        }
102                                //down_down
103                                if (inspin < 1.0 && outspin < 1.0){
104                                        comp_sld = cassign(b_sld.dd, 0.0);
105                                        temp_fi = cplx_mult(comp_sld, ephase);
106                                        sumj_dd = cplx_add(sumj_dd, temp_fi);
107                                        }
108                                //up_down
109                                if (inspin > 0.0 && outspin < 1.0){
110                                        comp_sld = cassign(b_sld.re_ud, b_sld.im_ud);
111                                        temp_fi = cplx_mult(comp_sld, ephase);
112                                        sumj_ud = cplx_add(sumj_ud, temp_fi);
113                                        }
114                                //down_up
115                                if (inspin < 1.0 && outspin > 0.0){
116                                        comp_sld = cassign(b_sld.re_du, b_sld.im_du);
117                                        temp_fi = cplx_mult(comp_sld, ephase);
118                                        sumj_du = cplx_add(sumj_du, temp_fi);
119                                        }
120                                if (i == 0){
121                                        count += 1;
122                                }
123                        }
124                }
125                //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);
126                I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im);
127                I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im);
128                I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im);
129                I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im);
130                I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix
131        }
132        //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]);
133}
Note: See TracBrowser for help on using the repository browser.