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

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since aab23e3 was a1daf86, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

hack around broken isfinite/isnan in tinycc

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