bessel.i0 <- function(x)
{
# Modified Bessel function of order zero
# GKS  1 Oct 96
	p1 <- 1
	p2 <- 3.5156229
	p3 <- 3.0899424
	p4 <- 1.2067492
	p5 <- 0.2659732
	p6 <- 0.360768/10
	p7 <- 0.45813/100
	q1 <- 0.39894228
	q2 <- 0.1328592/10
	q3 <- 0.225319/100
	q4 <- -0.157565/100
	q5 <- 0.916281/100
	q6 <- -0.2057706/10
	q7 <- 0.2635537/10
	q8 <- -0.1647633/10
	q9 <- 0.392377/100
	y <- (x/3.75)^2
	ax <- abs(x)
	z <- 3.75/ax
	ifelse(abs(x) < 3.75, p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (
		p6 + y * p7))))), (exp(ax)/sqrt(ax)) * (q1 + z * (q2 + z * (q3 + 
		z * (q4 + z * (q5 + z * (q6 + z * (q7 + z * (q8 + z * q9))))))))
		)
}