1 /*
2 * $Id: dawson.i,v 1.2 2007-11-10 20:03:49 dhmunro Exp $
3 * Dawson's integral and error functions after SLATEC
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
dawson(x)11 func dawson(x)
12 /* DOCUMENT dawson(x)
13 * return Dawson's integral, exp(-x^2)*integral[0 to x](exp(t^2)*dt)
14 * maximum is dawson(0.9241388730) = 0.5410442246
15 * inflection point is dawson(1.5019752682) = 0.4276866160
16 * SEE ALSO: erf, erfc
17 */
18 {
19 { local mask1,mask2,mask3,xx,yy; }
20 if (structof(x)==complex) error, "dawson function not valid for complex";
21 y = abs(x);
22 if (any_in(,y,1., mask1, yy, x,xx))
23 d1 = xx * (0.75 + _cheby_eval(_dawson_1, 2.*yy*yy-1.));
24 if (any_in(1.,y,4., mask2, yy, x,xx))
25 d2 = xx * (0.25 + _cheby_eval(_dawson_2, 0.125*yy*yy-1.));
26 if (any_in(4.,y,, mask3, yy, x,xx))
27 d3 = (0.50 + _cheby_eval(_dawson_3, 32./(yy*yy)-1.)) / xx;
28 return merge_n(d1,mask1, d2,mask2, d3,mask3);
29 }
30
31 _dawson_1 = [-.006351734375145949, -.229407147967738690, .022130500939084764,
32 -.001549265453892985, .000084973277156849, -.000003828266270972,
33 .000000146285480625, -.000000004851982381, .000000000142146357,
34 -.000000000003728836, .000000000000088549, -.000000000000001920,
35 .000000000000000038];
36 _dawson_2 = [-.056886544105215527, -.318113469961681310, .208738454136422370,
37 -.124754099137791310, .067869305186676777, -.033659144895270940,
38 .015260781271987972, -.006348370962596214, .002432674092074852,
39 -.000862195414910650, .000283765733363216, -.000087057549874170,
40 .000024986849985481, -.000006731928676416, .000001707857878557,
41 -.000000409175512264, .000000092828292216, -.000000019991403610,
42 .000000004096349064, -.000000000800324095, .000000000149385031,
43 -.000000000026687999, .000000000004571221, -.000000000000751873,
44 .000000000000118931, -.000000000000018116, .000000000000002661,
45 -.000000000000000377, .000000000000000051];
46 _dawson_3 = [ .01690485637765704, .00868325227840695, .00024248640424177,
47 .00001261182399572, .00000106645331463, .00000013581597947,
48 .00000002171042356, .00000000286701050, -.00000000019013363,
49 -.00000000030977804, -.00000000010294148, -.00000000000626035,
50 .00000000000856313, .00000000000303304, -.00000000000025236,
51 -.00000000000042106, -.00000000000004431, .00000000000004911,
52 .00000000000001235, -.00000000000000578, -.00000000000000228,
53 .00000000000000076, .00000000000000038, -.00000000000000011,
54 -.00000000000000006, .00000000000000002];
55
erf(x)56 func erf(x)
57 /* DOCUMENT erf(x)
58 * return erf(x), 2./sqrt(pi) * integral[0 to x](exp(-t^2)*dt)
59 * SEE ALSO: erfc, dawson
60 */
61 {
62 { local mask1,mask2,mask3,mask4,xx,yy,zz; }
63 if (structof(x)==complex) error, "erf function not valid for complex";
64 y = abs(x);
65 if (any_in(,y,1., mask1, yy, x,xx))
66 e1 = xx * (1.0 + _cheby_eval(_erf_1, 2.*yy*yy-1.));
67 if (any_in(1.,y,, mask2, yy, x,xx)) {
68 z = yy*yy;
69 yy = exp(-z)/yy;
70 if (any_in(,z,4., mask3, zz, yy,y))
71 e3 = y * (0.5 + _cheby_eval(_erfc_1, (8./zz-5.)/3.));
72 if (any_in(4.,z,, mask4, zz, yy,y))
73 e4 = y * (0.5 + _cheby_eval(_erfc_2, 8./zz-1.));
74 e2 = sign(xx) * (1.0 - merge_n(e3,mask3, e4,mask4));
75 }
76 return merge_n(e1,mask1, e2,mask2);
77 }
78
erfc(x)79 func erfc(x)
80 /* DOCUMENT erfc(x)
81 * return erfc(x), 2./sqrt(pi) * integral[x to infinity](exp(-t^2)*dt)
82 * SEE ALSO: erf, dawson
83 */
84 {
85 { local mask1,mask2,mask3,mask4,xx,yy,zz; }
86 if (structof(x)==complex) error, "erfc function not valid for complex";
87 y = abs(x);
88 if (any_in(,y,1., mask1, yy, x,xx))
89 e1 = 1.0 - xx * (1.0 + _cheby_eval(_erf_1, 2.*yy*yy-1.));
90 if (any_in(1.,y,, mask2, yy, x,xx)) {
91 z = yy*yy;
92 yy = exp(-z)/yy;
93 if (any_in(,z,4., mask3, zz, yy,y))
94 e3 = y * (0.5 + _cheby_eval(_erfc_1, (8./zz-5.)/3.));
95 if (any_in(4.,z,, mask4, zz, yy,y))
96 e4 = y * (0.5 + _cheby_eval(_erfc_2, 8./zz-1.));
97 e2 = 2.*(xx<0.) + sign(xx)*merge_n(e3,mask3, e4,mask4);
98 }
99 return merge_n(e1,mask1, e2,mask2);
100 }
101
102 _erf_1 = [-.049046121234691808, -.142261205103713640, .010035582187599796,
103 -.000576876469976748, .000027419931252196, -.000001104317550734,
104 .000000038488755420, -.000000001180858253, .000000000032334215,
105 -.000000000000799101, .000000000000017990, -.000000000000000371,
106 .000000000000000007];
107 _erfc_1 = [-.069601346602309501, -.041101339362620893, .003914495866689626,
108 -.000490639565054897, .000071574790013770, -.000011530716341312,
109 .000001994670590201, -.000000364266647159, .000000069443726100,
110 -.000000013712209021, .000000002788389661, -.000000000581416472,
111 .000000000123892049, -.000000000026906391, .000000000005942614,
112 -.000000000001332386, .000000000000302804, -.000000000000069666,
113 .000000000000016208, -.000000000000003809, .000000000000000904,
114 -.000000000000000216, .000000000000000052];
115 _erfc_2 = [ .071517931020292500, -.026532434337606719, .001711153977920853,
116 -.000163751663458512, .000019871293500549, -.000002843712412769,
117 .000000460616130901, -.000000082277530261, .000000015921418724,
118 -.000000003295071356, .000000000722343973, -.000000000166485584,
119 .000000000040103931, -.000000000010048164, .000000000002608272,
120 -.000000000000699105, .000000000000192946, -.000000000000054704,
121 .000000000000015901, -.000000000000004729, .000000000000001432,
122 -.000000000000000439, .000000000000000138, -.000000000000000048];
123
_cheby_eval(cs,x)124 func _cheby_eval(cs, x)
125 {
126 n = numberof(cs);
127 b0 = b1 = 0.;
128 x *= 2.;
129 for (i=0 ; i<n ; ++i) {
130 b2 = b1;
131 b1 = b0;
132 b0 = x*b1 - b2 + cs(n-i);
133 }
134 return 0.5*(b0-b2);
135 }
136
137 func any_in(left,x,right, &mask, &xx, a,&aa, b,&bb, c,&cc)
138 /* DOCUMENT any_in(left,x,right, mask, xx, a,aa, b,bb, c,cc
139 * return the number of elements of the array X which are in the
140 * interval LEFT < X <= RIGHT. Also return MASK, which has the
141 * shape of X and is 1 where X is in the interval and 0 otherwise,
142 * and XX = X(where(MASK)). Up to three optional arrays A, B, and C
143 * of the same shape as X may be supplied; the arrays AA, BB, and CC
144 * analogous to XX are returned. LEFT or RIGHT may be [] for the
145 * interval to extend to infinity on the corresponding side.
146 * LEFT and/or RIGHT may be arrays as long as they are conformable
147 * with X.
148 * SEE ALSO: merge_n, merge
149 */
150 {
151 if (is_void(left)) mask = array(1n, dimsof(x));
152 else mask = (x>left);
153 if (!is_void(right)) mask &= (x<=right);
154 list = where(mask);
155 n = numberof(list);
156 if (n) {
157 xx = x(list);
158 if (!is_void(a)) {
159 aa = a(list);
160 if (!is_void(b)) {
161 bb = b(list);
162 if (!is_void(c)) cc = c(list);
163 }
164 }
165 }
166 return n;
167 }
168
merge_n(val,mask,..)169 func merge_n(val,mask, ..)
170 /* DOCUMENT merge_n(val1,mask1, val2,mask2, ...)
171 * return array with shape of MASKi (which must all have same shape)
172 * and values VALi where MASKi is true. Unspecified values will be
173 * zero; the data type of the result is the data type of the first
174 * non-nil VALi. Each VALi must be a 1D array of length sum(MASKi!=0).
175 * SEE ALSO: any_in, merge
176 */
177 {
178 rslt = dimsof(mask);
179 while (is_void(val)) {
180 val = next_arg();
181 mask = next_arg();
182 if (is_void(mask)) return array(0., rslt);
183 }
184 rslt = array(structof(val), rslt);
185 for (;;) {
186 rslt(where(mask)) = val;
187 do {
188 val = next_arg();
189 mask = next_arg();
190 if (is_void(mask)) return rslt;
191 } while (is_void(val));
192 }
193 }
194