# please change directory
setwd("C:/work/R/competition")

A  <- 10        # initial concentration resource A
B  <- 15        # initial concentration resourceB
X1 <- 1 # initial biomass species X1
X2 <- 1         # initial biomass species X2
As <- 3         # supply point resource A
Bs <- 3 # supply point recource B


d  <- 0.5               # dilution rate

# uptake vectors
# coexistence possible or not ?
if (T){ # T (true) or F (false)
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

# parameters species one and two
mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

# calculation of the zero net growth isoclines
A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

# intersection of ZNGI's
A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value

# time vector for calulation
dt <- 0.05              # time step
t.max <- 20             # time max
time <- seq(0,t.max,dt)         # time vector

# colors for plot
col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2


# plot
graphics.off()  # clear graphic table
windows(10,10,xpos=1)   # window 1 (nutrients, ZNGI, supply point)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
# nutrient point

windows(4,4)    # window 2 (biomass)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)   # break in seconds
# loop
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))    # actual growth rate spec. 1
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))    # actual growth rate spec. 2

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2     # change rate nutrient A
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2     # change rate nutrient B
dX1 <- mu1*X1 - d*X1                                    # change rate species 1
dX2 <- mu2*X2 - d*X2                                    # change rate species 2

# change of the state variables
A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)      # switch to window 1
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
 # plot a new point
Sys.sleep(0.05) # break

dev.set(3)      # switch to window 2
points(t,X1,col=col1,cex=0.5)   # plot a new point
points(t,X2,col=col2,cex=0.5)   # plot a new point

}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 3         # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 20             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)


windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 6         # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 20             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 8         # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 50             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 10        # initial B
X1 <- 2 # initial X1
X2 <- 2         # initial X2
Y  <- 0.1       # initial Y
As <- 18        # supply point A
Bs <- 13        # supply point B


d  <- 0.5               # dilution rate

a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2

mu.max1 <- 1    # max growth rate species 1
mu.max2 <- 3    # max growth rate species 2
ka1 <- 2                # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 2                # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 100            # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point


windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")


for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

ing1 <- 0.0*X1/(X1+2)*Y
ing2 <- 4  *X2/(X2+2)*Y

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1 - ing1
dX2 <- mu2*X2 - d*X2 - ing2
dY  <- (ing1 + ing2) - 2*Y

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt
Y  <- Y  + dY*dt

A1. <- ka1 * (m1+ing1) / (mu.max1 - (m1+ing1))  # treshold value
resource A spec 2
B1. <- kb1 * (m1+ing1) / (mu.max1 - (m1+ing1))  # treshold value
resource A spec 2
A2. <- ka2 * (m2+ing2/X2) / (mu.max2 - (m2+ing2/X2))    # treshold
value resource A spec 2
B2. <- kb2 * (m2+ing2/X2) / (mu.max2 - (m2+ing2/X2))    # treshold
value resource A spec 2


A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value



dev.set(2)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point

points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2),pch=19)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
points(t,Y ,cex=0.5)

}

setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 10        # initial B
X1 <- 5 # initial X1
X2 <- 5         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (F){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 1    # max growth rate species 2
ka1 <- 10               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 10               # half saturation constante resource B species 1
kb2 <- 2                # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 100            # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point


windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,20),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")


for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))



dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2


A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

if (t %% 5 == 1){
#A <- A + runif(1)*10
#B <- B + runif(1)*10
X1 <- X1 * 0.3 + runif(1) - 0.5
X2 <- X2 * 0.3 + runif(1) - 0.5
}

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=0.25,pch=19)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 10        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 50             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point


windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")


for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 11        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 50             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 25             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point


windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")


for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (T){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 25             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 15        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (F){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 25             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,20),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 15        # initial A
B  <- 10        # initial B
X1 <- 1 # initial X1
X2 <- 1         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (F){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 25             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,20),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 15        # initial A
B  <- 10        # initial B
X1 <- 1 # initial X1
X2 <- 3         # initial X2
As <- 18        # supply point A
Bs <- 18        # supply point B


d  <- 0.5               # dilution rate

#coexistence ?
if (F){
a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2
}else{
a1 <- 0.5               # uptake rate resource A by species 1
a2 <- 1.5               # uptake rate resource A by species 2
b1 <- 1.5               # uptake rate resource B by species 1
b2 <- 0.5               # uptake rate resource B by species 2
}

mu.max1 <- 2    # max growth rate species 1
mu.max2 <- 2    # max growth rate species 2
ka1 <- 15               # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 15               # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 25             # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,20),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(10)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1
dX2 <- mu2*X2 - d*X2

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt

dev.set(2)
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2)/2,pch=19)
Sys.sleep(0.05)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
}


