# Sine integral (two routines)

#pragma rtGlobals=1     // Use modern global access method.
// Calculate the sine integral based on formulas from Abramowitz and Stegun.
// See equations 5.2.14 and 5.2.38 and 5.2.39. Si(x)=Integral from 0 to x of sin(x')/x'.
// Peak relative error is approximately 5e-7. This routine is approximately 10x faster
// than the higher precision routine below.
// Tested using IgorPro 6.02 in a Windows XP box.
// Version 11-30-07-01 C. DeJoseph, Jr.
function Si(xin)
variable xin

variable x2,ff,gg,c8,c9,tt,t,i
variable p2=1.57079632679489
variable term1,term2
x2=xin*xin
if(x2>=3.8)
term1=38.102495+x2*(335.677320+x2*(265.187033+x2*(38.027264+x2)))
term2=xin*(157.105423+x2*(570.236280+x2*(322.624911+x2*(40.021433+x2))))
ff=term1/term2
term1=21.821899+x2*(352.018498+x2*(302.757865+x2*(42.242855+x2)))
term2=x2*(449.690326+x2*(1114.978885+x2*(482.485984+x2*(48.196927+x2))))
gg=term1/term2
return p2*sign(xin)-ff*cos(xin)-gg*sin(xin)
else
c8=1.
c9=1.
tt=xin
t=0
i=1
do
t=t+tt/(c8*c9)
tt=-tt*x2
c8=2.*i+1.
c9=c9*c8*2.*i
i+=1
while(i<7)
return t
endif
end
// High precision version of the sine integral based on direct numerical integration
// of the sinc function. Si(x)=Integral from 0 to x of sin(x')/x'. The integration uses
// gaussian quadrature. Peak relative error is approximately 2e-15 and speed is
// about 10% of the lower precision routine above.
// Tested using IgorPro 6.02 in a Windows XP box.
// Version 11-30-07-01 C. DeJoseph, Jr.
function Si_hp(xin)
variable xin
return integrate1D(specialsinc, 0., xin,2)
end
function specialsinc(xin)
variable xin
return sinc(xin)
end

Forum

Support

Gallery