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