source: sasview/sansmodels/src/c_models/libfunc.c @ c2e5898

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 c2e5898 was f8644dd, checked in by Kieran Campbell <kieranrcampbell@…>, 12 years ago

Fixed calculation functions to pass unit tests
  • Property mode set to 100644
File size: 1.6 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{
16
17        int k, j;
18
19        if (i<2){
20
21                return 1;
22
23        }
24
25        k=1;
26
27        for(j=1;j<i;j++)
28
29        {
30
31                k=k*(j+1);
32
33        }
34
35        return k;
36
37}
38
39
40
41// Used in pearl nec model
42
43// Sine integral function: approximated within 1%!!!
44
45// integral of sin(x)/x up to namx term nmax=6 looks the best.
46
47double Si(double x)
48
49{
50
51        int i;
52
53        int nmax=6;
54
55        double out;
56
57        long double power;
58
59        double pi = 4.0*atan(1.0);
60
61        if (x >= pi*6.2/4.0){
62
63
64
65                double out_sin = 0.0;
66
67                double out_cos = 0.0;
68
69                out = pi/2.0;
70
71                for (i=0; i<nmax-2; i+=1){
72
73                        out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
74
75                        out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
76
77                }
78
79                out -= cos(x) * out_cos;
80
81                out -= sin(x) * out_sin;
82
83                return out;
84
85        }
86
87        out = 0.0;
88
89        for (i=0; i<nmax; i+=1)
90
91        {
92
93                if (i==0){
94
95                        out += x;
96
97                        continue;
98
99                }
100
101                power = pow(x,(2 * i + 1));
102
103                out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));
104
105                //printf ("Si=%g %g %d\n", x, out, i);
106
107        }
108
109        return out;
110
111}
112
113
114
115double sinc(double x)
116{
117  if (x==0.0){
118    return 1.0;
119  }
120  return sin(x)/x;
121}
122
123
124double gamln(double xx) {
125
126    double x,y,tmp,ser;
127    static double cof[6]={76.18009172947146,-86.50532032941677,
128    24.01409824083091,-1.231739572450155,
129    0.1208650973866179e-2,-0.5395239384953e-5};
130    int j;
131
132    y=x=xx;
133    tmp=x+5.5;
134    tmp -= (x+0.5)*log(tmp);
135    ser=1.000000000190015;
136    for (j=0;j<=5;j++) ser += cof[j]/++y;
137    return -tmp+log(2.5066282746310005*ser/x);
138}
139
Note: See TracBrowser for help on using the repository browser.