source: sasview/sansmodels/src/sans/models/c_extensions/libmultifunc/librefl.c @ 2b0f822

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 2b0f822 was 0164899a, checked in by Jae Cho <jhjcho@…>, 14 years ago

new models and some bug fixes

  • Property mode set to 100644
File size: 6.1 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 "libmultifunc/librefl.h"
7#include <stdio.h>
8#include <stdlib.h>
9
10
11complex cassign(real, imag)
12        double real, imag;
13{
14        complex x;
15        x.re = real;
16        x.im = imag;
17        return x;
18}
19
20
21complex cadd(x,y)
22        complex x,y;
23{
24        complex z;
25        z.re = x.re + y.re;
26        z.im = x.im + y.im;
27        return z;
28}
29
30complex rcmult(x,y)
31        double x;
32    complex y;
33{
34        complex z;
35        z.re = x*y.re;
36        z.im = x*y.im;
37        return z;
38}
39
40complex csub(x,y)
41        complex x,y;
42{
43        complex z;
44        z.re = x.re - y.re;
45        z.im = x.im - y.im;
46        return z;
47}
48
49
50complex cmult(x,y)
51        complex x,y;
52{
53        complex z;
54        z.re = x.re*y.re - x.im*y.im;
55        z.im = x.re*y.im + x.im*y.re;
56        return z;
57}
58
59complex cdiv(x,y)
60        complex x,y;
61{
62        complex z;
63        z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
64        z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
65        return z;
66}
67
68complex cexp(b)
69        complex b;
70{
71        complex z;
72        double br,bi;
73        br=b.re;
74        bi=b.im;
75        z.re = exp(br)*cos(bi);
76        z.im = exp(br)*sin(bi);
77        return z;
78}
79
80
81complex csqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
82        complex z;
83{
84        complex c;
85        double zr,zi,x,y,r,w;
86
87        zr=z.re;
88        zi=z.im;
89
90        if (zr==0.0 && zi==0.0)
91        {
92    c.re=0.0;
93        c.im=0.0;
94    return c;
95        }
96        else
97        {
98                x=fabs(zr);
99                y=fabs(zi);
100                if (x>y)
101                {
102                        r=y/x;
103                        w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
104                }
105                else
106                {
107                        r=x/y;
108                        w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
109                }
110                if (zr >=0.0)
111                {
112                        c.re=w;
113                        c.im=zi/(2.0*w);
114                }
115                else
116                {
117                        c.im=(zi >= 0) ? w : -w;
118                        c.re=zi/(2.0*c.im);
119                }
120                return c;
121        }
122}
123
124complex ccos(b)
125        complex b;
126{
127        complex neg,negb,zero,two,z,i,bi,negbi;
128        zero = cassign(0.0,0.0);
129        two = cassign(2.0,0.0);
130        i = cassign(0.0,1.0);
131        bi = cmult(b,i);
132        negbi = csub(zero,bi);
133        z = cdiv(cadd(cexp(bi),cexp(negbi)),two);
134        return z;
135}
136
137// normalized and modified erf
138//   |
139// 1 +                __  - - - -
140//   |             _
141//       |            _
142//   |        __
143// 0 + - - -
144//   |-------------+------------+--
145//   0           center       n_sub    --->
146//                                     ind
147//
148// n_sub = total no. of bins(or sublayers)
149// ind = x position: 0 to max
150// nu = max x to integration
151double err_mod_func(n_sub, ind, nu)
152        double n_sub,nu, ind;
153{
154        if (nu == 0.0)
155                nu = 1e-14;
156        if (n_sub == 0.0)
157                n_sub = 1.0;
158
159        double center, func;
160        //ind = (n_sub-1.0)/2.0-1.0 +ind;
161        center = n_sub/2.0;
162        // transform it so that min(ind) = 0
163        ind -= center;
164        // normalize by max limit
165        ind /= center;
166        // divide by sqrt(2) to get Gaussian func
167        nu /= sqrt(2.0);
168        ind *= nu;
169        // re-scale and normalize it so that max(erf)=1, min(erf)=0
170        func = erf(ind)/erf(nu)/2.0;
171        // shift it by +0.5 in y-direction so that min(erf) = 0
172        func += 0.5;
173
174        return func;
175}
176double linearfunc(n_sub, ind, nu)
177        double n_sub,nu, ind;
178{
179        if (n_sub == 0.0)
180                n_sub = 1.0;
181
182        double bin_size, func;
183        bin_size = 1.0/n_sub;  //size of each sub-layer
184        // rescale
185        ind *= bin_size;
186        func = ind;
187
188        return func;
189}
190// use the right hand side from the center of power func
191double power_r(n_sub, ind, nu)
192        double n_sub, nu, ind;
193{
194        if (nu == 0.0)
195                nu = 1e-14;
196        if (n_sub == 0.0)
197                n_sub = 1.0;
198
199        double bin_size,func;
200        bin_size = 1.0/n_sub;  //size of each sub-layer
201        // rescale
202        ind *= bin_size;
203        func = pow(ind, nu);
204
205        return func;
206}
207// use the left hand side from the center of power func
208double power_l(n_sub, ind, nu)
209        double n_sub, nu, ind;
210{
211        if (nu == 0.0)
212                nu = 1e-14;
213        if (n_sub == 0.0)
214                n_sub = 1.0;
215
216        double bin_size, func;
217        bin_size = 1.0/n_sub;  //size of each sub-layer
218        // rescale
219        ind *= bin_size;
220        func = 1.0-pow((1.0-ind),nu);
221
222        return func;
223}
224// use 1-exp func from x=0 to x=1
225double exp_r(n_sub, ind, nu)
226        double n_sub, nu, ind;
227{
228        if (nu == 0.0)
229                nu = 1e-14;
230        if (n_sub == 0.0)
231                n_sub = 1.0;
232
233        double bin_size, func;
234        bin_size = 1.0/n_sub;  //size of each sub-layer
235        // rescale
236        ind *= bin_size;
237        // modify func so that func(0) =0 and func(max)=1
238        func = 1.0-exp(-nu*ind);
239        // normalize by its max
240        func /= (1.0-exp(-nu));
241
242        return func;
243}
244
245// use the left hand side mirror image of exp func
246double exp_l(n_sub, ind, nu)
247        double n_sub, nu, ind;
248{
249        if (nu == 0.0)
250                nu = 1e-14;
251        if (n_sub == 0.0)
252                n_sub = 1.0;
253
254        double bin_size, func;
255        bin_size = 1.0/n_sub;  //size of each sub-layer
256        // rescale
257        ind *= bin_size;
258        // modify func
259        func = exp(-nu*(1.0-ind))-exp(-nu);
260        // normalize by its max
261        func /= (1.0-exp(-nu));
262
263        return func;
264}
265
266// To select function called
267// At nu = 0 (singular point), call line function
268double intersldfunc(fun_type, n_sub, i, nu, sld_l, sld_r)
269        double n_sub, nu, sld_l, sld_r,i;
270        int fun_type;
271{
272        double sld_i, func;
273        // this condition protects an error from the singular point
274        if (nu == 0.0){
275                nu = 1e-13;
276        }
277        // select func
278        switch(fun_type){
279                case 1 :
280                        func = power_r(n_sub, i, nu);
281                        break;
282                case 2 :
283                        func = power_l(n_sub, i, nu);
284                        break;
285                case 3 :
286                        func = exp_r(n_sub, i, nu);
287                        break;
288                case 4 :
289                        func = exp_l(n_sub, i, nu);
290                        break;
291                case 5 :
292                        func = linearfunc(n_sub, i, nu);
293                        break;
294                default:
295                        func = err_mod_func(n_sub, i, nu);
296                        break;
297        }
298        // compute sld
299        if (sld_r>sld_l){
300                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
301        }
302        else if (sld_r<sld_l){
303                func = 1.0-func;
304                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
305        }
306        else{
307                sld_i = sld_r;
308        }
309        return sld_i;
310}
311
312
313// used by refl.c
314double interfunc(fun_type, n_sub, i, sld_l, sld_r)
315        double n_sub, sld_l, sld_r, i;
316        int fun_type;
317{
318        double sld_i, func;
319        switch(fun_type){
320                case 0 :
321                        func = err_mod_func(n_sub, i, 2.5);
322                        break;
323                default:
324                        func = linearfunc(n_sub, i, 1.0);
325                        break;
326        }
327        if (sld_r>sld_l){
328                sld_i = (sld_r-sld_l)*func+sld_l; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
329        }
330        else if (sld_r<sld_l){
331                func = 1.0-func;
332                sld_i = (sld_l-sld_r)*func+sld_r; //sld_cal(sld[i],sld[i+1],n_sub,dz,thick);
333        }
334        else{
335                sld_i = sld_r;
336        }
337        return sld_i;
338}
Note: See TracBrowser for help on using the repository browser.