#r0 function
SIRR0<-function(b,k,g,p,beta,delta,c,S0,omega,aH,gammaH,h){
R0<-(b*k*(g*(-1+ p)^2+p)*beta*delta+b*c*(-g*(-1+p)^3+p)*S0*omega+(-1+g)*(-1+p)^2*(aH*h+gammaH-h*gammaH)*(-k*beta*delta+c*(-1+p)*S0*omega))/(b*k*(b*g-(-1+g)*(aH*h+gammaH-h*gammaH))*delta)
} 
#c values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

cval<-seq(0,50,by=5)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=cval,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~cval,pch=19)

#g values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

gval<-seq(0,1,by=.1)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=gval,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~gval,pch=19)

#h values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

hval<-seq(0,1,by=0.1)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=hval,omega=0,delta=1,k=1000000)
plot(R0vals~hval,pch=19)

#omega values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

omegaval<-seq(0,1,by=.1)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=omegaval,delta=1,k=1000000)
plot(R0vals~omegaval,pch=19)

#k values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

kval<-seq(0,1000000,by=100000)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=kval)
plot(R0vals~kval,pch=19)

#aH values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

aHval<-seq(0,20,by=2)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/aHval,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~aHval,pch=19)

#b values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

bval<-seq(0,5,by=.5)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/bval,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~bval,pch=19)

#delta values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

deltaval<-seq(0,5,by=.5)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1/deltaval,k=1000000)
plot(R0vals~deltaval,pch=19)

#S0 values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

S0val<-seq(0,100000,by=10000)
R0vals<-SIRR0(beta=13,S0=S0val,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~S0val,pch=19)

#gammaH values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

gammaHval<-seq(0,20,by=2)
R0vals<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/gammaHval,aH=1/2.7,h=0.00082,omega=0,delta=1,k=kval)
plot(R0vals~gammaHval,pch=19)

#beta values
R0<-SIRR0(beta=13,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)

betaval<-seq(0,20,by=2)
R0vals<-SIRR0(beta=betaval,S0=100000,p=.2,c=0,b=1/0.77,g=0,gammaH=1/2.7,aH=1/2.7,h=0.00082,omega=0,delta=1,k=1000000)
plot(R0vals~betaval,pch=19)