# Gauss peak curve function for calculating daily precipitation
# probability and mean rain volume
# numbering of months begins with August = 0
gauss <- function (day, amplitude, location, width, shape=2) 
{	day = (day + 182 - location) %% 365 + location - 182
	G = amplitude*exp(-(day-location)^shape/(2*width^2))
	return(G)
}

# Calculation of parameters based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
rmp<-function (MAP, relCh=0) 
{	A = 0.13 + 0.00041 * MAP # rain occurance
	A = A*(1-relCh +1.33*relCh^2-(0.61+1.57*A)*relCh^3)
	L = 177
	W = 52 + 0.007 * MAP
	a = -14 + 4 * log(MAP) # rain volume
	a = a * (1+relCh)
	l = 170
	w = 10488
	return(data.frame(P.amplitude=A, P.location=L, P.width=W, V.amplitude=a, V.location=l, V.width=w))
}

############# Production of stochastic time series ############################
# Parameters can be entered directly or as 'scenario' calculated by function  #
# 'rmp'.                                                                      #
# ReGen returns a matrix with 365 colummns and one row for each year.         #

ReGen<-function (years, P.amplitude=0, P.location=0, P.width=0, R.amplitude=0, R.location=0, R.width=0, scenario=0) 
{	DM = matrix(rep(1:365, years),nrow=years, byrow=T)
	if(sum(scenario)>0)
	{	P.amplitude=as.numeric(scenario[1])
		P.location=as.numeric(scenario[2])
		P.width=as.numeric(scenario[3])
		R.amplitude=as.numeric(scenario[4])
		R.location=as.numeric(scenario[5])
		R.width=as.numeric(scenario[6])
	}

	# probability of a rainy day
	G<-gauss(DM, P.amplitude, P.location, P.width)
	RT<-ifelse(runif(DM)0.05, 1, 0)
	# daily mean rain volume:
	RV<- gauss(DM, R.amplitude, R.location, R.width, 4)
	# actual rain volume:
	R<- RT * rexp(DM, 1/RV)
	return(R)
}