Mercurial > hg > octave-nkf

comparison scripts/specfun/expint.m @ 16586:f423873d3275

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Style fixes for ellipj.cc, ellipke.m, and expint.m * ellipj.cc, ellipke.m, expint.m: Style fixes.
author Mike Miller <mtmiller@ieee.org>
date Sun, 28 Apr 2013 18:50:51 -0400
parents 1a3bfb14b5da
children a3fdd6041e64
comparison
equal deleted inserted replaced
16585:1a3bfb14b5da 16586:f423873d3275
1 ## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com> 1 ## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com>
2 ## 2 ##
3 ## This program is free software; you can redistribute it and/or modify it under 3 ## This file is part of Octave.
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
6 ## version.
7 ## 4 ##
8 ## This program is distributed in the hope that it will be useful, but WITHOUT 5 ## Octave is free software; you can redistribute it and/or modify it
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 6 ## under the terms of the GNU General Public License as published by
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 7 ## the Free Software Foundation; either version 3 of the License, or (at
11 ## details. 8 ## your option) any later version.
12 ## 9 ##
13 ## You should have received a copy of the GNU General Public License along with 10 ## Octave is distributed in the hope that it will be useful, but
14 ## this program; if not, see <http://www.gnu.org/licenses/>. 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
15 18
16 ## -*- texinfo -*- 19 ## -*- texinfo -*-
17 ## @deftypefn {Function File} {@var{y} =} expint (@var{x}) 20 ## @deftypefn {Function File} {} expint (@var{x})
18 ## Compute the exponential integral, 21 ## Compute the exponential integral,
19 ## @verbatim 22 ## @verbatim
20 ## infinity 23 ## infinity
21 ## / 24 ## /
22 ## expint(x) = | exp(t)/t dt 25 ## expint(x) = | exp(t)/t dt
23 ## / 26 ## /
24 ## x 27 ## x
25 ## @end verbatim 28 ## @end verbatim
26 ## @seealso{expint_E1, expint_Ei}
27 ## @end deftypefn 29 ## @end deftypefn
28 30
29 function y = expint(x) 31 function y = expint (x)
32
30 if (nargin != 1) 33 if (nargin != 1)
31 print_usage; 34 print_usage ();
32 endif 35 endif
33 y = expint_E1(x); 36
37 y = expint_E1 (x);
38
34 endfunction 39 endfunction
35 40
36 ## -*- texinfo -*- 41 ## -*- texinfo -*-
37 ## @deftypefn {Function File} {@var{y} =} expint_E1 (@var{x}) 42 ## @deftypefn {Function File} {@var{y} =} expint_E1 (@var{x})
38 ## Compute the exponential integral, 43 ## Compute the exponential integral,
41 ## / 46 ## /
42 ## expint(x) = | exp(t)/t dt 47 ## expint(x) = | exp(t)/t dt
43 ## / 48 ## /
44 ## x 49 ## x
45 ## @end verbatim 50 ## @end verbatim
46 ## @seealso{expint, expint_Ei}
47 ## @end deftypefn 51 ## @end deftypefn
48 52
49 function y = expint_E1(x) 53 function y = expint_E1 (x)
54
50 if (nargin != 1) 55 if (nargin != 1)
51 print_usage; 56 print_usage ();
52 endif 57 endif
58
53 y = x; 59 y = x;
54 y(imag(x) > 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) > 0 & imag(x) != 0)) -i.*pi; 60
55 y(imag(x) < 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) < 0 & imag(x) != 0)) +i.*pi; 61 idx = (imag (x) > 0 & imag (x) != 0);
56 y(real(x) >= 0 & imag(x)==0) = -expint_Ei(-y(real(x) >= 0 & imag(x)==0)); 62 y(idx) = -expint_Ei (-y(idx)) - i.*pi;
57 y(real(x) < 0 & imag(x)==0) = -expint_Ei(-y(real(x) < 0 & imag(x)==0)) -i.*pi; 63
64 idx = (imag (x) < 0 & imag (x) != 0);
65 y(idx) = -expint_Ei (-y(idx)) + i.*pi;
66
67 idx = (real (x) >= 0 & imag (x) == 0);
68 y(idx) = -expint_Ei (-y(idx));
69
70 idx = (real (x) < 0 & imag (x) == 0);
71 y(idx) = -expint_Ei (-y(idx)) - i.*pi;
72
58 endfunction 73 endfunction
59 74
60 ## -*- texinfo -*- 75 ## -*- texinfo -*-
61 ## @deftypefn {Function File} {@var{y} =} expint_Ei (@var{x}) 76 ## @deftypefn {Function File} {@var{y} =} expint_Ei (@var{x})
62 ## Compute the exponential integral, 77 ## Compute the exponential integral,
65 ## / 80 ## /
66 ## expint_Ei(x) = - | exp(t)/t dt 81 ## expint_Ei(x) = - | exp(t)/t dt
67 ## / 82 ## /
68 ## -x 83 ## -x
69 ## @end verbatim 84 ## @end verbatim
70 ## @seealso{expint, expint_E1}
71 ## @end deftypefn 85 ## @end deftypefn
72 86
73 function y = expint_Ei(x) 87 function y = expint_Ei (x)
88
74 if (nargin != 1) 89 if (nargin != 1)
75 print_usage; 90 print_usage ();
76 endif 91 endif
77 y = zeros(size(x)); 92
78 F = @(x) exp(-x)./x; 93 y = zeros (size (x));
79 s = prod(size(x)); 94 F = @(x) exp (-x)./x;
95 s = prod (size (x));
96
80 for t = 1:s; 97 for t = 1:s;
81 if(x(t)<0 && imag(x(t)) == 0) 98 if (x(t) < 0 && imag (x(t)) == 0)
82 y(t) = -quad(F,-x(t),Inf); 99 y(t) = -quad (F, -x(t), Inf);
83 else 100 else
84 if(abs(x(t)) > 2 && imag(x(t)) == 0) 101 if (abs (x(t)) > 2 && imag (x(t)) == 0)
85 y(t) = expint_Ei(2) - quad(F,-x(t),-2); 102 y(t) = expint_Ei (2) - quad (F, -x(t), -2);
86 else 103 else
87 if(abs(x(t)) >= 10) 104 if (abs (x(t)) >= 10)
88 if(imag(x(t)) <= 0) 105 if (imag (x(t)) <= 0)
89 a1 = 4.03640; 106 a1 = 4.03640;
90 a2 = 1.15198; 107 a2 = 1.15198;
91 b1 = 5.03637; 108 b1 = 5.03637;
92 b2 = 4.19160; 109 b2 = 4.19160;
93 y(t) = -(x(t).^2 - a1.*x(t) + a2)./((x(t).^2-b1.*x(t)+b2).*(-x(t)).*exp(-x(t)))-i.*pi; 110 y(t) = -(x(t).^2 - a1.*x(t) + a2) ...
111 ./ ((x(t).^2 - b1.*x(t) + b2) .* (-x(t)) .* exp (-x(t))) ...
112 - i.*pi;
94 else 113 else
95 y(t) = conj(expint_Ei(conj(x(t)))); 114 y(t) = conj (expint_Ei (conj (x(t))));
96 endif; 115 endif;
97 ## Serie Expansion 116 ## Serie Expansion
98 else 117 else
99 for k = 1:100; 118 for k = 1:100;
100 y(t) = y(t) + x(t).^k./(k.*factorial(k)); 119 y(t) = y(t) + x(t).^k ./ (k.*factorial (k));
101 endfor 120 endfor
102 y(t) = 0.577215664901532860606512090082402431 + log(x(t)) + y(t); 121 y(t) = 0.577215664901532860606512090082402431 + log (x(t)) + y(t);
103 endif 122 endif
104 endif 123 endif
105 endif 124 endif
106 endfor 125 endfor
107 endfunction 126 endfunction