Download: Windows application, OS X Console application, C++ code
# ========= NEW VERSION: ReGen 2.0 ================================
# Parameters derived from stations in Israel and western Jordan
# Time units of parameters can be set to days or months
# Simulation of hourly values based on analysis of data from four field stations
# More flexible Gauss function (variable 'shape' parameter)
# 
# Martin Köchy, 2008-09-02
plot of ReGen 2.0 vs ReGen 1.0 predictions

# Gauss peak curve function for calculating daily precipitation
# probability and mean rain volume
# NUMBERING OF MONTHS AND DAYS STARTS WITH August 1 = 0

# The Gauss curve for time measured in DAYS
gauss.d <- function (day, amplitude, location, width, shape=2) 
{	day = (day + 182 - location) %% 365 + location - 182
	G = amplitude*exp(-abs(day-location)^shape/(2*width^shape))
	return(G)
}

# The Gauss curve for time measured in MONTHS
gauss.m <- function (month, amplitude, location, width, shape=2) 
{	month = (month + 6 - location) %% 12 + location - 6
	G = amplitude*exp(-abs(month-location)^shape/(2*width^shape))
	return(G)
}

# Generic Gauss curve, picks the appropriate function
gauss <- function (time.units, amplitude, location, width, shape=2)
{
	if (max(time.units)>300) return(gauss.d(time.units, amplitude, location, width, shape) else 
	return(gauss.m(time.units, amplitude, location, width, shape)
	}


# Calculation of PARAMETERS based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL and eastern JORDAN.
# THIS IS THE VERSION FOR DAYS
rmp.ILWJO.d<-function(MAP, relCh=0) {
	if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
	warning("This function has been validated only for MAP between 90 and 800 mm 
	and a factor between -0.3 and +0.3 for the rain volume") }
		MAP = (MAP+6)/1.03 # fine-tuning of output
		A = 0.1082772 + 0.0004189*MAP # R2 = 0.78, n=64
		A2 = A*1*(1-relCh + (0.3615545 - 0.0579928*log(MAP))*relCh + relCh^2 -0.4*relCh^3 )
		L = 182.35245 - 26.135966*A   # R2 = 0.44, n=64
		W = 31.24488 + 4.7929879*log(MAP) # R2 = 0.43, n=63
		S = -0.360935 -2.55397*A + 0.05713082*W + 0.000465807*MAP # R2 = 0.82, n=64

		a = -5.3822417 + 0.01660176* MAP + 4.17701974*log(MAP) -40.243339*A -0.0666435*W # R2 = 0.94, n=64
		a2 = a*(1+relCh)
		l = 191.510119 - 48.62739*A # R2 = 0.43, n=62
		w = -3.501397 + 0.0231263*MAP - 2.1062518*a + 0.67337071*l # R2 = 0.34, n=63
		s = 72.8194938 -0.3949018*W + 6.20975069*S - 1.4508933*w + 0.00861436*w^2 # R2 = 0.90, n=54
			
 return(data.frame(P.amplitude=A2, P.location=L, P.width=W, P.shape=S, V.amplitude=a2, V.location=l, V.width=w, V.shape=s))	}
 

# =========== PUBLISHED VERSION: ReGen 1.0 Parameters ==============================================
# Calculation of parameters based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL only. Corresponds to the published version.
# THIS IS THE VERSION FOR DAYS
rmp.iIL.d<-function (MAP, relCh=0) 
{		if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
			warning("This function has been validated only for MAP between 90 and 800 mm 
            and a factor between -0.3 and +0.3 for the rain volume") }

	A = 0.130034 + 0.000408 * 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 = 102.4109 # differs from published version to reflect that 'shape' in the Gauss function 
	             # is now a variable but produces identical results
	return(data.frame(P.amplitude=A, P.location=L, P.width=W, P.shape=2, V.amplitude=a, V.location=l, V.width=w, V.shape=4))
}

# Calculation of parameters based on mean annual precipitation (MAP)
# and change of daily mean rain (relCh)
# Parameters derived from stations in ISRAEL only. 
# THIS IS THE ORIGINAL VERSION FOR MONTHS.
rmp.iIL.m<-function (MAP, relCh=0) {
		if(min(MAP) < 80 | max(MAP) > 960 | min(relCh) < -0.3 | max(relCh) > 0.3) {
		warning("This function has been validated only for MAP between 90 and 800 mm 
		and a factor between -0.3 and +0.3 for the rain volume") }

 A = 0.130034 + 0.000408 * MAP # rain occurance
 A = A*(1-relCh +1.33*relCh^2-(0.61+1.57*A)*relCh^3)
 L = 5.735
 W = 1.704746 + 0.000242 * MAP
 a = -13.7178 + 4.0508301 * log(MAP) # rain volume
 a = a * (1+relCh)
 l = 5.586
 w = 11.336 # differs from published version to reflect that 'shape' in the Gauss function is now 
            # a variable but produces identical results
 return(data.frame(P.amplitude=A, P.location=L, P.width=W, V.amplitude=a, V.location=l, V.width=w))}

# ======================End ReGen 1.0 Parameters ====================================================
 

# Wrapper function to produce the parameters for the stochastic time series. 
# The user can choose time units (default is 'days')
# and the region (interior Israel or Jordan River valley) used for the parameterization;
# 'interior Israel' corresponds to the publsihed version.

rmp<-function(MAP, relCh=0, region='jordanvalley', time.units="d") {
	switch(region,
		'jordanvalley' = {
			switch(time.units,
			'd' = rmp.ILWJO.d(MAP, relCh), 
			warning("No function for this time units for this region"))},
		'interiorIsrael' = {
			switch(time.units,
			 'd' = rmp.iIL.d(MAP, relCh),
			 'm' = rmp.iIL.m(MAP, relCh),
			 warning("No function for this time unit for this region"))},
		warning("No function for this region."))
}
		
	
############# 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, P.shape=2, 
                        amplitude=0, location=0, width=0, shape=4, scenario=0) 
{	DM = matrix(rep(1:365, years),nrow=years, byrow=T)
	if(is.data.frame(scenario))
	{	P.amplitude=as.numeric(scenario$P.amplitude)
		P.location=as.numeric(scenario$P.location)
		P.width=as.numeric(scenario$P.width)
		P.shape=as.numeric(scenario$P.shape)
		V.amplitude=as.numeric(scenario$V.amplitude)
		V.location=as.numeric(scenario$V.location)
		V.width=as.numeric(scenario$V.width)
		V.shape=as.numeric(scenario$V.shape)
	}

	# probability of a rainy day
	G<-gauss.d(DM, P.amplitude, P.location, P.width, P.shape)
	RT<-ifelse(runif(DM)0.05, 1, 0)
	# daily mean rain volume:
	RV<- gauss.d(DM, V.amplitude, V.location, V.width, V.shape)
	# actual rain volume:
	R<- RT * ifelse(RV>0.01,rexp(DM, 1/pmax(RV,0.01)),0) # pmax: prevent division by zero
		return(R)
}


# =========   Break down the daily rain amount into hourly amounts ===============================
ReGen.Tag <- function(daily.amount) {
	# derived from measurements at 4 Israeli field sites
	a = 0.07*daily.amount # a = p/(1-p), p = rainy hours/d
	h = round(a/(1+a)*24)
	d1 = seq(0, floor((24-h)/2))
	d2 = 24-d1
	d = c(d1, seq(daily.amount/h, h), seq(0, d2))
	}