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

magnetic_scatt
Last change on this file since b36e7c7 was b36e7c7, checked in by dirk, 5 years ago

clean up obsolete? code of magnetic scattering calculation. Addresses #1086.

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