[9e531f2] | 1 | /** |
---|
| 2 | Computes the (magnetic) scattering form sld (n and m) profile |
---|
| 3 | */ |
---|
| 4 | #include <stdio.h> |
---|
| 5 | #include <math.h> |
---|
[3010f68] | 6 | #include "sld2i.h" |
---|
| 7 | #include "libfunc.h" |
---|
| 8 | #include "librefl.h" |
---|
[9e531f2] | 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 | */ |
---|
[3010f68] | 27 | void initGenI(GenI* this, int npix, double* x, double* y, double* z, double* sldn, |
---|
[9e531f2] | 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; |
---|
[3010f68] | 43 | } |
---|
[9e531f2] | 44 | |
---|
| 45 | /** |
---|
| 46 | * Compute 2D anisotropic |
---|
| 47 | */ |
---|
[3010f68] | 48 | void genicomXY(GenI* this, int npoints, double *qx, double *qy, double *I_out){ |
---|
| 49 | //npoints is given negative for angular averaging |
---|
[9e531f2] | 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; |
---|
[f54e82cf] | 55 | Cplx iqr; |
---|
| 56 | Cplx ephase; |
---|
| 57 | Cplx comp_sld; |
---|
[9e531f2] | 58 | |
---|
[f54e82cf] | 59 | Cplx sumj_uu; |
---|
| 60 | Cplx sumj_ud; |
---|
| 61 | Cplx sumj_du; |
---|
| 62 | Cplx sumj_dd; |
---|
| 63 | Cplx temp_fi; |
---|
| 64 | |
---|
| 65 | int i, j; |
---|
[9e531f2] | 66 | |
---|
| 67 | double count = 0.0; |
---|
| 68 | //check if this computation is for averaging |
---|
| 69 | |
---|
[f54e82cf] | 70 | cassign(&iqr, 0.0, 0.0); |
---|
| 71 | cassign(&ephase, 0.0, 0.0); |
---|
| 72 | cassign(&comp_sld, 0.0, 0.0); |
---|
| 73 | |
---|
[9e531f2] | 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 |
---|
[3010f68] | 78 | |
---|
[9e531f2] | 79 | // Loop over q-values and multiply apply matrix |
---|
[e6f2009] | 80 | for(i=0; i<npoints; i++){ |
---|
[9e531f2] | 81 | //I_out[i] = 0.0; |
---|
[f54e82cf] | 82 | cassign(&sumj_uu, 0.0, 0.0); |
---|
| 83 | cassign(&sumj_ud, 0.0, 0.0); |
---|
| 84 | cassign(&sumj_du, 0.0, 0.0); |
---|
| 85 | cassign(&sumj_dd, 0.0, 0.0); |
---|
| 86 | //printf("i: %d\n", i); |
---|
[9e531f2] | 87 | //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); |
---|
[e6f2009] | 88 | for(j=0; j<this->n_pix; j++){ |
---|
[3010f68] | 89 | if (this->sldn_val[j]!=0.0 |
---|
| 90 | ||this->mx_val[j]!=0.0 |
---|
| 91 | ||this->my_val[j]!=0.0 |
---|
| 92 | ||this->mz_val[j]!=0.0) |
---|
| 93 | { |
---|
[f54e82cf] | 94 | // printf("i,j: %d,%d\n", i,j); |
---|
[9e531f2] | 95 | //anisotropic |
---|
[f54e82cf] | 96 | cassign(&temp_fi, 0.0, 0.0); |
---|
| 97 | cal_msld(&b_sld, 0, qx[i], qy[i], this->sldn_val[j], |
---|
[3010f68] | 98 | this->mx_val[j], this->my_val[j], this->mz_val[j], |
---|
| 99 | this->inspin, this->outspin, this->stheta); |
---|
| 100 | qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); |
---|
[f54e82cf] | 101 | cassign(&iqr, 0.0, qr); |
---|
| 102 | cplx_exp(&ephase, iqr); |
---|
[3010f68] | 103 | |
---|
[9e531f2] | 104 | //Let's multiply pixel(atomic) volume here |
---|
[f54e82cf] | 105 | rcmult(&ephase, this->vol_pix[j], ephase); |
---|
[9e531f2] | 106 | //up_up |
---|
[3010f68] | 107 | if (this->inspin > 0.0 && this->outspin > 0.0){ |
---|
[f54e82cf] | 108 | cassign(&comp_sld, b_sld.uu, 0.0); |
---|
| 109 | cplx_mult(&temp_fi, comp_sld, ephase); |
---|
| 110 | cplx_add(&sumj_uu, sumj_uu, temp_fi); |
---|
[9e531f2] | 111 | } |
---|
| 112 | //down_down |
---|
[3010f68] | 113 | if (this->inspin < 1.0 && this->outspin < 1.0){ |
---|
[f54e82cf] | 114 | cassign(&comp_sld, b_sld.dd, 0.0); |
---|
| 115 | cplx_mult(&temp_fi, comp_sld, ephase); |
---|
| 116 | cplx_add(&sumj_dd, sumj_dd, temp_fi); |
---|
[9e531f2] | 117 | } |
---|
| 118 | //up_down |
---|
[3010f68] | 119 | if (this->inspin > 0.0 && this->outspin < 1.0){ |
---|
[f54e82cf] | 120 | cassign(&comp_sld, b_sld.re_ud, b_sld.im_ud); |
---|
| 121 | cplx_mult(&temp_fi, comp_sld, ephase); |
---|
| 122 | cplx_add(&sumj_ud, sumj_ud, temp_fi); |
---|
[9e531f2] | 123 | } |
---|
| 124 | //down_up |
---|
[3010f68] | 125 | if (this->inspin < 1.0 && this->outspin > 0.0){ |
---|
[f54e82cf] | 126 | cassign(&comp_sld, b_sld.re_du, b_sld.im_du); |
---|
| 127 | cplx_mult(&temp_fi, comp_sld, ephase); |
---|
| 128 | cplx_add(&sumj_du, sumj_du, temp_fi); |
---|
[9e531f2] | 129 | } |
---|
| 130 | |
---|
| 131 | if (i == 0){ |
---|
[3010f68] | 132 | count += this->vol_pix[j]; |
---|
[9e531f2] | 133 | } |
---|
| 134 | } |
---|
| 135 | } |
---|
| 136 | //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); |
---|
| 137 | |
---|
| 138 | I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im); |
---|
| 139 | I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im); |
---|
| 140 | I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im); |
---|
| 141 | I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im); |
---|
| 142 | |
---|
| 143 | I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix |
---|
| 144 | } |
---|
| 145 | //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); |
---|
| 146 | } |
---|
| 147 | /** |
---|
| 148 | * Compute 1D isotropic |
---|
| 149 | * Isotropic: Assumes all slds are real (no magnetic) |
---|
| 150 | * Also assumes there is no polarization: No dependency on spin |
---|
| 151 | */ |
---|
[3010f68] | 152 | void genicom(GenI* this, int npoints, double *q, double *I_out){ |
---|
| 153 | //npoints is given negative for angular averaging |
---|
[9e531f2] | 154 | // Assumes that q doesn't have qz component and sld_n is all real |
---|
| 155 | //double Pi = 4.0*atan(1.0); |
---|
[3010f68] | 156 | int is_sym = this->n_pix < 0; |
---|
[9e531f2] | 157 | double qr = 0.0; |
---|
| 158 | double sumj; |
---|
| 159 | double sld_j = 0.0; |
---|
| 160 | double count = 0.0; |
---|
[3010f68] | 161 | int n_pix = is_sym ? -this->n_pix : this->n_pix; |
---|
[9e531f2] | 162 | //Assume that pixel volumes are given in vol_pix in A^3 unit |
---|
| 163 | // Loop over q-values and multiply apply matrix |
---|
[f54e82cf] | 164 | int i, j, k; |
---|
[e6f2009] | 165 | for(i=0; i<npoints; i++){ |
---|
[3010f68] | 166 | sumj =0.0; |
---|
[e6f2009] | 167 | for(j=0; j<n_pix; j++){ |
---|
[9e531f2] | 168 | //Isotropic: Assumes all slds are real (no magnetic) |
---|
| 169 | //Also assumes there is no polarization: No dependency on spin |
---|
| 170 | if (is_sym == 1){ |
---|
| 171 | // approximation for a spherical symmetric particle |
---|
[3010f68] | 172 | 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]; |
---|
[9e531f2] | 173 | if (qr > 0.0){ |
---|
| 174 | qr = sin(qr) / qr; |
---|
[3010f68] | 175 | sumj += this->sldn_val[j] * this->vol_pix[j] * qr; |
---|
[9e531f2] | 176 | } |
---|
| 177 | else{ |
---|
[3010f68] | 178 | sumj += this->sldn_val[j] * this->vol_pix[j]; |
---|
[9e531f2] | 179 | } |
---|
| 180 | } |
---|
| 181 | else{ |
---|
| 182 | //full calculation |
---|
| 183 | //pragma omp parallel for |
---|
[e6f2009] | 184 | for(k=0; k<n_pix; k++){ |
---|
[3010f68] | 185 | sld_j = this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k]; |
---|
| 186 | qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+ |
---|
| 187 | (this->y_val[j]-this->y_val[k])*(this->y_val[j]-this->y_val[k])+ |
---|
| 188 | (this->z_val[j]-this->z_val[k])*(this->z_val[j]-this->z_val[k]); |
---|
[9e531f2] | 189 | qr = sqrt(qr) * q[i]; |
---|
| 190 | if (qr > 0.0){ |
---|
| 191 | sumj += sld_j*sin(qr)/qr; |
---|
| 192 | } |
---|
| 193 | else{ |
---|
| 194 | sumj += sld_j; |
---|
| 195 | } |
---|
| 196 | } |
---|
| 197 | } |
---|
| 198 | if (i == 0){ |
---|
[3010f68] | 199 | count += this->vol_pix[j]; |
---|
[9e531f2] | 200 | } |
---|
| 201 | } |
---|
| 202 | I_out[i] = sumj; |
---|
| 203 | if (is_sym == 1){ |
---|
| 204 | I_out[i] *= sumj; |
---|
| 205 | } |
---|
| 206 | I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix |
---|
| 207 | } |
---|
| 208 | //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); |
---|
| 209 | } |
---|