| 1 | // by jcho |
|---|
| 2 | |
|---|
| 3 | #include <math.h> |
|---|
| 4 | |
|---|
| 5 | #include "libfunc.h" |
|---|
| 6 | |
|---|
| 7 | #include <stdio.h> |
|---|
| 8 | |
|---|
| 9 | |
|---|
| 10 | |
|---|
| 11 | //used in Si func |
|---|
| 12 | |
|---|
| 13 | int factorial(int i) { |
|---|
| 14 | |
|---|
| 15 | int k, j; |
|---|
| 16 | if (i<2){ |
|---|
| 17 | return 1; |
|---|
| 18 | } |
|---|
| 19 | |
|---|
| 20 | k=1; |
|---|
| 21 | |
|---|
| 22 | for(j=1;j<i;j++) { |
|---|
| 23 | k=k*(j+1); |
|---|
| 24 | } |
|---|
| 25 | |
|---|
| 26 | return k; |
|---|
| 27 | |
|---|
| 28 | } |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | |
|---|
| 32 | // Used in pearl nec model |
|---|
| 33 | |
|---|
| 34 | // Sine integral function: approximated within 1%!!! |
|---|
| 35 | |
|---|
| 36 | // integral of sin(x)/x up to namx term nmax=6 looks the best. |
|---|
| 37 | |
|---|
| 38 | double Si(double x) |
|---|
| 39 | |
|---|
| 40 | { |
|---|
| 41 | int i; |
|---|
| 42 | int nmax=6; |
|---|
| 43 | double out; |
|---|
| 44 | long double power; |
|---|
| 45 | double pi = 4.0*atan(1.0); |
|---|
| 46 | |
|---|
| 47 | if (x >= pi*6.2/4.0){ |
|---|
| 48 | double out_sin = 0.0; |
|---|
| 49 | double out_cos = 0.0; |
|---|
| 50 | out = pi/2.0; |
|---|
| 51 | |
|---|
| 52 | for (i=0; i<nmax-2; i+=1) { |
|---|
| 53 | out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1); |
|---|
| 54 | out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2); |
|---|
| 55 | } |
|---|
| 56 | |
|---|
| 57 | out -= cos(x) * out_cos; |
|---|
| 58 | out -= sin(x) * out_sin; |
|---|
| 59 | return out; |
|---|
| 60 | } |
|---|
| 61 | |
|---|
| 62 | out = 0.0; |
|---|
| 63 | |
|---|
| 64 | for (i=0; i<nmax; i+=1) { |
|---|
| 65 | if (i==0) { |
|---|
| 66 | out += x; |
|---|
| 67 | continue; |
|---|
| 68 | } |
|---|
| 69 | |
|---|
| 70 | power = pow(x,(2 * i + 1)); |
|---|
| 71 | out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1)); |
|---|
| 72 | |
|---|
| 73 | //printf ("Si=%g %g %d\n", x, out, i); |
|---|
| 74 | } |
|---|
| 75 | |
|---|
| 76 | return out; |
|---|
| 77 | } |
|---|
| 78 | |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | double sinc(double x) |
|---|
| 82 | { |
|---|
| 83 | if (x==0.0){ |
|---|
| 84 | return 1.0; |
|---|
| 85 | } |
|---|
| 86 | return sin(x)/x; |
|---|
| 87 | } |
|---|
| 88 | |
|---|
| 89 | |
|---|
| 90 | double gamln(double xx) { |
|---|
| 91 | |
|---|
| 92 | double x,y,tmp,ser; |
|---|
| 93 | static double cof[6]={76.18009172947146,-86.50532032941677, |
|---|
| 94 | 24.01409824083091,-1.231739572450155, |
|---|
| 95 | 0.1208650973866179e-2,-0.5395239384953e-5}; |
|---|
| 96 | int j; |
|---|
| 97 | |
|---|
| 98 | y=x=xx; |
|---|
| 99 | tmp=x+5.5; |
|---|
| 100 | tmp -= (x+0.5)*log(tmp); |
|---|
| 101 | ser=1.000000000190015; |
|---|
| 102 | for (j=0;j<=5;j++) ser += cof[j]/++y; |
|---|
| 103 | return -tmp+log(2.5066282746310005*ser/x); |
|---|
| 104 | } |
|---|
| 105 | |
|---|
| 106 | // calculate magnetic sld and return total sld |
|---|
| 107 | // bn : contrast (not just sld of the layer) |
|---|
| 108 | // m0: max mag of M; mtheta: angle from x-z plane; |
|---|
| 109 | // mphi: angle (anti-clock-wise)of x-z projection(M) from x axis |
|---|
| 110 | // spinfraci: the fraction of UP among UP+Down (before sample) |
|---|
| 111 | // spinfracf: the fraction of UP among UP+Down (after sample and before detector) |
|---|
| 112 | // spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis |
|---|
| 113 | // Note: all angles are in degrees. |
|---|
| 114 | polar_sld cal_msld(int isangle, double qx, double qy, double bn, |
|---|
| 115 | double m01, double mtheta1, double mphi1, |
|---|
| 116 | double spinfraci, double spinfracf, double spintheta) |
|---|
| 117 | { |
|---|
| 118 | //locals |
|---|
| 119 | double q_x = qx; |
|---|
| 120 | double q_y = qy; |
|---|
| 121 | double sld = bn; |
|---|
| 122 | int is_angle = isangle; |
|---|
| 123 | double pi = 4.0*atan(1.0); |
|---|
| 124 | double s_theta = spintheta * pi/180.0; |
|---|
| 125 | double m_max = m01; |
|---|
| 126 | double m_phi = mphi1; |
|---|
| 127 | double m_theta = mtheta1; |
|---|
| 128 | double in_spin = spinfraci; |
|---|
| 129 | double out_spin = spinfracf; |
|---|
| 130 | |
|---|
| 131 | double m_perp = 0.0; |
|---|
| 132 | double m_perp_z = 0.0; |
|---|
| 133 | double m_perp_y = 0.0; |
|---|
| 134 | double m_perp_x = 0.0; |
|---|
| 135 | double m_sigma_x = 0.0; |
|---|
| 136 | double m_sigma_z = 0.0; |
|---|
| 137 | double m_sigma_y = 0.0; |
|---|
| 138 | //double b_m = 0.0; |
|---|
| 139 | double q_angle = 0.0; |
|---|
| 140 | double mx = 0.0; |
|---|
| 141 | double my = 0.0; |
|---|
| 142 | double mz = 0.0; |
|---|
| 143 | polar_sld p_sld; |
|---|
| 144 | p_sld.uu = sld; |
|---|
| 145 | p_sld.dd = sld; |
|---|
| 146 | p_sld.re_ud = 0.0; |
|---|
| 147 | p_sld.im_ud = 0.0; |
|---|
| 148 | p_sld.re_du = 0.0; |
|---|
| 149 | p_sld.im_du = 0.0; |
|---|
| 150 | |
|---|
| 151 | //No mag means no further calculation |
|---|
| 152 | if (isangle>0){ |
|---|
| 153 | if (m_max < 1.0e-32){ |
|---|
| 154 | p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; |
|---|
| 155 | p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; |
|---|
| 156 | return p_sld; |
|---|
| 157 | } |
|---|
| 158 | } |
|---|
| 159 | else{ |
|---|
| 160 | if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){ |
|---|
| 161 | p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; |
|---|
| 162 | p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; |
|---|
| 163 | return p_sld; |
|---|
| 164 | } |
|---|
| 165 | } |
|---|
| 166 | |
|---|
| 167 | //These are needed because of the precision of inputs |
|---|
| 168 | if (in_spin < 0.0) in_spin = 0.0; |
|---|
| 169 | if (in_spin > 1.0) in_spin = 1.0; |
|---|
| 170 | if (out_spin < 0.0) out_spin = 0.0; |
|---|
| 171 | if (out_spin > 1.0) out_spin = 1.0; |
|---|
| 172 | |
|---|
| 173 | if (q_x == 0.0) q_angle = pi / 2.0; |
|---|
| 174 | else q_angle = atan(q_y/q_x); |
|---|
| 175 | if (q_y < 0.0 && q_x < 0.0) q_angle -= pi; |
|---|
| 176 | else if (q_y > 0.0 && q_x < 0.0) q_angle += pi; |
|---|
| 177 | |
|---|
| 178 | q_angle = pi/2.0 - q_angle; |
|---|
| 179 | if (q_angle > pi) q_angle -= 2.0 * pi; |
|---|
| 180 | else if (q_angle < -pi) q_angle += 2.0 * pi; |
|---|
| 181 | |
|---|
| 182 | if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){ |
|---|
| 183 | m_perp = 0.0; |
|---|
| 184 | } |
|---|
| 185 | else { |
|---|
| 186 | m_perp = m_max; |
|---|
| 187 | } |
|---|
| 188 | if (is_angle > 0){ |
|---|
| 189 | m_phi *= pi/180.0; |
|---|
| 190 | m_theta *= pi/180.0; |
|---|
| 191 | mx = m_perp * cos(m_theta) * cos(m_phi); |
|---|
| 192 | my = m_perp * sin(m_theta); |
|---|
| 193 | mz = -(m_perp * cos(m_theta) * sin(m_phi)); |
|---|
| 194 | } |
|---|
| 195 | else{ |
|---|
| 196 | mx = m_perp; |
|---|
| 197 | my = m_phi; |
|---|
| 198 | mz = m_theta; |
|---|
| 199 | } |
|---|
| 200 | //ToDo: simplify these steps |
|---|
| 201 | // m_perp1 -m_perp2 |
|---|
| 202 | m_perp_x = (mx) * cos(q_angle); |
|---|
| 203 | m_perp_x -= (my) * sin(q_angle); |
|---|
| 204 | m_perp_y = m_perp_x; |
|---|
| 205 | m_perp_x *= cos(-q_angle); |
|---|
| 206 | m_perp_y *= sin(-q_angle); |
|---|
| 207 | m_perp_z = mz; |
|---|
| 208 | |
|---|
| 209 | m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta)); |
|---|
| 210 | m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta)); |
|---|
| 211 | m_sigma_z = (m_perp_z); |
|---|
| 212 | |
|---|
| 213 | //Find b |
|---|
| 214 | p_sld.uu -= m_sigma_x; |
|---|
| 215 | p_sld.dd += m_sigma_x; |
|---|
| 216 | p_sld.re_ud = m_sigma_y; |
|---|
| 217 | p_sld.re_du = m_sigma_y; |
|---|
| 218 | p_sld.im_ud = m_sigma_z; |
|---|
| 219 | p_sld.im_du = -m_sigma_z; |
|---|
| 220 | |
|---|
| 221 | p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; |
|---|
| 222 | p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; |
|---|
| 223 | |
|---|
| 224 | p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud; |
|---|
| 225 | p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud; |
|---|
| 226 | p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du; |
|---|
| 227 | p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du; |
|---|
| 228 | |
|---|
| 229 | return p_sld; |
|---|
| 230 | } |
|---|
| 231 | |
|---|
| 232 | |
|---|
| 233 | /** Modifications below by kieranrcampbell@gmail.com |
|---|
| 234 | Institut Laue-Langevin, July 2012 |
|---|
| 235 | **/ |
|---|
| 236 | |
|---|
| 237 | /** |
|---|
| 238 | Implements eq 6.2.5 (small gamma) of Numerical Recipes in C, essentially |
|---|
| 239 | the incomplete gamma function multiplied by the gamma function. |
|---|
| 240 | Required for implementation of fast error function (erf) |
|---|
| 241 | **/ |
|---|
| 242 | |
|---|
| 243 | |
|---|
| 244 | #define ITMAX 100 |
|---|
| 245 | #define EPS 3.0e-7 |
|---|
| 246 | #define FPMIN 1.0e-30 |
|---|
| 247 | |
|---|
| 248 | void gser(float *gamser, float a, float x, float *gln) { |
|---|
| 249 | int n; |
|---|
| 250 | float sum,del,ap; |
|---|
| 251 | |
|---|
| 252 | *gln = gamln(a); |
|---|
| 253 | if(x <= 0.0) { |
|---|
| 254 | if (x < 0.0) printf("Error: x less than 0 in routine gser"); |
|---|
| 255 | *gamser = 0.0; |
|---|
| 256 | return; |
|---|
| 257 | } else { |
|---|
| 258 | ap = a; |
|---|
| 259 | del = sum = 1.0/a; |
|---|
| 260 | |
|---|
| 261 | for(n=1;n<=ITMAX;n++) { |
|---|
| 262 | ++ap; |
|---|
| 263 | del *= x/ap; |
|---|
| 264 | sum += del; |
|---|
| 265 | if(fabs(del) < fabs(sum)*EPS) { |
|---|
| 266 | *gamser = sum * exp(-x + a * log(x) - (*gln)); |
|---|
| 267 | return; |
|---|
| 268 | } |
|---|
| 269 | } |
|---|
| 270 | printf("a too large, ITMAX too small in routine gser"); |
|---|
| 271 | return; |
|---|
| 272 | |
|---|
| 273 | } |
|---|
| 274 | |
|---|
| 275 | |
|---|
| 276 | } |
|---|
| 277 | |
|---|
| 278 | /** |
|---|
| 279 | Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction |
|---|
| 280 | representation |
|---|
| 281 | **/ |
|---|
| 282 | |
|---|
| 283 | void gcf(float *gammcf, float a, float x, float *gln) { |
|---|
| 284 | int i; |
|---|
| 285 | float an,b,c,d,del,h; |
|---|
| 286 | |
|---|
| 287 | *gln = gamln(a); |
|---|
| 288 | b = x+1.0-a; |
|---|
| 289 | c = 1.0/FPMIN; |
|---|
| 290 | d = 1.0/b; |
|---|
| 291 | h=d; |
|---|
| 292 | for (i=1;i <= ITMAX; i++) { |
|---|
| 293 | an = -i*(i-a); |
|---|
| 294 | b += 2.0; |
|---|
| 295 | d = an*d + b; |
|---|
| 296 | if (fabs(d) < FPMIN) d = FPMIN; |
|---|
| 297 | c = b+an/c; |
|---|
| 298 | if (fabs(c) < FPMIN) c = FPMIN; |
|---|
| 299 | d = 1.0/d; |
|---|
| 300 | del = d*c; |
|---|
| 301 | h += del; |
|---|
| 302 | if (fabs(del-1.0) < EPS) break; |
|---|
| 303 | } |
|---|
| 304 | if (i > ITMAX) printf("a too large, ITMAX too small in gcf"); |
|---|
| 305 | *gammcf = exp(-x+a*log(x)-(*gln))*h; |
|---|
| 306 | return; |
|---|
| 307 | } |
|---|
| 308 | |
|---|
| 309 | |
|---|
| 310 | /** |
|---|
| 311 | Represents incomplete error function, P(a,x) |
|---|
| 312 | **/ |
|---|
| 313 | float gammp(float a, float x) { |
|---|
| 314 | float gamser,gammcf,gln; |
|---|
| 315 | if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp"); |
|---|
| 316 | if (x < (a+1.0)) { |
|---|
| 317 | gser(&gamser,a,x,&gln); |
|---|
| 318 | return gamser; |
|---|
| 319 | } else { |
|---|
| 320 | gcf(&gammcf,a,x,&gln); |
|---|
| 321 | return 1.0 - gammcf; |
|---|
| 322 | } |
|---|
| 323 | } |
|---|
| 324 | |
|---|
| 325 | /** |
|---|
| 326 | Implementation of the error function, erf(x) |
|---|
| 327 | **/ |
|---|
| 328 | |
|---|
| 329 | float erff(float x) { |
|---|
| 330 | return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x); |
|---|
| 331 | } |
|---|
| 332 | |
|---|
