source: sasview/src/sas/sascalc/calculator/c_extensions/sld2i.c @ f00d3fd

ESS_GUIESS_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 f00d3fd was e6f2009, checked in by trnielsen, 6 years ago

C99 compliant

Changes to be committed:

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