| 1 | /** |
|---|
| 2 | Computes the (magnetic) scattering form sld (n and m) profile |
|---|
| 3 | */ |
|---|
| 4 | #include "sld2i.hh" |
|---|
| 5 | #include <stdio.h> |
|---|
| 6 | #include <math.h> |
|---|
| 7 | using namespace std; |
|---|
| 8 | extern "C" { |
|---|
| 9 | #include "libmultifunc/libfunc.h" |
|---|
| 10 | #include "libmultifunc/librefl.h" |
|---|
| 11 | } |
|---|
| 12 | /** |
|---|
| 13 | * Constructor for GenI |
|---|
| 14 | * |
|---|
| 15 | * binning |
|---|
| 16 | * //@param qx: array of Qx values |
|---|
| 17 | * //@param qy: array of Qy values |
|---|
| 18 | * //@param qz: array of Qz values |
|---|
| 19 | * @param x: array of x values |
|---|
| 20 | * @param y: array of y values |
|---|
| 21 | * @param z: array of z values |
|---|
| 22 | * @param sldn: array of sld n |
|---|
| 23 | * @param mx: array of sld mx |
|---|
| 24 | * @param my: array of sld my |
|---|
| 25 | * @param mz: array of sld mz |
|---|
| 26 | * @param in_spin: ratio of up spin in Iin |
|---|
| 27 | * @param out_spin: ratio of up spin in Iout |
|---|
| 28 | * @param s_theta: angle (from x-axis) of the up spin in degree |
|---|
| 29 | */ |
|---|
| 30 | GenI :: GenI(int npix, double* x, double* y, double* z, double* sldn, |
|---|
| 31 | double* mx, double* my, double* mz, double* voli, |
|---|
| 32 | double in_spin, double out_spin, |
|---|
| 33 | double s_theta) { |
|---|
| 34 | //this->qx_val = qx; |
|---|
| 35 | //this->qy_val = qy; |
|---|
| 36 | //this->qz_val = qz; |
|---|
| 37 | this->n_pix = npix; |
|---|
| 38 | this->x_val = x; |
|---|
| 39 | this->y_val = y; |
|---|
| 40 | this->z_val = z; |
|---|
| 41 | this->sldn_val = sldn; |
|---|
| 42 | this->mx_val = mx; |
|---|
| 43 | this->my_val = my; |
|---|
| 44 | this->mz_val = mz; |
|---|
| 45 | this->vol_pix = voli; |
|---|
| 46 | this->inspin = in_spin; |
|---|
| 47 | this->outspin = out_spin; |
|---|
| 48 | this->stheta = s_theta; |
|---|
| 49 | }; |
|---|
| 50 | |
|---|
| 51 | /** |
|---|
| 52 | * Compute |
|---|
| 53 | */ |
|---|
| 54 | void GenI :: genicom(int npoints, double *qx, double *qy, double *I_out){ |
|---|
| 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 | complex iqr = cassign(0.0, 0.0); |
|---|
| 61 | complex ephase = cassign(0.0, 0.0); |
|---|
| 62 | complex comp_sld = cassign(0.0, 0.0); |
|---|
| 63 | complex sumj; |
|---|
| 64 | complex sumj_uu; |
|---|
| 65 | complex sumj_ud; |
|---|
| 66 | complex sumj_du; |
|---|
| 67 | complex sumj_dd; |
|---|
| 68 | complex temp_fi; |
|---|
| 69 | |
|---|
| 70 | double count = 0.0; |
|---|
| 71 | //Assume that pixel volumes are given in vol_pix in A^3 unit |
|---|
| 72 | //int x_size = 0; //in Ang |
|---|
| 73 | //int y_size = 0; //in Ang |
|---|
| 74 | //int z_size = 0; //in Ang |
|---|
| 75 | // Loop over q-values and multiply apply matrix |
|---|
| 76 | |
|---|
| 77 | for(int i=0; i<npoints; i++){ |
|---|
| 78 | //I_out[i] = 0.0; |
|---|
| 79 | sumj = cassign(0.0, 0.0); |
|---|
| 80 | sumj_uu = cassign(0.0, 0.0); |
|---|
| 81 | sumj_ud = cassign(0.0, 0.0); |
|---|
| 82 | sumj_du = cassign(0.0, 0.0); |
|---|
| 83 | sumj_dd = cassign(0.0, 0.0); |
|---|
| 84 | //printf ("%d ", i); |
|---|
| 85 | q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); |
|---|
| 86 | |
|---|
| 87 | for(int j=0; j<n_pix; j++){ |
|---|
| 88 | |
|---|
| 89 | if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0) |
|---|
| 90 | { |
|---|
| 91 | temp_fi = cassign(0.0, 0.0); |
|---|
| 92 | b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j], |
|---|
| 93 | mx_val[j], my_val[j], mz_val[j], |
|---|
| 94 | inspin, outspin, stheta); |
|---|
| 95 | qr = (qx[i]*x_val[j] + qy[i]*y_val[j]); |
|---|
| 96 | iqr = cassign(0.0, qr); |
|---|
| 97 | ephase = cplx_exp(iqr); |
|---|
| 98 | //up_up |
|---|
| 99 | if (inspin > 0.0 && outspin > 0.0){ |
|---|
| 100 | comp_sld = cassign(b_sld.uu, 0.0); |
|---|
| 101 | temp_fi = cplx_mult(comp_sld, ephase); |
|---|
| 102 | temp_fi = rcmult(vol_pix[j],temp_fi); |
|---|
| 103 | sumj_uu = cplx_add(sumj_uu, temp_fi); |
|---|
| 104 | } |
|---|
| 105 | //down_down |
|---|
| 106 | if (inspin < 1.0 && outspin < 1.0){ |
|---|
| 107 | comp_sld = cassign(b_sld.dd, 0.0); |
|---|
| 108 | temp_fi = cplx_mult(comp_sld, ephase); |
|---|
| 109 | temp_fi = rcmult(vol_pix[j],temp_fi); |
|---|
| 110 | sumj_dd = cplx_add(sumj_dd, temp_fi); |
|---|
| 111 | } |
|---|
| 112 | //up_down |
|---|
| 113 | if (inspin > 0.0 && outspin < 1.0){ |
|---|
| 114 | comp_sld = cassign(b_sld.re_ud, b_sld.im_ud); |
|---|
| 115 | temp_fi = cplx_mult(comp_sld, ephase); |
|---|
| 116 | temp_fi = rcmult(vol_pix[j],temp_fi); |
|---|
| 117 | sumj_ud = cplx_add(sumj_ud, temp_fi); |
|---|
| 118 | } |
|---|
| 119 | //down_up |
|---|
| 120 | if (inspin < 1.0 && outspin > 0.0){ |
|---|
| 121 | comp_sld = cassign(b_sld.re_du, b_sld.im_du); |
|---|
| 122 | temp_fi = cplx_mult(comp_sld, ephase); |
|---|
| 123 | temp_fi = rcmult(vol_pix[j],temp_fi); |
|---|
| 124 | sumj_du = cplx_add(sumj_du, temp_fi); |
|---|
| 125 | } |
|---|
| 126 | if (i == 0){ |
|---|
| 127 | count += vol_pix[j]; |
|---|
| 128 | } |
|---|
| 129 | } |
|---|
| 130 | } |
|---|
| 131 | //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); |
|---|
| 132 | I_out[i] = (sumj_uu.re*sumj_uu.re + sumj_uu.im*sumj_uu.im); |
|---|
| 133 | I_out[i] += (sumj_ud.re*sumj_ud.re + sumj_ud.im*sumj_ud.im); |
|---|
| 134 | I_out[i] += (sumj_du.re*sumj_du.re + sumj_du.im*sumj_du.im); |
|---|
| 135 | I_out[i] += (sumj_dd.re*sumj_dd.re + sumj_dd.im*sumj_dd.im); |
|---|
| 136 | I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix |
|---|
| 137 | } |
|---|
| 138 | //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); |
|---|
| 139 | } |
|---|
