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

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 4c29e4d was 4c29e4d, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

remove compiler warnings from mac build

  • Property mode set to 100644
File size: 8.5 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
105complex cassign(real, imag)
106        double real, imag;
107{
108        complex x;
109        x.re = real;
110        x.im = imag;
111        return x;
112}
113
114
115complex cplx_add(x,y)
116        complex x,y;
117{
118        complex z;
119        z.re = x.re + y.re;
120        z.im = x.im + y.im;
121        return z;
122}
123
124complex rcmult(x,y)
125        double x;
126    complex y;
127{
128        complex z;
129        z.re = x*y.re;
130        z.im = x*y.im;
131        return z;
132}
133
134complex cplx_sub(x,y)
135        complex x,y;
136{
137        complex z;
138        z.re = x.re - y.re;
139        z.im = x.im - y.im;
140        return z;
141}
142
143
144complex cplx_mult(x,y)
145        complex x,y;
146{
147        complex z;
148        z.re = x.re*y.re - x.im*y.im;
149        z.im = x.re*y.im + x.im*y.re;
150        return z;
151}
152
153complex cplx_div(x,y)
154        complex x,y;
155{
156        complex z;
157        z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
158        z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
159        return z;
160}
161
162complex cplx_exp(b)
163        complex b;
164{
165        complex z;
166        double br,bi;
167        br=b.re;
168        bi=b.im;
169        z.re = exp(br)*cos(bi);
170        z.im = exp(br)*sin(bi);
171        return z;
172}
173
174
175complex cplx_sqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
176        complex z;
177{
178        complex c;
179        double zr,zi,x,y,r,w;
180
181        zr=z.re;
182        zi=z.im;
183
184        if (zr==0.0 && zi==0.0)
185        {
186    c.re=0.0;
187        c.im=0.0;
188    return c;
189        }
190        else
191        {
192                x=fabs(zr);
193                y=fabs(zi);
194                if (x>y)
195                {
196                        r=y/x;
197                        w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
198                }
199                else
200                {
201                        r=x/y;
202                        w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
203                }
204                if (zr >=0.0)
205                {
206                        c.re=w;
207                        c.im=zi/(2.0*w);
208                }
209                else
210                {
211                        c.im=(zi >= 0) ? w : -w;
212                        c.re=zi/(2.0*c.im);
213                }
214                return c;
215        }
216}
217
218complex cplx_cos(b)
219        complex b;
220{
221        complex zero,two,z,i,bi,negbi;
222        zero = cassign(0.0,0.0);
223        two = cassign(2.0,0.0);
224        i = cassign(0.0,1.0);
225        bi = cplx_mult(b,i);
226        negbi = cplx_sub(zero,bi);
227        z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two);
228        return z;
229}
230
231// normalized and modified erf
232//   |
233// 1 +                __  - - - -
234//   |             _
235//       |            _
236//   |        __
237// 0 + - - -
238//   |-------------+------------+--
239//   0           center       n_sub    --->
240//                                     ind
241//
242// n_sub = total no. of bins(or sublayers)
243// ind = x position: 0 to max
244// nu = max x to integration
245double err_mod_func(double n_sub, double ind, double nu)
246{
247  double center, func;
248  if (nu == 0.0)
249                nu = 1e-14;
250        if (n_sub == 0.0)
251                n_sub = 1.0;
252
253
254        //ind = (n_sub-1.0)/2.0-1.0 +ind;
255        center = n_sub/2.0;
256        // transform it so that min(ind) = 0
257        ind -= center;
258        // normalize by max limit
259        ind /= center;
260        // divide by sqrt(2) to get Gaussian func
261        nu /= sqrt(2.0);
262        ind *= nu;
263        // re-scale and normalize it so that max(erf)=1, min(erf)=0
264        func = erf(ind)/erf(nu)/2.0;
265        // shift it by +0.5 in y-direction so that min(erf) = 0
266        func += 0.5;
267
268        return func;
269}
270double linearfunc(double n_sub, double ind, double nu)
271{
272  double bin_size, func;
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 = ind;
280
281        return func;
282}
283// use the right hand side from the center of power func
284double power_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        func = pow(ind, nu);
296
297        return func;
298}
299// use the left hand side from the center of power func
300double power_l(double n_sub, double ind, double nu)
301{
302  double bin_size, func;
303        if (nu == 0.0)
304                nu = 1e-14;
305        if (n_sub == 0.0)
306                n_sub = 1.0;
307
308        bin_size = 1.0/n_sub;  //size of each sub-layer
309        // rescale
310        ind *= bin_size;
311        func = 1.0-pow((1.0-ind),nu);
312
313        return func;
314}
315// use 1-exp func from x=0 to x=1
316double exp_r(double n_sub, double ind, double nu)
317{
318  double bin_size, func;
319        if (nu == 0.0)
320                nu = 1e-14;
321        if (n_sub == 0.0)
322                n_sub = 1.0;
323
324        bin_size = 1.0/n_sub;  //size of each sub-layer
325        // rescale
326        ind *= bin_size;
327        // modify func so that func(0) =0 and func(max)=1
328        func = 1.0-exp(-nu*ind);
329        // normalize by its max
330        func /= (1.0-exp(-nu));
331
332        return func;
333}
334
335// use the left hand side mirror image of exp func
336double exp_l(double n_sub, double ind, double nu)
337{
338  double bin_size, func;
339        if (nu == 0.0)
340                nu = 1e-14;
341        if (n_sub == 0.0)
342                n_sub = 1.0;
343
344        bin_size = 1.0/n_sub;  //size of each sub-layer
345        // rescale
346        ind *= bin_size;
347        // modify func
348        func = exp(-nu*(1.0-ind))-exp(-nu);
349        // normalize by its max
350        func /= (1.0-exp(-nu));
351
352        return func;
353}
354
355// To select function called
356// At nu = 0 (singular point), call line function
357double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r)
358{
359        double sld_i, func;
360        // this condition protects an error from the singular point
361        if (nu == 0.0){
362                nu = 1e-13;
363        }
364        // select func
365        switch(fun_type){
366                case 1 :
367                        func = power_r(n_sub, i, nu);
368                        break;
369                case 2 :
370                        func = power_l(n_sub, i, nu);
371                        break;
372                case 3 :
373                        func = exp_r(n_sub, i, nu);
374                        break;
375                case 4 :
376                        func = exp_l(n_sub, i, nu);
377                        break;
378                case 5 :
379                        func = linearfunc(n_sub, i, nu);
380                        break;
381                default:
382                        func = err_mod_func(n_sub, i, nu);
383                        break;
384        }
385        // compute sld
386        if (sld_r>sld_l){
387                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
388        }
389        else if (sld_r<sld_l){
390                func = 1.0-func;
391                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
392        }
393        else{
394                sld_i = sld_r;
395        }
396        return sld_i;
397}
398
399
400// used by refl.c
401double interfunc(int fun_type, double n_sub, double i, double sld_l, double sld_r)
402{
403        double sld_i, func;
404        switch(fun_type){
405                case 0 :
406                        func = err_mod_func(n_sub, i, 2.5);
407                        break;
408                default:
409                        func = linearfunc(n_sub, i, 1.0);
410                        break;
411        }
412        if (sld_r>sld_l){
413                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
414        }
415        else if (sld_r<sld_l){
416                func = 1.0-func;
417                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
418        }
419        else{
420                sld_i = sld_r;
421        }
422        return sld_i;
423}
Note: See TracBrowser for help on using the repository browser.