#FIRST THE CODE FOR FIGURE 1
#2 alleles in x with a reduction due to males not mated

sr <- function(s,c,r){ 3 a function that gives eq 2, given s, c and r
   a<-(1-c)*(1-r)
   b<- r+c-c*r
   x<-2   #just draw it for two alleles
   ab<-(2*a*s + b*x + sqrt(x*(4*a*s+b*b*x)))/(2*a*(x-1))  #eq2 in the same variables
   ab} #send the answer back to the request


r <- (50:100)/100
s <- 1
MF2s0 <- matrix(0,nrow = 51, ncol = 2)
MF2s5 <- matrix(0,nrow = 51, ncol = 2)
MF2s1 <- matrix(0,nrow = 51, ncol = 2)
MFi <- matrix(0,nrow = 51, ncol = 2)
for (cc in 1:2){
   if (cc==1) c<- 0            #for two levels of constrained oviposition = 0, 0.387
   if (cc==2) c<- 0.387  
   for (i in 1:51){
      MF2s0[i,cc] <- sr(0,c,r[i])   # for three levels of s = 0, 0.5 and 1
      MF2s5[i,cc] <- sr(0.5,c,r[i])
      MF2s1[i,cc] <- sr(1,c,r[i]) 
      MFi[i,cc] <- (c+(1-c)*r[i])/((1-c)*(1-r[i]))  #eq 1 correct
      }
   }

par(mar=c(5.5,5,0.5,5.5))
plot(r,MF2s5[,2],xlab ="",ylab = "Expected male:female ratio",type="l",lwd=2,xlim <- c(0.5,1),ylim<-c(0,12),col="red",cex=1.4,las=2,cex.lab=1.4, cex.axis=1.3)
title(xlab=expression(Maternal~investment~"in"~sons~(italic(r))), line=4.2, cex.lab=1.4) 
lines(r,MF2s5[,1],lty = "dashed",lwd="2",col="red")

lines(r,MFi[,2],lty = "solid",lwd="2",col="blue")
lines(r,MFi[,1],lty = "dashed",lwd="2",col="blue")

lines(r,MF2s0[,2],lty = "solid",lwd="2",col="brown")
lines(r,MF2s0[,1],lty = "dashed",lwd="2",col="brown")

lines(r,MF2s1[,2],lty = "solid",lwd="2",col="green")
lines(r,MF2s1[,1],lty = "dashed",lwd="2",col="green")


abline(h=10,lwd="2",col="grey")
abline(v=0.66,lwd="2",col="grey")
maleprop_label <- seq(0, 12, 2)
axis(4, at=maleprop_label, labels=round(maleprop_label/(1+maleprop_label), 2), las=2, cex.axis=1.3)
box(lwd=2)
mtext(side=4,  cex=1.4, "Male proportion", line=4)



####SECOND THE CODE FOR FIGURE S1

#now for diploid male frequency
dip <- function(s,c,r){
   a<-(1-c)*(1-r)
   b<- r+c-c*r
   x<-2
   df<-(2*a*s + b*x - sqrt(x*(4*a*s+x*b*b)))/(2*a*(s-x))   #eq 3
   df}   # pass back the frequency of diploid males


r <- (50:100)/100
MF2s0 <- matrix(0,nrow = 51, ncol = 2)
MF2s5 <- matrix(0,nrow = 51, ncol = 2)
MF2s1 <- matrix(0,nrow = 51, ncol = 2)
for (cc in 1:2){
   if (cc==1) c<- 0
   if (cc==2) c<- 0.387  
   for (i in 1:51){
      MF2s0[i,cc] <- dip(0,c,r[i]) 
      MF2s5[i,cc] <- dip(0.5,c,r[i])
      MF2s1[i,cc] <- dip(1,c,r[i]) 
      }
   }
par(mar=c(5.5,6,0.5,0.5))
plot(r,MF2s5[,2],xlab ="",ylab = "",type="l",lwd=2,xlim <- c(0.5,1),ylim<-c(0,0.25),col="red",cex=1.4,las=2,cex.lab=1.4, cex.axis=1.3)
title(xlab=expression(Maternal~investment~"in"~sons~(italic(r))), line=4.2, cex.lab=1.4) 
title(ylab = "fraction of males that are diploid",mgp=c(4.5,1,0), cex.lab=1.4)

lines(r,MF2s5[,1],lty = "dashed",lwd="2",col="red")

lines(r,MF2s0[,2],lty = "solid",lwd="2",col="brown")
lines(r,MF2s0[,1],lty = "dashed",lwd="2",col="brown")

lines(r,MF2s1[,2],lty = "solid",lwd="2",col="green")
lines(r,MF2s1[,1],lty = "dashed",lwd="2",col="green")
box(lwd=2)


s<-0.5
x<-2
2*(-1+c)*(-1+r)*s + x*(r+c-r*c) 
sqrt(x*(4*(-1+c)*(-1+r)*s+x*((c + r - c*r)^2)))