EvolReproEffort = function(asa,alphasa,af,alphaf,mutarate,sdev,iniadults,theta,maxspec,nostep,iffig) {
# asa - multiplicative constant of adult mortality
# alphasa - power of adult mortality
# af - multiplicative constant of adult fecundity 
# alphaf - power of adult fecundity 
# mutarate - speciation rate
# sdev - standard deviation of the distribution from which the theta of the new species is drawn
# iniadults - initial proportions of adults
# theta - reproductive efforts for the initial species
# maxpec - technical limit of maximum species number
# nostep - number of steps
# iffig - whether to make figures

  eps = 0.001

  if (iffig) {
    graphics.off()
    a<-system("wmic desktopmonitor get screenwidth",intern=T) 
    scrw = as.numeric(a[2])  
  }

  outad = matrix(0,nrow=nostep+1,ncol=maxspec)
  fe = af * theta ^ alphaf
  sa = asa * (1 - theta) ^ alphasa
  adults = iniadults / sum(iniadults,na.rm=T)
  outad[1,1:length(adults)] = adults 
  besttheta = NULL
  wtheta = NULL
  for (i in 1:nostep) {
    newadults = adults * sa + adults * fe
    adults = newadults / sum(newadults,na.rm=T)
    fe = ifelse(adults<eps,0,fe)
    sa = ifelse(adults<eps,0,sa)
    theta = ifelse(adults<eps,NA,theta)
    adults = ifelse(adults<eps,0,adults)
    adults = adults / sum(adults,na.rm=T)
    outad[i+1,1:length(adults)] = adults 
    if (runif(1) < mutarate) {
      if (sum(adults==0) > 0)
        ind = min(which(adults==0))
      else
        ind = length(adults) + 1
      tmp = theta[!is.na(theta)]
      ispec = ceiling(runif(1)*length(tmp))
      mn = tmp[ispec]
      theta[ind] = min(max(rnorm(1,mn, sdev),0),1)  
      adults[ind] = eps * 3
      fe[ind] = af * theta[ind] ^ alphaf
      sa[ind] = asa * (1 - theta[ind]) ^ alphasa
      adults = adults / sum(adults,na.rm=T)
    }
    ii = which(adults == max(adults,na.rm=T))
    besttheta = c(besttheta,theta[ii])
    wtheta = c(wtheta,sum(theta*adults,na.rm=T)/(sum(adults,na.rm=T)))
  }

  dev.new(2,xpos=10,ypos=10,width=6,height=6)
  matplot(outad,col=rainbow(length(adults)),type="l",lwd=2,lty=1,xlab="Time step",ylab="Species proportion")
  text(format(theta,digits=2),y=adults,x=nostep+1)

  dev.new(3,xpos=10+scrw/2,ypos=10,width=6,height=6)
  matplot(1:nostep,cbind(besttheta,wtheta),type="l",lty=1,lwd=2,xlab="Time step",ylab="Reproductive effort")
  legend("topleft",col=c("black","red"),lty=1,lwd=2,legend=c("Dominant species","Community mean"))

  return(list(adults[adults>0],theta[adults>0]))
  
}



PlotFitnessTheta = function(asa,alphasa,af,alphaf) {
  theta = seq(0,1,by=0.01)
  fe = af * theta ^ alphaf
  sa = asa * (1 - theta) ^ alphasa
  fitness = fe+sa
  matplot(theta,cbind(fe,sa,fitness),lty=1,lwd=2,col=c("blue","green","red"),type="l",
         ylim=c(0,max(fe,sa,fitness)+0.4),xlab="Reproductive effort",ylab="Fitness components")
  legend("top",lty=1,pch=20,col=c("blue","green","red"),
       legend=c("Fecundity","Adult survival","Fitness"),ncol=2)
  points(theta[fitness==max(fitness)],fitness[fitness==max(fitness)],pch=20,cex=2,col="red")
  text(paste(theta[fitness==max(fitness)]),x=theta[fitness==max(fitness)],y=fitness[fitness==max(fitness)],pos=3)
}