source: sasview/src/sas/models/c_extension/c_models/libfunc.c @ be0c318

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since be0c318 was 79492222, checked in by krzywon, 10 years ago

Changed the file and folder names to remove all SANS references.
  • Property mode set to 100644
File size: 7.0 KB
Line 
1// by jcho
2
3#include <math.h>
4
5#include "libmultifunc/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
90double 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.
114polar_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
248void 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
283void 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**/
313float 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
329float erff(float x) {
330  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
331}
332
Note: See TracBrowser for help on using the repository browser.