# Parameter für V
A = 2.0
D = 180
S = 2.5

# Verdunst. erz.
V = function(A=2,D=90,S=2.5,Ta){
  Ep = A*sin((Ta+D)/365*2*pi)+S
  return(Ep)
}

#init
Ta = seq(1,365,1)
N = numeric(length(Ta))
Q = numeric(length(Ta))
I = numeric(length(Ta))
A = numeric(length(Ta))
E = numeric(length(Ta))

# Parameter
Imax = 100
Imin = 5
I0 = 20
A0 = 0.5

lambda = 1
eta = 1
beta = 1
Qs = 0

# initialisieren
n = 0

for (ti in Ta) {
  
  n = n+1
  N[n] = runif(1,0.1,4)
  Ep=V(A=2,D=90,S=2.5,ti) 

  # Defizit
  if(n==1){
    I[n] = I0
    Def = Imax-I[n]
    A[n] = A0
    Q[n] = 0
    E[n] = Ep
  } else{
      Def <- Imax*(1-(1-A[n])^{1/beta})-I[n]
  }

  # Änd. der Feuchte
  dS = N[n]-E[n]
  
  if (dS<=Def){
    I[n] = I[n]+dS
  } 
  else{
  }
  
  # neuer Flächenanteil
  A[n] <- 1-(1-(I[n]/Imax))^beta
  
  # Abfluss rechnen
  Q[n] = A[n]*N[n]
  
  # Verd rechnen
  E[n] = Ep*A[n]

}

hydro = data.frame(N,I,A,Q,E)