[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; |
---|
| 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 |
---|
[3010f68] | 72 | |
---|
[9e531f2] | 73 | // Loop over q-values and multiply apply matrix |
---|
[e6f2009] | 74 | int i; |
---|
| 75 | for(i=0; i<npoints; i++){ |
---|
[9e531f2] | 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); |
---|
[3010f68] | 80 | sumj_dd = cassign(0.0, 0.0); |
---|
[9e531f2] | 81 | //printf ("%d ", i); |
---|
| 82 | //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); |
---|
[e6f2009] | 83 | int j; |
---|
| 84 | for(j=0; j<this->n_pix; j++){ |
---|
[3010f68] | 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 | { |
---|
[9e531f2] | 90 | //anisotropic |
---|
| 91 | temp_fi = cassign(0.0, 0.0); |
---|
[3010f68] | 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]); |
---|
[9e531f2] | 96 | iqr = cassign(0.0, qr); |
---|
| 97 | ephase = cplx_exp(iqr); |
---|
[3010f68] | 98 | |
---|
[9e531f2] | 99 | //Let's multiply pixel(atomic) volume here |
---|
[3010f68] | 100 | ephase = rcmult(this->vol_pix[j], ephase); |
---|
[9e531f2] | 101 | //up_up |
---|
[3010f68] | 102 | if (this->inspin > 0.0 && this->outspin > 0.0){ |
---|
[9e531f2] | 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 |
---|
[3010f68] | 108 | if (this->inspin < 1.0 && this->outspin < 1.0){ |
---|
[9e531f2] | 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 |
---|
[3010f68] | 114 | if (this->inspin > 0.0 && this->outspin < 1.0){ |
---|
[9e531f2] | 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 |
---|
[3010f68] | 120 | if (this->inspin < 1.0 && this->outspin > 0.0){ |
---|
[9e531f2] | 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){ |
---|
[3010f68] | 128 | count += this->vol_pix[j]; |
---|
[9e531f2] | 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 | */ |
---|
[3010f68] | 148 | void genicom(GenI* this, int npoints, double *q, double *I_out){ |
---|
| 149 | //npoints is given negative for angular averaging |
---|
[9e531f2] | 150 | // Assumes that q doesn't have qz component and sld_n is all real |
---|
| 151 | //double Pi = 4.0*atan(1.0); |
---|
[3010f68] | 152 | int is_sym = this->n_pix < 0; |
---|
[9e531f2] | 153 | double qr = 0.0; |
---|
| 154 | double sumj; |
---|
| 155 | double sld_j = 0.0; |
---|
| 156 | double count = 0.0; |
---|
[3010f68] | 157 | int n_pix = is_sym ? -this->n_pix : this->n_pix; |
---|
[9e531f2] | 158 | //Assume that pixel volumes are given in vol_pix in A^3 unit |
---|
| 159 | // Loop over q-values and multiply apply matrix |
---|
[e6f2009] | 160 | int i; |
---|
| 161 | for(i=0; i<npoints; i++){ |
---|
[3010f68] | 162 | sumj =0.0; |
---|
[e6f2009] | 163 | int j; |
---|
| 164 | for(j=0; j<n_pix; j++){ |
---|
[9e531f2] | 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 |
---|
[3010f68] | 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]; |
---|
[9e531f2] | 170 | if (qr > 0.0){ |
---|
| 171 | qr = sin(qr) / qr; |
---|
[3010f68] | 172 | sumj += this->sldn_val[j] * this->vol_pix[j] * qr; |
---|
[9e531f2] | 173 | } |
---|
| 174 | else{ |
---|
[3010f68] | 175 | sumj += this->sldn_val[j] * this->vol_pix[j]; |
---|
[9e531f2] | 176 | } |
---|
| 177 | } |
---|
| 178 | else{ |
---|
| 179 | //full calculation |
---|
| 180 | //pragma omp parallel for |
---|
[e6f2009] | 181 | int k; |
---|
| 182 | for(k=0; k<n_pix; k++){ |
---|
[3010f68] | 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]); |
---|
[9e531f2] | 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){ |
---|
[3010f68] | 197 | count += this->vol_pix[j]; |
---|
[9e531f2] | 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 | } |
---|