source: sasview/src/sas/sascalc/calculator/c_extensions/libfunc.c @ 1342f6a

magnetic_scatt
Last change on this file since 1342f6a was 1342f6a, checked in by dirk, 19 months ago

Make polarisation efficiency weights consistently the same as in SASmodels

  • Property mode set to 100644
File size: 7.5 KB
Line 
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
13int 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
38double 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
81double sinc(double x)
82{
83  if (x==0.0){
84    return 1.0;
85  }
86  return sin(x)/x;
87}
88
89//The code calculating magnetic sld below seems to be not used.
90//Keeping it in order to check if it is not breaking anything.
91// calculate magnetic sld and return total sld
92// bn : contrast (not just sld of the layer)
93// m0: max mag of M; mtheta: angle from x-z plane;
94// mphi: angle (anti-clock-wise)of x-z projection(M) from x axis
95// spinfraci: the fraction of UP among UP+Down (before sample)
96// spinfracf: the fraction of UP among UP+Down (after sample and before detector)
97// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis
98// Note: all angles are in degrees.
99// void cal_msld(polar_sld *p_sld, int isangle, double qx, double qy, double bn,
100//                              double m01, double mtheta1, double mphi1,
101//                              double spinfraci, double spinfracf, double spintheta)
102// {
103//      //locals
104//      double q_x = qx;
105//      double q_y = qy;
106//      double sld = bn;
107//      int is_angle = isangle;
108//      double pi = 4.0*atan(1.0);
109//      double s_theta = spintheta * pi/180.0;
110//      double m_max = m01;
111//      double m_phi = mphi1;
112//      double m_theta = mtheta1;
113//      double in_spin = spinfraci;
114//      double out_spin = spinfracf;
115//
116//      double m_perp = 0.0;
117//      double m_perp_z = 0.0;
118//      double m_perp_y = 0.0;
119//      double m_perp_x = 0.0;
120//      double m_sigma_x = 0.0;
121//      double m_sigma_z = 0.0;
122//      double m_sigma_y = 0.0;
123//      //double b_m = 0.0;
124//      double q_angle = 0.0;
125//      double mx = 0.0;
126//      double my = 0.0;
127//      double mz = 0.0;
128//      double uu = sld;
129//      double dd = sld;
130//      double re_ud = 0.0;
131//      double im_ud = 0.0;
132//      double re_du = 0.0;
133//      double im_du = 0.0;
134//  const double norm=outspin;
135
136
137// The norm is needed to make sure that the scattering cross sections are
138//correctly weighted, such that the sum of spin-resolved measurements adds up to
139// the unpolarised or half-polarised scattering cross section. No intensity weighting
140// needed on the incoming polariser side (assuming that a user), has normalised
141// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.
142
143//if (out_spin < 0.5){norm=1-out_spin;}
144//else{norm=out_spin;}
145
146//      //No mag means no further calculation
147//      if (isangle>0) {
148//              if (m_max < 1.0e-32){
149//                      uu = sqrt(in_spin * out_spin/ norm) * uu ;
150//                      dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd ;
151//              }
152//      }
153//      else if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
154//                      uu = sqrt(in_spin * out_spin/ norm) * uu;
155//                      dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd;
156//      } else {
157//
158//              //These are needed because of the precision of inputs
159//              if (in_spin < 0.0) in_spin = 0.0;
160//              if (in_spin > 1.0) in_spin = 1.0;
161//              if (out_spin < 0.0) out_spin = 0.0;
162//              if (out_spin > 1.0) out_spin = 1.0;
163//
164//              if (q_x == 0.0) q_angle = pi / 2.0;
165//              else q_angle = atan(q_y/q_x);
166//              if (q_y < 0.0 && q_x < 0.0) q_angle -= pi;
167//              else if (q_y > 0.0 && q_x < 0.0) q_angle += pi;
168//
169//              q_angle = pi/2.0 - q_angle;
170//              if (q_angle > pi) q_angle -= 2.0 * pi;
171//              else if (q_angle < -pi) q_angle += 2.0 * pi;
172//
173//              if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){
174//                      m_perp = 0.0;
175//                      }
176//              else {
177//                      m_perp = m_max;
178//                      }
179//              if (is_angle > 0){
180//                      m_phi *= pi/180.0;
181//                      m_theta *= pi/180.0;
182//                      mx = m_perp * cos(m_theta) * cos(m_phi);
183//                      my = m_perp * sin(m_theta);
184//                      mz = -(m_perp * cos(m_theta) * sin(m_phi));
185//              }
186//              else{
187//                      mx = m_perp;
188//                      my = m_phi;
189//                      mz = m_theta;
190//              }
191//              //ToDo: simplify these steps
192//              // m_perp1 -m_perp2
193//              m_perp_x = (mx) *  cos(q_angle);
194//              m_perp_x -= (my) * sin(q_angle);
195//              m_perp_y = m_perp_x;
196//              m_perp_x *= cos(-q_angle);
197//              m_perp_y *= sin(-q_angle);
198//              m_perp_z = mz;
199//
200//              m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta));
201//              m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta));
202//              m_sigma_z = (m_perp_z);
203//
204//              //Find b
205//              uu += m_sigma_x;
206//              dd -= m_sigma_x;
207//              re_ud = m_sigma_y;
208//              re_du = m_sigma_y;
209//              im_ud = m_sigma_z;
210//              im_du = -m_sigma_z;
211//
212//              uu = sqrt((in_spin) * (out_spin)/ norm) * uu;
213//              dd = sqrt((1.0 - in_spin) * (1.0 - out_spin)/ norm) * dd;
214//
215//              re_ud = sqrt(in_spin * (1.0 - out_spin)/ norm) * re_ud;
216//              im_ud = sqrt(in_spin * (1.0 - out_spin)/ norm) * im_ud;
217//              re_du = sqrt((1.0 - in_spin) * out_spin/ norm) * re_du;
218//              im_du = sqrt((1.0 - in_spin) * out_spin/ norm) * im_du;
219//      }
220//      p_sld->uu = uu;
221//      p_sld->dd = dd;
222//      p_sld->re_ud = re_ud;
223//      p_sld->im_ud = im_ud;
224//      p_sld->re_du = re_du;
225//      p_sld->im_du = im_du;
226// }
227
228
229/** Modifications below by kieranrcampbell@gmail.com
230    Institut Laue-Langevin, July 2012
231**/
232
233/**
234Wojtek's comment Mar 22 2016: The remaing code can mostly likely be deleated
235Keeping it in order to check if it is not breaking anything
236**/
237
238/*
239#define ITMAX 100
240#define EPS 3.0e-7
241#define FPMIN 1.0e-30
242
243void gser(float *gamser, float a, float x, float *gln) {
244  int n;
245  float sum,del,ap;
246
247  *gln = lgamma(a);
248  if(x <= 0.0) {
249    if (x < 0.0) printf("Error: x less than 0 in routine gser");
250    *gamser = 0.0;
251    return;
252  } else {
253    ap = a;
254    del = sum = 1.0/a;
255
256    for(n=1;n<=ITMAX;n++) {
257      ++ap;
258      del *= x/ap;
259      sum += del;
260      if(fabs(del) < fabs(sum)*EPS) {
261        *gamser = sum * exp(-x + a * log(x) - (*gln));
262        return;
263      }
264    }
265    printf("a too large, ITMAX too small in routine gser");
266    return;
267
268  }
269
270}
271*/
272/**
273   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction
274   representation
275**/
276/*
277void gcf(float *gammcf, float a, float x, float *gln) {
278  int i;
279  float an,b,c,d,del,h;
280
281  *gln = lgamma(a);
282  b = x+1.0-a;
283  c = 1.0/FPMIN;
284  d = 1.0/b;
285  h=d;
286  for (i=1;i <= ITMAX; i++) {
287    an = -i*(i-a);
288    b += 2.0;
289    d = an*d + b;
290    if (fabs(d) < FPMIN) d = FPMIN;
291    c = b+an/c;
292    if (fabs(c) < FPMIN) c = FPMIN;
293    d = 1.0/d;
294    del = d*c;
295    h += del;
296    if (fabs(del-1.0) < EPS) break;
297  }
298  if (i > ITMAX) printf("a too large, ITMAX too small in gcf");
299  *gammcf = exp(-x+a*log(x)-(*gln))*h;
300  return;
301}
302*/
303/**
304   Represents incomplete error function, P(a,x)
305**/
306/*
307float gammp(float a, float x) {
308  float gamser,gammcf,gln;
309  if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp");
310  if (x < (a+1.0)) {
311    gser(&gamser,a,x,&gln);
312    return gamser;
313  } else {
314    gcf(&gammcf,a,x,&gln);
315    return 1.0 - gammcf;
316  }
317}
318*/
319/**
320    Implementation of the error function, erf(x)
321**/
322/*
323float erff(float x) {
324  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
325}
326*/
Note: See TracBrowser for help on using the repository browser.