function(v.scale = 0.00015173, gamma = 3.8724, peak = 13.669)
{
#Program compute.gamma(v.scale, gamma, peak)
#
#Default argumentes based on gamma function fitted to empirical data (Basia).
#
	x <- 1:(aa - 10) - 0.5
	h.scale <- peak/(gamma - 1)
	y <- v.scale * x^(gamma - 1) * exp( - x/h.scale)
	phi <- c(rep(0, times = 10), y)
	phi <- phi[1:aa]	#truncate if necessary
	phi <- phi/max(phi)	#normalize to make maximum value equal to one
	names(phi) <- 1:aa - 0.5
	phi <- round(phi, 4)
	write.ascii(compute.phi)
	phi
}