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