unequalSimpson <- function(nodes, values)
{
	#this is a specialized simpson which assumes that each 
	#interval has a midpoint at which the function has been 
	#evaluated
	
	if(length(nodes) != length(values))
	{stop("Simpson: number of nodes unequal to number of values")}
	
	len <- length(nodes)
	pieces <- (len - 1)/2
	if(round(pieces) != pieces)
	{stop("Simpson: number of nodes must equal 2*n + 1")}

	#calculates simpsons rule for all pieces defined by three
	#successive points starting with the first one
	simpEvals <- vector(mode = "numeric", length = pieces)

	for(i in seq(1,pieces))
	{
		idx0 <- 2*i - 1
		idx1 <- 2*i
		idx2 <- 2*i + 1
		
		y0 <- values[idx0]
		y1 <- values[idx1]
		y2 <- values[idx2]
		
		delX1 <- nodes[idx1] - nodes[idx0]
		delX2 <- nodes[idx2] - nodes[idx1]	
	
		simpEvals[i] <- 1/6*(delX1 + delX2)*(
							(2 - delX2/delX1)*y0 +
							(delX1 + delX2)^2/(delX1*delX2)*y1 +
							(2 - delX1/delX2)*y2)
	}
	
	return(simpEvals)
}
