1 | /** |
---|
2 | Computes 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 | */ |
---|
27 | void 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 | */ |
---|
159 | void 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 | } |
---|