source: sasview/src/sas/sascalc/calculator/c_extensions/libfunc.c @ 45eddeb7

Last change on this file since 45eddeb7 was 144e032a, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

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