setwd("C:/work/R/competition")

# load initial values/parameters
#source("compini1")

A  <- 10        # initial A
B  <- 10        # initial B
X1 <- 2 # initial X1
X2 <- 2         # initial X2
Y  <- 0 # initial Y
Y <- 0.1
As <- 18        # supply point A
Bs <- 13        # supply point B


d  <- 0.5               # dilution rate

a1 <- 1.5               # uptake rate resource A by species 1
a2 <- 0.5               # uptake rate resource A by species 2
b1 <- 0.5               # uptake rate resource B by species 1
b2 <- 1.5               # uptake rate resource B by species 2

mu.max1 <- 1    # max growth rate species 1
mu.max2 <- 3    # max growth rate species 2
ka1 <- 2                # half saturation constante resource A species 1
ka2 <- 2                # half saturation constante resource A species 2
kb1 <- 2                # half saturation constante resource B species 1
kb2 <- 2                # half saturation constante resource B species 2
m1 <- 0.5               # mortality rate species 1
m2 <- 0.5               # mortality rate species 1

A1. <- ka1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 1
A2. <- ka2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2
B1. <- kb1 * m1 / (mu.max1 - m1)        # treshold value resource A spec 2
B2. <- kb2 * m2 / (mu.max2 - m2)        # treshold value resource A spec 2

A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value


dt <- 0.05              # time step
t.max <- 100            # time max
time <- seq(0,t.max,dt)         # time vector

col1 <- rgb(0.7,0,0)    # colour of spec 1
col2 <- rgb(0,0,0.7)    # colour of spec 2

graphics.off()
windows(10,10,xpos=1)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
lines(c(A..,A.. + 1000*a1),c(B..,B.. + 1000*b1),col=col1,lwd=2,lty=2)
 # consumption vector 1
lines(c(A..,A.. + 1000*a2),c(B..,B.. + 1000*b2),col=col2,lwd=2,lty=2)
 # consumption vector 2

points(A..,B..,cex=1.5)         # intersection point
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point
points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2),pch=19)

windows(4,4)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,t.max),ylim=c(0,10),xaxs="i",yaxs="i",xlab =
"time", ylab="biomass")

Sys.sleep(20)
for (t in time) {

mu1 <- mu.max1 * min( A/(A+ka1) , B/(B+kb1))
mu2 <- mu.max2 * min( A/(A+ka2) , B/(B+kb2))

ing1 <- 0.0*X1/(X1+2)*Y
ing2 <- 4  *X2/(X2+2)*Y

dA  <- d * (As - A) - a1*mu1*X1 - a2*mu2*X2
dB  <- d * (Bs - B) - b1*mu1*X1 - b2*mu2*X2
dX1 <- mu1*X1 - d*X1 - ing1
dX2 <- mu2*X2 - d*X2 - ing2
dY  <- (ing1 + ing2) - 2*Y

A  <- A  + dA*dt
B  <- B  + dB*dt
X1 <- X1 + dX1*dt
X2 <- X2 + dX2*dt
Y  <- Y  + dY*dt

A1. <- ka1 * (m1+ing1) / (mu.max1 - (m1+ing1))  # treshold value
resource A spec 2
B1. <- kb1 * (m1+ing1) / (mu.max1 - (m1+ing1))  # treshold value
resource A spec 2
A2. <- ka2 * (m2+ing2/X2) / (mu.max2 - (m2+ing2/X2))    # treshold
value resource A spec 2
B2. <- kb2 * (m2+ing2/X2) / (mu.max2 - (m2+ing2/X2))    # treshold
value resource A spec 2


A.. <- max(A1.,A2.)     # intersection point A value
B.. <- max(B1.,B2.)     # intersection point B value



dev.set(2)
par(mar=c(1.5,1.5,0,0)+0.5,mgp=c(0.75,0,0),yaxt="n",xaxt="n",cex.lab=1)
plot(-1,-1,xlim=c(0,20),ylim=c(0,20),xlab = "A", ylab="B",xaxs="i",yaxs="i")
lines(c(A1.,A1.,10000),c(10000,B1.,B1.),col=col1,lwd=2) # ZNGI1
lines(c(A2.,A2.,10000),c(10000,B2.,B2.),col=col2,lwd=2) # ZNGI2
points(As,Bs,cex=3,pch="x",col=rgb(0.25,0.25,0))                # supply point

points(A,B,col=rgb(X1/(X1+X2),0,X2/(X2+X1)),cex=sqrt(X1+X2),pch=19)

dev.set(3)
points(t,X1,col=col1,cex=0.5)
points(t,X2,col=col2,cex=0.5)
points(t,Y ,cex=0.5)

}