rho     <- 0.5
Sigma1f <- 0.7
Sigma1m <- 0.3
Sigma2f <- 0.3
Sigma2m <- 0.075
Gammaf  <- 0.6
Gammam  <- 0.5

beta <- 15
B <- function(n){
  Bp <- 0
  if ((n[2]+n[4]) > 0){Bp <- beta*n[2]*n[4]/(n[2]+n[4])}
  B <- Bp
  }

A <- function(n){
  {if (n[2] > 0){Ff <- 0.5*B(n)/n[2]}
      else{Ff <- 0}}
  {if (n[4] > 0){Fm <- 0.5*B(n)/n[4]}
      else{Fm <- 0}}
  Ap <- array(0,c(4,4))
  Ap[1,1] <- Sigma1f*(1-Gammaf)
  Ap[1,2] <- rho*Ff
  Ap[1,4] <- rho*Fm
  Ap[2,1] <- Sigma1f*Gammaf
  Ap[2,2] <- Sigma2f
  Ap[3,2] <- (1-rho)*Ff
  Ap[3,3] <- Sigma1m*(1-Gammam)
  Ap[3,4] <- (1-rho)*Fm
  Ap[4,3] <- Sigma1m*Gammam
  Ap[4,4] <- Sigma2m
  A <- Ap
  }

n0 <- c(1,0,1,0)
tmax <- 50

t <- seq(0,tmax,1)
ns <- array(NA,c(4,1))
nn <- array(NA,c(4,1))
Nf <- array(NA,tmax+1)
Nm <- array(NA,tmax+1)
S <- array(NA,c(4,tmax+1))

ns[,1] <- n0
Nf[1] <- ns[1]+ns[2]
Nm[1] <- ns[3]+ns[4]
s <- sum(ns[,1])
for (j in 1:4){S[j,1] <- sum(ns[1:j,1])/s}
for (i in 1:tmax){
    nn <- A(ns) %*% ns
    Nf[i+1] <- nn[1,1]+nn[2,1]
    Nm[i+1] <- nn[3,1]+nn[4,1]
    s <- sum(nn[,1])
    for (j in 1:4){S[j,i+1] <- sum(nn[1:j,1])/s}
    ns <- nn
    }



split.screen(c(2,1))
screen(1)
plot(t,Nf,col="red",type="p",pch=19,
     ylim=c(0,max(Nf,Nm)),xlab="time",ylab="",main="Census",bty="l")
points(t,Nm,col="blue",type="p",pch=19)
points(t,Nm,col="blue",type="l")
points(t,Nf,col="red",type="l")
screen(2)
plot(c(0,tmax),c(0,1),type="n",
     ylim=c(0,1.1),xlab="time",ylab="",main="Structure",bty="l")
polygon(c(0,t,seq(tmax,0,-1)),c(0,S[1,],array(0,tmax+1)),
        col="red",density=15)
polygon(c(0,t,seq(tmax,0,-1)),c(S[1,1],S[2,],S[1,seq(tmax+1,1,-1)]),
        col="red")
polygon(c(0,t,seq(tmax,0,-1)),c(S[2,1],S[3,],S[2,seq(tmax+1,1,-1)]),
        col="blue",density=15,angle=-45)
polygon(c(0,t,seq(tmax,0,-1)),c(S[3,1],S[4,],S[3,seq(tmax+1,1,-1)]),
        col="blue")
close.screen(all=TRUE)

    
print((Nf[tmax+1]+Nm[tmax+1])/(Nf[tmax]+Nm[tmax]))
print(nn/sum(nn))