source: sasview/src/sas/sascalc/calculator/c_extensions/librefl.c @ f54e82cf

ESS_GUIESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since f54e82cf was f54e82cf, checked in by Torin Cooper-Bennun <torin.cooper-bennun@…>, 6 years ago

Cherry-pick changes from master for tinycc compatibility (note: tinycc compilation not yet supported fully):

144e032af2 tinycc doesn't return structures, so must pass return structure as pointer
a1daf86c0d hack around broken isfinite/isnan in tinycc
7e82256ecb declare all variables at the start of the block so C89 compilers don't complain

  • Property mode set to 100644
File size: 8.6 KB
Line 
1// The original code, of which work was not DANSE funded,
2// was provided by J. Cho.
3// And modified to fit sansmodels/sansview: JC
4
5#include <math.h>
6#include "librefl.h"
7#include <stdio.h>
8#include <stdlib.h>
9#if defined(_MSC_VER)
10#define NEED_ERF
11#endif
12
13
14
15#if defined(NEED_ERF)
16/* erf.c  - public domain implementation of error function erf(3m)
17
18reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
19            (New Algorithm handbook in C language) (Gijyutsu hyouron
20            sha, Tokyo, 1991) p.227 [in Japanese]                 */
21
22
23#ifdef _WIN32
24# include <float.h>
25# if !defined __MINGW32__ || defined __NO_ISOCEXT
26#  ifndef isnan
27#   define isnan(x) _isnan(x)
28#  endif
29#  ifndef isinf
30#   define isinf(x) (!_finite(x) && !_isnan(x))
31#  endif
32#  ifndef finite
33#   define finite(x) _finite(x)
34#  endif
35# endif
36#endif
37
38static double q_gamma(double, double, double);
39
40/* Incomplete gamma function
41   1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  */
42static double p_gamma(double a, double x, double loggamma_a)
43{
44    int k;
45    double result, term, previous;
46
47    if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
48    if (x == 0)     return 0;
49    result = term = exp(a * log(x) - x - loggamma_a) / a;
50    for (k = 1; k < 1000; k++) {
51        term *= x / (a + k);
52        previous = result;  result += term;
53        if (result == previous) return result;
54    }
55    fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
56    return result;
57}
58
59/* Incomplete gamma function
60   1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  */
61static double q_gamma(double a, double x, double loggamma_a)
62{
63    int k;
64    double result, w, temp, previous;
65    double la = 1, lb = 1 + x - a;  /* Laguerre polynomial */
66
67    if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
68    w = exp(a * log(x) - x - loggamma_a);
69    result = w / lb;
70    for (k = 2; k < 1000; k++) {
71        temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
72        la = lb;  lb = temp;
73        w *= (k - 1 - a) / k;
74        temp = w / (la * lb);
75        previous = result;  result += temp;
76        if (result == previous) return result;
77    }
78    fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
79    return result;
80}
81
82#define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */
83
84double erf(double x)
85{
86    if (!finite(x)) {
87        if (isnan(x)) return x;      /* erf(NaN)   = NaN   */
88        return (x>0 ? 1.0 : -1.0);   /* erf(+-inf) = +-1.0 */
89    }
90    if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI_OVER_2);
91    else        return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
92}
93
94double erfc(double x)
95{
96    if (!finite(x)) {
97        if (isnan(x)) return x;      /* erfc(NaN)   = NaN      */
98        return (x>0 ? 0.0 : 2.0);    /* erfc(+-inf) = 0.0, 2.0 */
99    }
100    if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI_OVER_2);
101    else        return  1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
102}
103#endif // NEED_ERF
104
105void cassign(Cplx *x, double real, double imag)
106{
107        x->re = real;
108        x->im = imag;
109}
110
111
112void cplx_add(Cplx *z, Cplx x, Cplx y)
113{
114        z->re = x.re + y.re;
115        z->im = x.im + y.im;
116}
117
118void rcmult(Cplx *z, double x, Cplx y)
119{
120        z->re = x*y.re;
121        z->im = x*y.im;
122}
123
124void cplx_sub(Cplx *z, Cplx x, Cplx y)
125{
126        z->re = x.re - y.re;
127        z->im = x.im - y.im;
128}
129
130
131void cplx_mult(Cplx *z, Cplx x, Cplx y)
132{
133        z->re = x.re*y.re - x.im*y.im;
134        z->im = x.re*y.im + x.im*y.re;
135}
136
137void cplx_div(Cplx *z, Cplx x, Cplx y)
138{
139        z->re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
140        z->im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
141}
142
143void cplx_exp(Cplx *z, Cplx b)
144{
145        double br,bi;
146        br=b.re;
147        bi=b.im;
148        z->re = exp(br)*cos(bi);
149        z->im = exp(br)*sin(bi);
150}
151
152
153void cplx_sqrt(Cplx *c, Cplx z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
154{
155        double zr,zi,x,y,r,w;
156
157        zr=z.re;
158        zi=z.im;
159
160        if (zr==0.0 && zi==0.0)
161        {
162                c->re=0.0;
163                c->im=0.0;
164        } else {
165                x=fabs(zr);
166                y=fabs(zi);
167                if (x>y)
168                {
169                        r=y/x;
170                        w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
171                } else {
172                        r=x/y;
173                        w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
174                }
175                if (zr >=0.0)
176                {
177                        c->re=w;
178                        c->im=zi/(2.0*w);
179                } else {
180                        c->im=(zi >= 0) ? w : -w;
181                        c->re=zi/(2.0*c->im);
182                }
183        }
184}
185
186void cplx_cos(Cplx *z, Cplx b)
187{
188        // cos(b) = (e^bi + e^-bi)/2
189        //        = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2
190        //        = (e^b.im cos(-b.re) + e^b.im sin(-b.re) i)/2 + (e^-b.im cos(b.re) + e^-b.im sin(b.re) i)/2
191        //        = e^b.im cos(b.re)/2 - e^b.im sin(b.re)/2 i + 1/e^b.im cos(b.re)/2 + 1/e^b.im sin(b.re)/2 i
192        //        = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i
193        //        = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i
194        double exp_b_im = exp(b.im);
195        z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re);
196        z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re);
197}
198
199// normalized and modified erf
200//   |
201// 1 +                __  - - - -
202//   |             _
203//       |            _
204//   |        __
205// 0 + - - -
206//   |-------------+------------+--
207//   0           center       n_sub    --->
208//                                     ind
209//
210// n_sub = total no. of bins(or sublayers)
211// ind = x position: 0 to max
212// nu = max x to integration
213double err_mod_func(double n_sub, double ind, double nu)
214{
215  double center, func;
216  if (nu == 0.0)
217                nu = 1e-14;
218        if (n_sub == 0.0)
219                n_sub = 1.0;
220
221
222        //ind = (n_sub-1.0)/2.0-1.0 +ind;
223        center = n_sub/2.0;
224        // transform it so that min(ind) = 0
225        ind -= center;
226        // normalize by max limit
227        ind /= center;
228        // divide by sqrt(2) to get Gaussian func
229        nu /= sqrt(2.0);
230        ind *= nu;
231        // re-scale and normalize it so that max(erf)=1, min(erf)=0
232        func = erf(ind)/erf(nu)/2.0;
233        // shift it by +0.5 in y-direction so that min(erf) = 0
234        func += 0.5;
235
236        return func;
237}
238double linearfunc(double n_sub, double ind, double nu)
239{
240  double bin_size, func;
241        if (n_sub == 0.0)
242                n_sub = 1.0;
243
244        bin_size = 1.0/n_sub;  //size of each sub-layer
245        // rescale
246        ind *= bin_size;
247        func = ind;
248
249        return func;
250}
251// use the right hand side from the center of power func
252double power_r(double n_sub, double ind, double nu)
253{
254  double bin_size,func;
255        if (nu == 0.0)
256                nu = 1e-14;
257        if (n_sub == 0.0)
258                n_sub = 1.0;
259
260        bin_size = 1.0/n_sub;  //size of each sub-layer
261        // rescale
262        ind *= bin_size;
263        func = pow(ind, nu);
264
265        return func;
266}
267// use the left hand side from the center of power func
268double power_l(double n_sub, double ind, double nu)
269{
270  double bin_size, func;
271        if (nu == 0.0)
272                nu = 1e-14;
273        if (n_sub == 0.0)
274                n_sub = 1.0;
275
276        bin_size = 1.0/n_sub;  //size of each sub-layer
277        // rescale
278        ind *= bin_size;
279        func = 1.0-pow((1.0-ind),nu);
280
281        return func;
282}
283// use 1-exp func from x=0 to x=1
284double exp_r(double n_sub, double ind, double nu)
285{
286  double bin_size, func;
287        if (nu == 0.0)
288                nu = 1e-14;
289        if (n_sub == 0.0)
290                n_sub = 1.0;
291
292        bin_size = 1.0/n_sub;  //size of each sub-layer
293        // rescale
294        ind *= bin_size;
295        // modify func so that func(0) =0 and func(max)=1
296        func = 1.0-exp(-nu*ind);
297        // normalize by its max
298        func /= (1.0-exp(-nu));
299
300        return func;
301}
302
303// use the left hand side mirror image of exp func
304double exp_l(double n_sub, double ind, double nu)
305{
306  double bin_size, func;
307        if (nu == 0.0)
308                nu = 1e-14;
309        if (n_sub == 0.0)
310                n_sub = 1.0;
311
312        bin_size = 1.0/n_sub;  //size of each sub-layer
313        // rescale
314        ind *= bin_size;
315        // modify func
316        func = exp(-nu*(1.0-ind))-exp(-nu);
317        // normalize by its max
318        func /= (1.0-exp(-nu));
319
320        return func;
321}
322
323// To select function called
324// At nu = 0 (singular point), call line function
325double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r)
326{
327        double sld_i, func;
328        // this condition protects an error from the singular point
329        if (nu == 0.0){
330                nu = 1e-13;
331        }
332        // select func
333        switch(fun_type){
334                case 1 :
335                        func = power_r(n_sub, i, nu);
336                        break;
337                case 2 :
338                        func = power_l(n_sub, i, nu);
339                        break;
340                case 3 :
341                        func = exp_r(n_sub, i, nu);
342                        break;
343                case 4 :
344                        func = exp_l(n_sub, i, nu);
345                        break;
346                case 5 :
347                        func = linearfunc(n_sub, i, nu);
348                        break;
349                default:
350                        func = err_mod_func(n_sub, i, nu);
351                        break;
352        }
353        // compute sld
354        if (sld_r>sld_l){
355                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
356        }
357        else if (sld_r<sld_l){
358                func = 1.0-func;
359                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
360        }
361        else{
362                sld_i = sld_r;
363        }
364        return sld_i;
365}
366
367
368// used by refl.c
369double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r)
370{
371        double sld_i, func;
372        switch(fun_type){
373                case 0 :
374                        func = err_mod_func(n_sub, i, 2.5);
375                        break;
376                default:
377                        func = linearfunc(n_sub, i, 1.0);
378                        break;
379        }
380        if (sld_r>sld_l){
381                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
382        }
383        else if (sld_r<sld_l){
384                func = 1.0-func;
385                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
386        }
387        else{
388                sld_i = sld_r;
389        }
390        return sld_i;
391}
Note: See TracBrowser for help on using the repository browser.