## Formula for the rectangualr hyperbola

rh <- function(Q, Asat = 40, phi = 0.04, theta = 0.8, Rd = 1){

  num1 <- phi * Q + Asat
  num2 <- num1 ^ 2
  num3 <- 4 * theta * phi * Q * Asat

  den <- 2 * theta

  num <- num1 - sqrt(num2 - num3)

  ans <- num/den - Rd

  ans

}

## Testing the function

Qs <- seq(0,2000)
## Qs <- 0:2000

A.40 <- rh(Qs)
A.50 <- rh(Qs, Asat = 50)
A.60 <- rh(Qs, Asat = 60)

xyplot(A.40 + A.50 + A.60 ~ Qs ,
       auto.key=TRUE,
       xlab = "Quantum flux",
       ylab = "Assimilation")

phi.0.04 <- rh(Qs)
phi.0.05 <- rh(Qs, phi = 0.05)
phi.0.06 <- rh(Qs, phi = 0.06)

xyplot(phi.0.04 + phi.0.05 + phi.0.06 ~ Qs ,
       auto.key=TRUE,
       xlab = "Quantum flux",
       ylab = "Assimilation")

theta.0.3 <- rh(Qs, theta = 0.3)
theta.0.8 <- rh(Qs, theta = 0.8)
theta.1 <- rh(Qs, theta = 1)

xyplot(theta.0.3 + theta.0.8 + theta.1 ~ Qs ,
       auto.key=TRUE,
       xlab = "Quantum flux",
       ylab = "Assimilation")

## Adding reference lines

xyplot(theta.0.3 + theta.0.8 + theta.1 ~ Qs ,
       auto.key=TRUE,
       panel = function(x,y,...){
         panel.xyplot(x,y,...)
         panel.abline(h = 0)
       },
       xlab = "Quantum flux",
       ylab = "Assimilation")