## Copyright (C) 2006   Sylvain Pelissier   <sylvain.pelissier@gmail.com>
##               2013   Bastiaan Zapf       <bzapf@online.de>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

##SI compute the sine integral function define by:
##
##                    x
##                   /
##           Si(x) = | sin(t)/t dt
##                   /
##                   0
##

## In 2013, a part was added that uses Padé Expansions
## to efficiently compute real values.
##
## http://en.wikipedia.org/wiki/Trigonometric_integral

function y = Si(x)

	if (nargin != 1)
		usage ("Si(x)");
	endif

	if (all(real(x)))
		sx=x(find(abs(x)<4));
		lx=x(find(abs(x)>=4));

		y=zeros(size(x));
		
		if (length(sx))
			
		# Polynomials for values with magnitude < 4
		# Expansion of Si(x)/x
			
			sx2=sx.^2;
			
			NC=[ -6.05338212010422477e-16
			      7.08240282274875911e-13
			     -3.53201978997168357e-10
			      9.43280809438713025e-8
			     -1.41018536821330254e-5
			      1.15457225751016682e-3
			     -4.54393409816329991e-2
			      1. ]';
	
			DC=[ 3.21107051193712168e-16
			     4.5049097575386581e-13
			     3.28067571055789734e-10
			     1.55654986308745614e-7
			     4.99175116169755106e-5
			     1.01162145739225565e-2
			     1. ]';

			yN=polyval(NC,sx2);
			yD=polyval(DC,sx2);

			sy=sx.*(yN./yD);

			y(find(abs(x)<4))=sy;

		endif

		if (length(lx))

		# Polynomials for values with magnitude >= 4
		# These are actually two Padé expansions,
		# one of (Si(x)-pi/2)/(x*cos(x)), the other one
		# of (Si(x)-pi/2)/sin(x)
			
			lx_=lx.*sign(lx);

			fNC=[ -4.94701168645415959931e11
			       4.94816688199951963482e12
			       1.00795182980368574617e13
			       4.20968180571076940208e12
			       4.20968180571076940208e12
			       6.40533830574022022911e11
			       4.33736238870432522765e10
			       1.43073403821274636888e9
			       2.37750310125431834034e7
			       1.96396372895146869801e5
			       7.44437068161936700618e2
			       1. ];

			fDC=[ 1.11535493509914254097e13
			      1.43468549171581016479e13
			      5.06084464593475076774e12
			      7.08501308149515401563e11
			      4.58595115847765779830e10
			      1.47478952192985464958e9
			      2.41535670165126845144e7
			      1.97865247031583951450e5
			      7.46437068161927678031e2
			      1. ];

			gNC =[ -1.36517137670871689e12
			        6.43291613143049485e12
			        1.81004487464664575e13
			        7.57664583257834349e12
			        1.09049528450362786e12
			        6.83052205423625007e10
			        2.06297595146763354e9
			        3.12557570795778731e7
			        2.35239181626478200e5
			        8.1359520115168615e2
			        1. ];
    
			gDC=[ 3.99653257887490811e13
			      4.01839087307656620e13
			      1.17164723371736605e13
			      1.39866710696414565e12
			      7.87465017341829930e10
			      2.23355543278099360e9
			      3.26026661647090822e7
			      2.40036752835578777e5
			      8.19595201151451564e2
			      1. ];
    
			lx2=1./lx_.^2;

			f=1./lx_.*polyval(fNC,lx2)./polyval(fDC,lx2);
			g=lx2.*polyval(gNC,lx2)./polyval(gDC,lx2);
			
			ly=(pi/2)-f.*cos(lx_)-g.*sin(lx_);
			
			ly.*=sign(lx);
			
			y(find(abs(x)>=4))=ly;

		endif
	else
		y = zeros(size(x));
 		if prod(size(x)) < 101
  			for k = 1:prod(size(x))
   				y(k) = sum(besselj([0:100]+0.5,(x(k)/2)).^2);
  			endfor
			y = y.*pi;
 		else
  			for k=0:100
   				y += besselj(k+0.5,x/2).^2;
  			endfor
			y = y.*pi;
 		endif 
	endif
endfunction;