FitGo_cali <- function(a,InitV,TimeScale,C,psz1900,rmort,Vx){
  fchild<-a[1];vchild<-a[2];rchild<-a[3];w<-a[4];pchild<-a[5];pelderly<-a[6];
  padult<-a[7];n<-a[8];uiscaleC<-a[9];relderly<-a[10];radult<-a[11];x<-a[12];
  fadult<-a[13];felderly<-a[14];vadult<-a[15];CDRscale<-a[16];e<-a[17];
  nelderly<-a[18];velderly<-a[19];uiscaleA<-a[20];uiscaleE<-a[21];neta<-a[22];
  CDRscaleE<-a[23]
  
  C<-0
  InitV<-c(2,0.5,0.02,0.02,0.8,0.07)
  TimeScale<-c(1900,2050)
  psz1900<-530886
  rmort<-0
  
  year1<-TimeScale[1]
  yearend<-TimeScale[2]
  
  if (length(InitV) == 3) {
    for (i in 1:length(InitV)) {
      run<-InitV[1]
      dt<-InitV[2]
      prop<-InitV[3]
    }
  } else {
    for (i in 1:2) {
      run<-InitV[1]
      dt<-InitV[2]
    }
    prop<-InitV[3:6]
  }
  
  
  chiyrs <- 15
  aduyrs <- 50
  yaduyrs <- 40
  eldyrs <- (Mnage - chiyrs - aduyrs)
  
  g<-x 
  hchild<-fchild
  hadult<-fadult
  helderly<-felderly
  
  voldadu <- (vadult + velderly) / 2
  noldadu <- (n + nelderly) / 2
  roldadu <- (radult + relderly) / 2
  CDRscaleO <- (CDRscale + CDRscaleE) / 2
  
  p <- c((rep(pchild, l = chiyrs)), (rep(padult, l = aduyrs)), (rep(pelderly, l = eldyrs)))
  f <- c((rep(fchild, l = chiyrs)), (rep(fadult, l = aduyrs)), (rep(felderly, l = eldyrs)))
  h <- c((rep(hchild, l = chiyrs)), (rep(hadult, l = aduyrs)), (rep(helderly, l = eldyrs)))
  v <- c((rep(vchild, l = chiyrs)), (rep(vadult, l = yaduyrs)), (rep(voldadu, l = (aduyrs - yaduyrs))),(rep(velderly, l = eldyrs)))
  r <- c((rep(rchild, l = chiyrs)), (rep(radult, l = yaduyrs)), (rep(roldadu, l = (aduyrs - yaduyrs))), (rep(relderly, l = eldyrs)))
  n <- c((rep(n, l = chiyrs)), (rep(n, l = yaduyrs)), (rep(noldadu, l = (aduyrs - yaduyrs))), (rep(nelderly, l = eldyrs)))
  
  times <- seq(year1, (yearend + (1 - dt)), dt)
  steps <- length(times)
  
  # Country data
  probt <- 0.1
  mort<-as.matrix(mortage)
  bb <- births[,2]/1000
  cdr<-matrix(0,length(cdrm),Mnage)
  cdr[,1:Mnage]<-cdrm
  suctt<-suctm['China',]
  ps<-pstruc[1:101,'China']
  
  
  if (uiscaleC<0){
    uichild<-uiscaleC*(ui)+ui
    unichild<-uiscaleC*(uni)+uni
  } else {
    uichild<-uiscaleC*(1-ui)+ui
    unichild<-uiscaleC*(1-uni)+uni
  }
  
  if (uiscaleA<0){
    uiadult<-uiscaleA*(ui)+ui
    uniadult<-uiscaleA*(uni)+uni
  } else {
    uiadult<-uiscaleA*(1-ui)+ui
    uniadult<-uiscaleA*(1-uni)+uni
  }
  
  if (uiscaleE<0){
    uielderly<-uiscaleE*(ui)+ui
    unielderly<-uiscaleE*(uni)+uni
  } else {
    uielderly<-uiscaleE*(1-ui)+ui
    unielderly<-uiscaleE*(1-uni)+uni
  }
  
  
  ui<-c((rep(uichild, l=chiyrs)),(rep(uiadult, l=aduyrs)),(rep(uielderly, l=eldyrs)))
  uni<-c((rep(unichild, l=chiyrs)),(rep(uniadult, l=aduyrs)),(rep(unielderly, l=eldyrs)))
  
  ## Generate Vaccine specific parameters
  if(length(Vx)>1){
    assign('vaccine',Vx[1],envir = .GlobalEnv); assign('coverage',Vx[2],envir = .GlobalEnv); 
    assign('coverageM',Vx[3],envir = .GlobalEnv); assign('effI',Vx[4],envir = .GlobalEnv); 
    assign('effD',Vx[5],envir = .GlobalEnv); assign('D',Vx[6],envir = .GlobalEnv); 
    assign('fmass',Vx[7],envir = .GlobalEnv); assign('vxage_child',Vx[8],envir = .GlobalEnv);
    assign('vxage_elder',Vx[9],envir = .GlobalEnv);assign('vxage_type',Vx[10],envir = .GlobalEnv)
    
    source("vaccine_matrix_v3.R")
  } else {d<-matrix(0,steps,Mnage); 
  thetaS<-matrix(0,steps,Mnage);
  thetaL<-matrix(0,steps,Mnage);
  thetaR<-matrix(0,steps,Mnage);
  effI<-0;effD<-0;D<-0;vxage_type<-0;
  coverage<-0;coverageM<-0;fmass<-0;
  vxage_child<-0;vxage_elder<-0;vaccine<-0
  }
  
  ############################################################################################
  ## Build all matrices and initialise
  ## INITIALISING ##
  # Builds all matrices and initialises (both indicating matrices and output ones)
  times<-seq(year1,(yearend+(1-dt)),dt)
  steps<-length(times)
  
  # Susceptible
  S<-matrix(0,steps,Mnage); 
  # Latent
  L<-matrix(0,steps,Mnage);  
  # Infectious
  I<-matrix(0,steps,Mnage);
  # Non-Infectious, 
  NI<-matrix(0,steps,Mnage);  
  # Recovered, 
  R<-matrix(0,steps,Mnage);  
  #newly infected people
  new_infect<-matrix(0,steps,Mnage)
  # New Infectious, 
  new_I<-matrix(0,steps,Mnage); 
  # New non-Infectious
  new_NI<-matrix(0,steps,Mnage); 
  #new reactivation infectious/non infectious
  new_I_react<-matrix(0,steps,Mnage)
  new_I_noconv<-matrix(0,steps,Mnage)
  new_NI_react<-matrix(0,steps,Mnage)
  new_actv<-matrix(0,steps,Mnage)
  new_actv_react<-matrix(0,steps,Mnage)
  new_actv_inf<-matrix(0,steps,Mnage)
  new_actv_chk<-matrix(0,steps,Mnage)
  num_vac<-matrix(0,steps,Mnage)
  new_actvinf<-matrix(0,steps,Mnage)
  
  #New notifications
  new_notif<-matrix(0,steps,Mnage)
  
  # Vaccine matrices: Susceptible, Latent, activeinfectious, active non infectious, Recovered and vaccinated. 
  Sv<-matrix(0,steps,Mnage); 
  Lv<-matrix(0,steps,Mnage);
  Iv<-matrix(0,steps,Mnage);
  NIv<-matrix(0,steps,Mnage);
  Rv<-matrix(0,steps,Mnage); 
  new_Iv<-matrix(0,steps,Mnage); 
  new_NIv<-matrix(0,steps,Mnage);
  new_Iv_noconv<-matrix(0,steps,Mnage); 
  new_notifv<-matrix(0,steps,Mnage); 
  new_actvv<-matrix(0,steps,Mnage)
  
  
  ## Initialising
  if (run==1){
    # If first ever run 
    S0<-(1-prop)*psz1900*(ps/(sum(ps)))
    I0<-prop*psz1900*(ps/(sum(ps)))
    S[1,]<-S0; I[1,]<-I0; 
  } else if (run > 1){
    # If not first run then take input
    S0<-(1-sum(prop))*psz1900*(ps/(sum(ps)))
    I0<-prop[1]*psz1900*(ps/(sum(ps)))
    NI0<-prop[2]*psz1900*(ps/(sum(ps)))
    L0<-prop[3]*psz1900*(ps/(sum(ps)))
    R0<-prop[4]*psz1900*(ps/(sum(ps)))
    S[1,]<-S0; I[1,]<-I0; NI[1,]<-NI0; L[1,]<-L0; R[1,]<-R0;  
  }
  
  # Population size
  psize<-matrix(0,steps,1)
  psize014<-matrix(0,steps,1)
  psize1554<-matrix(0,steps,1)
  psize5564<-matrix(0,steps,1)
  psize65plus<-matrix(0,steps,1)
  psize1559<-matrix(0,steps,1)
  psize1529<-matrix(0,steps,1)
  psize3044<-matrix(0,steps,1)
  psize4559<-matrix(0,steps,1)
  psize60plus<-matrix(0,steps,1)
  psize55plus<-matrix(0,steps,1)
  psize55minus<-matrix(0,steps,1)
  psize0509<-matrix(0,steps,1)
  psize1019<-matrix(0,steps,1)
  psize2029<-matrix(0,steps,1)
  psize3039<-matrix(0,steps,1)
  psize4049<-matrix(0,steps,1)
  psize5059<-matrix(0,steps,1)
  psize6069<-matrix(0,steps,1)
  psize70plus<-matrix(0,steps,1)
  psize5574<-matrix(0,steps,1)
  psize75plus<-matrix(0,steps,1)
  psize1524<-matrix(0,steps,1)
  psize2554<-matrix(0,steps,1)
  psize15plus<-matrix(0,steps,1)
  psizeVX<-matrix(0,steps,1)
  psize15<-matrix(0,steps,1)
  psize60<-matrix(0,steps,1)
  psize6174<-matrix(0,steps,1)
  
  #psize and I for contact matirx
  psizematrix<-matrix(0,steps,16)
  Imatrix<-matrix(0,steps,16)
  BIRTHS<-vector('numeric',steps)
  
  
  #filling in first time step
  psize[1]<-sum(S[1,],L[1,],R[1,],I[1,],NI[1,],Lv[1,],Iv[1,],NIv[1,],Rv[1,])
  psize014[1]<-sum(S[1,1:15],L[1,1:15],R[1,1:15],I[1,1:15],NI[1,1:15],Lv[1,1:15],Iv[1,1:15],NIv[1,1:15],Rv[1,1:15])
  psize1554[1]<-sum(S[1,16:55],L[1,16:55],R[1,16:55],I[1,16:55],NI[1,16:55],Lv[1,16:55],Iv[1,16:55],NIv[1,16:55],Rv[1,16:55])
  psize5564[1]<-sum(S[1,56:65],L[1,56:65],R[1,56:65],I[1,56:65],NI[1,56:65],Lv[1,56:65],Iv[1,56:65],NIv[1,56:65],Rv[1,56:65])
  psize65plus[1]<-sum(S[1,66:Mnage],L[1,66:Mnage],R[1,66:Mnage],I[1,66:Mnage],NI[1,66:Mnage],Lv[1,66:Mnage],Rv[1,66:Mnage],Iv[1,66:Mnage],NIv[1,66:Mnage])
  psize1559[1]<-sum(S[1,16:60],L[1,16:60],R[1,16:60],I[1,16:60],NI[1,16:60],Lv[1,16:60],Rv[1,16:60],Iv[1,16:60],NIv[1,16:60])
  psize1529[1]<-sum(S[1,16:30],L[1,16:30],R[1,16:30],I[1,16:30],NI[1,16:30],Lv[1,16:30],Rv[1,16:30],Iv[1,16:30],NIv[1,16:30])
  psize3044[1]<-sum(S[1,31:45],L[1,31:45],R[1,31:45],I[1,31:45],NI[1,31:45],Lv[1,31:45],Rv[1,31:45],Iv[1,31:45],NIv[1,31:45])
  psize4559[1]<-sum(S[1,46:60],L[1,46:60],R[1,46:60],I[1,46:60],NI[1,46:60],Lv[1,46:60],Rv[1,46:60],Iv[1,46:60],NIv[1,46:60])
  psize60plus[1]<-sum(S[1,61:Mnage],L[1,61:Mnage],R[1,61:Mnage],I[1,61:Mnage],NI[1,61:Mnage],Lv[1,61:Mnage],Rv[1,61:Mnage],Iv[1,61:Mnage],NIv[1,61:Mnage])
  psize55plus[1]<-sum(S[1,56:Mnage],L[1,56:Mnage],R[1,56:Mnage],I[1,56:Mnage],NI[1,56:Mnage],Lv[1,56:Mnage],Rv[1,56:Mnage],Iv[1,56:Mnage],NIv[1,56:Mnage])
  psize55minus[1]<-sum(S[1,1:55],L[1,1:55],R[1,1:55],I[1,1:55],NI[1,1:55],Lv[1,1:55],Rv[1,1:55],Iv[1,1:55],NIv[1,1:55])
  psize0509[1]<-sum(S[1,6:10],L[1,6:10],R[1,6:10],I[1,6:10],NI[1,6:10],Lv[1,6:10],Rv[1,6:10],Iv[1,6:10],NIv[1,6:10])
  psize1019[1]<-sum(S[1,11:20],L[1,11:20],R[1,11:20],I[1,11:20],NI[1,11:20],Lv[1,11:20],Rv[1,11:20],Iv[1,11:20],NIv[1,11:20])
  psize2029[1]<-sum(S[1,21:30],L[1,21:30],R[1,21:30],I[1,21:30],NI[1,21:30],Lv[1,21:30],Rv[1,21:30],Iv[1,21:30],NIv[1,21:30])
  psize3039[1]<-sum(S[1,31:40],L[1,31:40],R[1,31:40],I[1,31:40],NI[1,31:40],Lv[1,31:40],Rv[1,31:40],Iv[1,31:40],NIv[1,31:40])
  psize4049[1]<-sum(S[1,41:50],L[1,41:50],R[1,41:50],I[1,41:50],NI[1,41:50],Lv[1,41:50],Rv[1,41:50],Iv[1,41:50],NIv[1,41:50])
  psize5059[1]<-sum(S[1,51:60],L[1,51:60],R[1,51:60],I[1,51:60],NI[1,51:60],Lv[1,51:60],Rv[1,51:60],Iv[1,51:60],NIv[1,51:60])
  psize6069[1]<-sum(S[1,61:70],L[1,61:70],R[1,61:70],I[1,61:70],NI[1,61:70],Lv[1,61:70],Rv[1,61:70],Iv[1,61:70],NIv[1,61:70])
  psize70plus[1]<-sum(S[1,71:Mnage],L[1,71:Mnage],R[1,71:Mnage],I[1,71:Mnage],NI[1,71:Mnage],Lv[1,71:Mnage],Rv[1,71:Mnage],Iv[1,71:Mnage],NIv[1,71:Mnage])
  psize5574[1]<-sum(S[1,56:75],L[1,56:75],R[1,56:75],I[1,56:75],NI[1,56:75],Lv[1,56:75],Rv[1,56:75],Iv[1,56:75],NIv[1,56:75])
  psize75plus[1]<-sum(S[1,76:Mnage],L[1,76:Mnage],R[1,76:Mnage],I[1,76:Mnage],NI[1,76:Mnage],Lv[1,76:Mnage],Rv[1,76:Mnage],Iv[1,76:Mnage],NIv[1,76:Mnage])
  psize1524[1]<-sum(S[1,16:25],L[1,16:25],R[1,16:25],I[1,16:25],NI[1,16:25],Lv[1,16:25],Rv[1,16:25],Iv[1,16:25],NIv[1,16:25])
  psize2554[1]<-sum(S[1,26:55],L[1,26:55],R[1,26:55],I[1,26:55],NI[1,26:55],Lv[1,26:55],Rv[1,26:55],Iv[1,26:55],NIv[1,26:55])
  psize15plus[1]<-sum(S[1,16:Mnage],L[1,16:Mnage],R[1,16:Mnage],I[1,16:Mnage],NI[1,16:Mnage],Lv[1,16:Mnage],Rv[1,16:Mnage],Iv[1,16:Mnage],NIv[1,16:Mnage])
  psize60[1]<-sum(S[1,61],L[1,61],R[1,61],I[1,61],NI[1,61],Lv[1,61],Rv[1,61],Iv[1,61],NIv[1,61])
  psize15[1]<-sum(S[1,16],L[1,16],R[1,16],I[1,16],NI[1,16],Lv[1,16],Rv[1,16],Iv[1,16],NIv[1,16])
  psize6174[1]<-sum(S[1,62:75],L[1,62:75],R[1,62:75],I[1,62:75],NI[1,62:75],Lv[1,62:75],Rv[1,62:75],Iv[1,62:75],NIv[1,62:75])
  
  #initalise contact matrices
  psizematrix[1,1]<-sum(S[1,1:5],L[1,1:5],R[1,1:5],I[1,1:5],NI[1,1:5],Lv[1,1:5],Rv[1,1:5],Iv[1,1:5],NIv[1,1:5])
  psizematrix[1,2]<-sum(S[1,6:10],L[1,6:10],R[1,6:10],I[1,6:10],NI[1,6:10],Lv[1,6:10],Rv[1,6:10],Iv[1,6:10],NIv[1,6:10])
  psizematrix[1,3]<-sum(S[1,11:15],L[1,11:15],R[1,11:15],I[1,11:15],NI[1,11:15],Lv[1,11:15],Rv[1,11:15],Iv[1,11:15],NIv[1,11:15])
  psizematrix[1,4]<-sum(S[1,16:20],L[1,16:20],R[1,16:20],I[1,16:20],NI[1,16:20],Lv[1,16:20],Rv[1,16:20],Iv[1,16:20],NIv[1,16:20])
  psizematrix[1,5]<-sum(S[1,21:25],L[1,21:25],R[1,21:25],I[1,21:25],NI[1,21:25],Lv[1,21:25],Rv[1,21:25],Iv[1,21:25],NIv[1,21:25])
  psizematrix[1,6]<-sum(S[1,26:30],L[1,26:30],R[1,26:30],I[1,26:30],NI[1,26:30],Lv[1,26:30],Rv[1,26:30],Iv[1,26:30],NIv[1,26:30])
  psizematrix[1,7]<-sum(S[1,31:35],L[1,21:65],R[1,21:65],I[1,21:65],NI[1,21:65],Lv[1,21:65],Rv[1,21:65],Iv[1,21:65],NIv[1,21:65])
  psizematrix[1,8]<-sum(S[1,36:40],L[1,36:40],R[1,36:40],I[1,36:40],NI[1,36:40],Lv[1,36:40],Rv[1,36:40],Iv[1,36:40],NIv[1,36:40])
  psizematrix[1,9]<-sum(S[1,41:45],L[1,41:45],R[1,41:45],I[1,41:45],NI[1,41:45],Lv[1,41:45],Rv[1,41:45],Iv[1,41:45],NIv[1,41:45])
  psizematrix[1,10]<-sum(S[1,46:50],L[1,46:50],R[1,46:50],I[1,46:50],NI[1,46:50],Lv[1,46:50],Rv[1,46:50],Iv[1,46:50],NIv[1,46:50])
  psizematrix[1,11]<-sum(S[1,51:55],L[1,51:55],R[1,51:55],I[1,51:55],NI[1,51:55],Lv[1,51:55],Rv[1,51:55],Iv[1,51:55],NIv[1,51:55])
  psizematrix[1,12]<-sum(S[1,56:60],L[1,56:60],R[1,56:60],I[1,56:60],NI[1,56:60],Lv[1,56:60],Rv[1,56:60],Iv[1,56:60],NIv[1,56:60])
  psizematrix[1,13]<-sum(S[1,61:65],L[1,61:65],R[1,61:65],I[1,61:65],NI[1,61:65],Lv[1,61:65],Rv[1,61:65],Iv[1,61:65],NIv[1,61:65])
  psizematrix[1,14]<-sum(S[1,66:70],L[1,66:70],R[1,66:70],I[1,66:70],NI[1,66:70],Lv[1,66:70],Rv[1,66:70],Iv[1,66:70],NIv[1,66:70])
  psizematrix[1,15]<-sum(S[1,71:75],L[1,71:75],R[1,71:75],I[1,71:75],NI[1,71:75],Lv[1,71:75],Rv[1,71:75],Iv[1,71:75],NIv[1,71:75])
  psizematrix[1,16]<-sum(S[1,76:Mnage],L[1,76:Mnage],R[1,76:Mnage],I[1,76:Mnage],NI[1,76:Mnage],Lv[1,76:Mnage],Rv[1,76:Mnage],Iv[1,76:Mnage],NIv[1,76:Mnage])
  
  
  Imatrix[1,1]<-sum(I[1,1:5],Iv[1,1:5])
  Imatrix[1,2]<-sum(I[1,6:10],Iv[1,6:10])
  Imatrix[1,3]<-sum(I[1,11:15],Iv[1,11:15])
  Imatrix[1,4]<-sum(I[1,16:20],Iv[1,16:20])
  Imatrix[1,5]<-sum(I[1,21:25],Iv[1,21:25])
  Imatrix[1,6]<-sum(I[1,26:30],Iv[1,26:30])
  Imatrix[1,7]<-sum(I[1,31:35],Iv[1,31:35])
  Imatrix[1,8]<-sum(I[1,36:40],Iv[1,36:40])
  Imatrix[1,9]<-sum(I[1,41:45],Iv[1,41:45])
  Imatrix[1,10]<-sum(I[1,46:50],Iv[1,46:50])
  Imatrix[1,11]<-sum(I[1,51:55],Iv[1,51:55])
  Imatrix[1,12]<-sum(I[1,56:60],Iv[1,56:60])
  Imatrix[1,13]<-sum(I[1,61:65],Iv[1,61:65])
  Imatrix[1,14]<-sum(I[1,66:70],Iv[1,66:70])
  Imatrix[1,15]<-sum(I[1,71:75],Iv[1,71:75])
  Imatrix[1,16]<-sum(I[1,76:Mnage],Iv[1,76:Mnage])
  
  # FOI
  lambda<-matrix(0,steps,Mnage);
  
  #initialise lambda
  lambda[1,1:Mnage]<- t(neta * (1-exp(colSums(-(myneta[1:16,1:Mnage]) * z * ((Imatrix[1,1:16])/(psizematrix[1,1:16]))))))
  
  # Number of births,b rate
  birthsnum<-matrix(0,steps,1); bv<-c(); brate<-c();
  
  # Age of death
  ADeaths<-matrix(0,steps,Mnage); 
  # Number of TB deaths, 
  TBDeaths<-matrix(0,steps,Mnage);
  # Number vaccinated
  VX<-matrix(0,steps,3); vaccgive<-matrix(0,1,5); 
  vaccgiveyr<-matrix(0,1,26);
  
  
  # TB Incidence, TB mortality, etc
  TBI<-matrix(0,steps,10);
  TBM<-matrix(0,steps,14);
  TBP<-matrix(0,steps,8);
  TBPb<-matrix(0,steps,9);
  TBPI<-matrix(0,steps,17);
  PSIZEy<-matrix(0,steps,27);
  TBRa<-matrix(0,steps,6);
  TBRi<-matrix(0,steps,6);
  TBRinew<-matrix(0,steps,6);
  TBInew<-matrix(0,steps,6);
  TBN<-matrix(0,steps,9) #添加15-64
  TBAc<-matrix(0,(yearend-year1+1),1)
  TBAc_age<-matrix(0,(yearend-year1+1),6)
  TBMo<-matrix(0,(yearend-year1+1),1)
  NV<-matrix(0,(yearend-year1+1),1)
  TBcases<-matrix(0,steps,14)
  TBdead<-matrix(0,steps,Mnage)
  TBcases_age<-matrix(0,steps,Mnage)
  
  colnames(TBI)<-c("All ages","0-14", "15-54", "55-64", "65+", "55+", "<55","15-24","25-54","15+")
  colnames(TBN)<-c("All ages","0-14", "15-54", "55-64", "65+","55+", "<55","15+", "15-64")
  colnames(TBM)<-c("All ages", "0-14", "15-54", "55-64", "65+", "15-59", "60+", "55+", "<55","15-24","25-54", "15-29", "30-44", "45-59")
  colnames(TBP)<-c("All ages","0-14", "15-29", "30-44", "45-59", "60+", "55+",'30-59')
  colnames(TBPb)<-c("All ages","0-14", "15-29", "30-44", "45-59", "60+", "55+", "15+", "15-44")
  colnames(TBPI)<-c("All ages", "0-14", "15-54", "55-64", "65+", "55+", "5-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70+","6574", "75+", "1524")
  colnames(PSIZEy)<-c("All ages", "0-14", "15-54", "55-64", "65+", "15-59", "15-29", "30-44", "45-59", "60+", "55+", "5-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70+","55-74", "75+","15-24", "25-54","15+",'15','60','61-74')
  
  #others
  I2050<-matrix(0,5,1)
  CDR_av<-matrix(0,(yearend-year1),1)
  #############################################################################################################
  #### RUN MODEL #### 
  for (k in year1:(yearend)){
    
    #if (k < 1997){ z <- 0.2 } else if(k > 2004) { z <- 0.1 } else{z<-0.15}
    #### •••••••••••••••••• Yearly parameters #••••••••••••••••••
    if (k <= 1950){ yr <- 1950 } else { yr <- k }
    # Only start marking years after 2009 
    #yr <- k
    mortdrop<-1950
    
    if (k <= mortdrop){if (rmort < 0) {u<-as.vector(rmort*(mort[1,2:102]) + mort[1,2:102])
    } else {u<-as.vector(rmort*(1-(mort[1,2:102])) + mort[1,2:102])}
    } else {if (rmort < 0) {u<-as.vector(rmort*(mort[1+yr-1950,2:102]) + mort[1+yr-1950,2:102])  
    } else {u<-as.vector(rmort*(1-(mort[1+yr-1950,2:102])) + mort[1+yr-1950,2:102])}
    }
    #### MORTALITY. Runs from age 1 to age 101 (equiv 0-100yo)
    # FIT: Rmort multiplies background death rates. Range [-1,1]
    
    #### CDR & TREATMENT SUCCESS (a proportion of cases that are found and successfully treated)
    if (k <= 1990){ CDR_yr <- 1990 } else { CDR_yr <- k }
    if (k <= 1994){ CoT_yr <- 1994 } else { CoT_yr <- k }
    
    if (CDRscale < 0) {
      CDRscaled <- rep(min(((CDRscale * cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), (chiyrs + yaduyrs))
    } else {
      CDRscaled <- rep(min((CDRscale * (1 - cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), (chiyrs + yaduyrs))
    }
    
    if (CDRscaleO < 0) {
      CDRscaledO <- rep(min(((CDRscaleO * cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), (Mnage - (eldyrs + chiyrs + yaduyrs)))
    } else {
      CDRscaledO <- rep(min((CDRscaleO * (1 - cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), (Mnage - (eldyrs + chiyrs + yaduyrs)))
    }
    
    if (CDRscaleE < 0) {
      CDRscaledE <- rep(min(((CDRscaleE * cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), eldyrs)
    } else {
      CDRscaledE <- rep(min((CDRscaleE * (1 - cdr[1 + CDR_yr - 1990, ]) + cdr[1 + CDR_yr - 1990, ]), 1), eldyrs)
    }
    
    CDR<-c(CDRscaled,CDRscaledO,CDRscaledE)
    
    #calls in 1994 data pre-1994
    CoT<-suctt[1+CoT_yr-1994]
    
    #### BIRTHS 
    ## USES UN BIRTH RATE from 1979
    
    fertdrop<-1950
    
    if (k < fertdrop){ br<-bb[1]
    if (k == year1){B<-round(br*psize[1])
    } else {B<-round(br*psize[((k-year1)*(1/dt))])} 
    } else {B<-round(bb[1+yr-1950]*psize[((k-year1)*(1/dt))])}
    
    BIRTHS[i]<-B
    
    #### END OF YEARLY PARAMETERS  ####
    
    
    #### Start of model run
    if (k > year1){ # Start of year 1 is the initial condition
      ## These are those where the time step is the first of the year
      i = ((1/dt)*(k-year1)+1) # Time zero is first row
      start <- 1 # whether first ts of year. All year dependent parameters must be last years.
      
      #transmission matrix
      lambda[i-1,1:Mnage] <- t(neta * (1-exp(colSums(-(myneta[1:16,1:Mnage]) * z * ((Imatrix[i-1,1:16])/(psizematrix[i-1,1:16]))))))
      
      
      ####= TB model ####
      
      ### First time step of year, aging and transitions of non-vax ###
      
      j = 1; S[i,j]<-B  # Age 1, first time step of the year, all births occur
      
      S[i,2:Mnage] = S[i-1,1:(Mnage-1)] - (u[1:(Mnage-1)]+lambda[i-1,1:(Mnage-1)])*S[i-1,1:(Mnage-1)]*dt 
      L[i,2:Mnage] = L[i-1,1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*(1 - p[1:(Mnage-1)])*(S[i-1,1:(Mnage-1)] + g*R[i-1,1:(Mnage-1)])*dt - (v[1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x + u[1:(Mnage-1)])*L[i-1,1:(Mnage-1)]*dt 
      
      new_infect[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*S[i-1,1:(Mnage-1)]*dt + lambda[i-1,1:(Mnage-1)]*(x*L[i-1,1:(Mnage-1)] + g*R[i-1,1:(Mnage-1)])*dt  
      new_I_react[i,2:Mnage] = v[1:(Mnage-1)]*f[1:(Mnage-1)]*(L[i-1,1:(Mnage-1)])*dt + r[1:(Mnage-1)]*h[1:(Mnage-1)]*R[i-1,1:(Mnage-1)]*dt 
      new_NI_react[i,2:Mnage] =  v[1:(Mnage-1)]*(1 - f[1:(Mnage-1)])*L[i-1,1:(Mnage-1)]*dt + r[1:(Mnage-1)]*(1 - h[1:(Mnage-1)])*R[i-1,1:(Mnage-1)]*dt  
      new_actv_react[i,2:Mnage] = v[1:(Mnage-1)]*(L[i-1,1:(Mnage-1)])*dt + r[1:(Mnage-1)]*R[i-1,1:(Mnage-1)]*dt 
      new_actv_inf[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*S[i-1,1:(Mnage-1)]*dt + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x*(L[i-1,1:(Mnage-1)])*dt + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*g*R[i-1,1:(Mnage-1)]*dt 
      new_I[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*f[1:(Mnage-1)]*(S[i-1,1:(Mnage-1)] + g*R[i-1,1:(Mnage-1)])*dt + (v[1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x)*f[1:(Mnage-1)]*L[i-1,1:(Mnage-1)]*dt + r[1:(Mnage-1)]*h[1:(Mnage-1)]*R[i-1,1:(Mnage-1)]*dt + w*NI[i-1,1:(Mnage-1)]*dt
      new_I_noconv[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*f[1:(Mnage-1)]*(S[i-1,1:(Mnage-1)] + g*R[i-1,1:(Mnage-1)])*dt + (v[1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x)*f[1:(Mnage-1)]*L[i-1,1:(Mnage-1)]*dt + r[1:(Mnage-1)]*h[1:(Mnage-1)]*R[i-1,1:(Mnage-1)]*dt
      new_NI[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*(1 - f[1:(Mnage-1)])*(S[i-1,1:(Mnage-1)] + g*R[i-1,1:(Mnage-1)])*dt + (v[1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x)*(1 - f[1:(Mnage-1)])*L[i-1,1:(Mnage-1)]*dt + r[1:(Mnage-1)]*(1 - h[1:(Mnage-1)])*R[i-1,1:(Mnage-1)]*dt  
      new_actv[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*S[i-1,1:(Mnage-1)]*dt + (v[1:(Mnage-1)] +lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*x)*L[i-1,1:(Mnage-1)]*dt + (r[1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*g)*R[i-1,1:(Mnage-1)]*dt
      new_notif[i,2:Mnage] = CDR[2:(Mnage)]*(new_I[i,2:Mnage] + e*new_NI[i,2:Mnage])
      new_actvinf[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*S[i-1,1:(Mnage-1)]*dt
      
      R[i,2:Mnage] = R[i-1,1:(Mnage-1)] + n[1:(Mnage-1)]*(I[i-1,1:(Mnage-1)] + NI[i-1,1:(Mnage-1)])*dt + CDR[2:(Mnage)]*CoT*(new_I[i,2:Mnage] + e*new_NI[i,2:Mnage]) - (r[1:(Mnage-1)] + g*lambda[i-1,1:(Mnage-1)] + u[1:(Mnage-1)])*R[i-1,1:(Mnage-1)]*dt 
      I[i,2:Mnage] = I[i-1,1:(Mnage-1)] + (1 - CDR[2:(Mnage)]*CoT)*(new_I[i,2:Mnage]) - (n[1:(Mnage-1)] + u[1:(Mnage-1)] + ui[1:(Mnage-1)])*I[i-1,1:(Mnage-1)]*dt
      NI[i,2:Mnage] = NI[i-1,1:(Mnage-1)] + ((1 - CDR[2:(Mnage)]*CoT*e)*new_NI[i,2:Mnage]) - (n[1:(Mnage-1)] + u[1:(Mnage-1)] + uni[1:(Mnage-1)] + w)*NI[i-1,1:(Mnage-1)]*dt                    
      
      
      #### VACCINE STRATA AGING + INFECTION ####
      ## Others have aged so need these to be reset to zero
      Sv[i,1] = 0;Lv[i,1] = 0;Rv[i,1] = 0      
      
      Sv[i,2:Mnage] = Sv[i-1,1:(Mnage-1)] - (u[1:(Mnage-1)]+(lambda[i-1,1:(Mnage-1)]))*Sv[i-1,1:(Mnage-1)]*dt 
      Lv[i,2:Mnage] = Lv[i-1,1:(Mnage-1)] + lambda[i-1,1:(Mnage-1)]*Sv[i-1,1:(Mnage-1)]*dt-u[1:(Mnage-1)]*Lv[i-1,1:(Mnage-1)]*dt+lambda[i-1,1:(Mnage-1)]*x*Rv[i-1,1:(Mnage-1)]*dt
      Rv[i,2:Mnage] = Rv[i-1,1:(Mnage-1)] -lambda[i-1,1:(Mnage-1)]*x*Rv[i-1,1:(Mnage-1)]*dt- u[1:(Mnage-1)]*Rv[i-1,1:(Mnage-1)]*dt
    
       
      #### Vaccination campaign: implement both vaccination and waning ####
      
      
      S2 = S[i,] + Sv[i,]*(d[i,]*(1-thetaS[i,])) - thetaS[i,]*S[i,]  
      L2 = L[i,] + Lv[i,]*(d[i,]*(1-thetaL[i,])) - thetaL[i,]*L[i,]
      R2 = R[i,] + Rv[i,]*(d[i,]*(1-thetaR[i,])) - thetaR[i,]*R[i,]
      
      Sv2 = Sv[i,] - Sv[i,]*(d[i,]*(1-thetaS[i,])) + thetaS[i,]*S[i,]
      Lv2 = Lv[i,] - Lv[i,]*(d[i,]*(1-thetaL[i,])) + thetaL[i,]*L[i,]
      Rv2 = Rv[i,] - Rv[i,]*(d[i,]*(1-thetaR[i,])) + thetaR[i,]*R[i,]

      
      
      ##calculate the number just vaccinated before the S gets replaced by S2 etc. Everyone except active cases included.
      if (vaccine==1){num_vac[i,] = (S[i,] + L[i,] + R[i,] + Sv[i,]*d[i,] + Lv[i,]*d[i,] + Rv[i,]*d[i,])*thetaV2[i,]/effD}
      if (vaccine==2){num_vac[i,] = (L[i,] + R[i,] + Lv[i,]*d[i,] + Rv[i,]*d[i,])*thetaV2[i,]/effD}
      if (vaccine==3){num_vac[i,] = (S[i,] + L[i,] + R[i,] + Sv[i,] + Lv[i,] + Rv[i,])*thetaV2[i,]/effD}
      
      #replace S etc with the new number post-vaccination campaign and waning
      S[i,]<-S2;L[i,]<-L2;R[i,]<-R2;   
      Sv[i,]<-Sv2;Lv[i,]<-Lv2;Rv[i,]<-Rv2;
      
      ## POPULATION SIZE ##
      psize[i]<-sum(S[i,],L[i,],R[i,],I[i,],NI[i,],Sv[i,],Lv[i,],Rv[i,])
      
      #ages needed to fit to incidence and population size
      psize014[i]<-sum(S[i,1:15],L[i,1:15],R[i,1:15],I[i,1:15],NI[i,1:15],Sv[i,1:15],Lv[i,1:15],Rv[i,1:15])
      psize1554[i]<-sum(S[i,16:55],L[i,16:55],R[i,16:55],I[i,16:55],NI[i,16:55],Sv[i,16:55],Lv[i,16:55],Rv[i,16:55])
      psize5564[i]<-sum(S[i,56:65],L[i,56:65],R[i,56:65],I[i,56:65],NI[i,56:65],Sv[i,56:65],Lv[i,56:65],Rv[i,56:65])
      psize65plus[i]<-sum(S[i,66:Mnage],L[i,66:Mnage],R[i,66:Mnage],I[i,66:Mnage],NI[i,66:Mnage],Sv[i,66:Mnage],Lv[i,66:Mnage],Rv[i,66:Mnage])
      psize55plus[i]<-sum(S[i,56:Mnage],L[i,56:Mnage],R[i,56:Mnage],I[i,56:Mnage],NI[i,56:Mnage],Sv[i,56:Mnage],Lv[i,56:Mnage],Rv[i,56:Mnage])
      psize55minus[i]<-sum(S[i,1:55],L[i,1:55],R[i,1:55],I[i,1:55],NI[i,1:55],Sv[i,1:55],Lv[i,1:55],Rv[i,1:55])
      
      #ages needed to fit to mort  and prevalence as have different groupings
      psize1559[i]<-sum(S[i,16:60],L[i,16:60],R[i,16:60],I[i,16:60],NI[i,16:60],Sv[i,16:60],Lv[i,16:60],Rv[i,16:60])
      psize1529[i]<-sum(S[i,16:30],L[i,16:30],R[i,16:30],I[i,16:30],NI[i,16:30],Sv[i,16:30],Lv[i,16:30],Rv[i,16:30])
      psize3044[i]<-sum(S[i,31:45],L[i,31:45],R[i,31:45],I[i,31:45],NI[i,31:45],Sv[i,31:45],Lv[i,31:45],Rv[i,31:45])
      psize4559[i]<-sum(S[i,46:60],L[i,46:60],R[i,46:60],I[i,46:60],NI[i,46:60],Sv[i,46:60],Lv[i,46:60],Rv[i,46:60])
      psize60plus[i]<-sum(S[i,61:Mnage],L[i,61:Mnage],R[i,61:Mnage],I[i,61:Mnage],NI[i,61:Mnage],Sv[i,61:Mnage],Lv[i,61:Mnage],Rv[i,61:Mnage])
      
      #other ages
      psize0509[i]<-sum(S[i,6:10],L[i,6:10],R[i,6:10],I[i,6:10],NI[i,6:10],Sv[i,6:10],Lv[i,6:10],Rv[i,6:10])
      psize1019[i]<-sum(S[i,11:20],L[i,11:20],R[i,11:20],I[i,11:20],NI[i,11:20],Sv[i,11:20],Lv[i,11:20],Rv[i,11:20])
      psize2029[i]<-sum(S[i,21:30],L[i,21:30],R[i,21:30],I[i,21:30],NI[i,21:30],Sv[i,21:30],Lv[i,21:30],Rv[i,21:30])
      psize3039[i]<-sum(S[i,31:40],L[i,31:40],R[i,31:40],I[i,31:40],NI[i,31:40],Sv[i,31:40],Lv[i,31:40],Rv[i,31:40])
      psize4049[i]<-sum(S[i,41:50],L[i,41:50],R[i,41:50],I[i,41:50],NI[i,41:50],Sv[i,41:50],Lv[i,41:50],Rv[i,41:50])
      psize5059[i]<-sum(S[i,51:60],L[i,51:60],R[i,51:60],I[i,51:60],NI[i,51:60],Sv[i,51:60],Lv[i,51:60],Rv[i,51:60])
      psize6069[i]<-sum(S[i,61:70],L[i,61:70],R[i,61:70],I[i,61:70],NI[i,61:70],Sv[i,61:70],Lv[i,61:70],Rv[i,61:70])
      psize70plus[i]<-sum(S[i,71:Mnage],L[i,71:Mnage],R[i,71:Mnage],I[i,71:Mnage],NI[i,71:Mnage],Sv[i,71:Mnage],Lv[i,71:Mnage],Rv[i,71:Mnage])
      psize5574[i]<-sum(S[i,56:75],L[i,56:75],R[i,56:75],I[i,56:75],NI[i,56:75],Sv[i,56:75],Lv[i,56:75],Rv[i,56:75])
      psize75plus[i]<-sum(S[i,76:Mnage],L[i,76:Mnage],R[i,76:Mnage],I[i,76:Mnage],NI[i,76:Mnage],Sv[i,76:Mnage],Lv[i,76:Mnage],Rv[i,76:Mnage])
      psize1524[i]<-sum(S[i,16:25],L[i,16:25],R[i,16:25],I[i,16:25],NI[i,16:25],Sv[i,76:Mnage],Lv[i,16:25],Rv[i,16:25])
      psize2554[i]<-sum(S[i,26:55],L[i,26:55],R[i,26:55],I[i,26:55],NI[i,26:55],Sv[i,26:55],Lv[i,26:55],Rv[i,26:55])
      psize15plus[i]<-sum(S[i,16:Mnage],L[i,16:Mnage],R[i,16:Mnage],I[i,16:Mnage],NI[i,16:Mnage],Sv[i,16:Mnage],Lv[i,16:Mnage],Rv[i,16:Mnage])
      psize60[i]<-sum(S[i,61],L[i,61],R[i,61],I[i,61],NI[i,61],Sv[i,61],Lv[i,61],Rv[i,61])
      psize15[i]<-sum(S[i,16],L[i,16],R[i,16],I[i,16],NI[i,16],Sv[i,16],Lv[i,16],Rv[i,16])
      psize6174[i]<-sum(S[i,62:75],L[i,62:75],R[i,62:75],I[i,62:75],NI[i,62:75],Lv[i,62:75],Rv[i,62:75],Iv[i,62:75],NIv[i,62:75])
      
      #Population size and incidence for contact matrices
      psizematrix[i,1]<-sum(S[i,1:5],L[i,1:5],R[i,1:5],I[i,1:5],NI[i,1:5],Sv[i,1:5],Lv[i,1:5],Rv[i,1:5])
      psizematrix[i,2]<-sum(S[i,6:10],L[i,6:10],R[i,6:10],I[i,6:10],NI[i,6:10],Sv[i,6:10],Lv[i,6:10],Rv[i,6:10])
      psizematrix[i,3]<-sum(S[i,11:15],L[i,11:15],R[i,11:15],I[i,11:15],NI[i,11:15],Sv[i,11:15],Lv[i,11:15],Rv[i,11:15])
      psizematrix[i,4]<-sum(S[i,16:20],L[i,16:20],R[i,16:20],I[i,16:20],NI[i,16:20],Sv[i,16:20],Lv[i,16:20],Rv[i,16:20])
      psizematrix[i,5]<-sum(S[i,21:25],L[i,21:25],R[i,21:25],I[i,21:25],NI[i,21:25],Sv[i,21:25],Lv[i,21:25],Rv[i,21:25])
      psizematrix[i,6]<-sum(S[i,26:30],L[i,26:30],R[i,26:30],I[i,26:30],NI[i,26:30],Sv[i,26:30],Lv[i,26:30],Rv[i,26:30])
      psizematrix[i,7]<-sum(S[i,31:35],L[i,21:65],R[i,21:65],I[i,21:65],NI[i,21:65],Sv[i,21:65],Lv[i,21:65],Rv[i,21:65])
      psizematrix[i,8]<-sum(S[i,36:40],L[i,36:40],R[i,36:40],I[i,36:40],NI[i,36:40],Sv[i,36:40],Lv[i,36:40],Rv[i,36:40])
      psizematrix[i,9]<-sum(S[i,41:45],L[i,41:45],R[i,41:45],I[i,41:45],NI[i,41:45],Sv[i,41:45],Lv[i,41:45],Rv[i,41:45])
      psizematrix[i,10]<-sum(S[i,46:50],L[i,46:50],R[i,46:50],I[i,46:50],NI[i,46:50],Sv[i,46:50],Lv[i,46:50],Rv[i,46:50])
      psizematrix[i,11]<-sum(S[i,51:55],L[i,51:55],R[i,51:55],I[i,51:55],NI[i,51:55],Sv[i,51:55],Lv[i,51:55],Rv[i,51:55])
      psizematrix[i,12]<-sum(S[i,56:60],L[i,56:60],R[i,56:60],I[i,56:60],NI[i,56:60],Sv[i,56:60],Lv[i,56:60],Rv[i,56:60])
      psizematrix[i,13]<-sum(S[i,61:65],L[i,61:65],R[i,61:65],I[i,61:65],NI[i,61:65],Sv[i,61:65],Lv[i,61:65],Rv[i,61:65])
      psizematrix[i,14]<-sum(S[i,66:70],L[i,66:70],R[i,66:70],I[i,66:70],NI[i,66:70],Sv[i,66:70],Lv[i,66:70],Rv[i,66:70])
      psizematrix[i,15]<-sum(S[i,71:75],L[i,71:75],R[i,71:75],I[i,71:75],NI[i,71:75],Sv[i,71:75],Lv[i,71:75],Rv[i,71:75])
      psizematrix[i,16]<-sum(S[i,76:Mnage],L[i,76:Mnage],R[i,76:Mnage],I[i,76:Mnage],NI[i,76:Mnage],Sv[i,76:Mnage],Lv[i,76:Mnage],Rv[i,76:Mnage])
      
      
      Imatrix[i,1]<-sum(I[i,1:5])
      Imatrix[i,2]<-sum(I[i,6:10])
      Imatrix[i,3]<-sum(I[i,11:15])
      Imatrix[i,4]<-sum(I[i,16:20])
      Imatrix[i,5]<-sum(I[i,21:25])
      Imatrix[i,6]<-sum(I[i,26:30])
      Imatrix[i,7]<-sum(I[i,31:35])
      Imatrix[i,8]<-sum(I[i,36:40])
      Imatrix[i,9]<-sum(I[i,41:45])
      Imatrix[i,10]<-sum(I[i,46:50])
      Imatrix[i,11]<-sum(I[i,51:55])
      Imatrix[i,12]<-sum(I[i,56:60])
      Imatrix[i,13]<-sum(I[i,61:65])
      Imatrix[i,14]<-sum(I[i,66:70])
      Imatrix[i,15]<-sum(I[i,71:75])
      Imatrix[i,16]<-sum(I[i,76:Mnage])
      
      # TB deaths
      TBDeaths[i,1:Mnage]=dt*((ui[1:Mnage]*I[i-1,1:Mnage])+(uni[1:Mnage]*NI[i-1,1:Mnage]));
      # All deaths
      ADeaths[i,]=dt*(u*S[i-1,]+u*L[i-1,]+(u+ui)*I[i-1,]+(u+uni)*NI[i-1,]+u*R[i-1,]+u*Sv[i-1,]+u*Lv[i-1,]+u*Rv[i-1,])
      
    } #### END OF START OF YEAR RUNS
    
    #### MIDDLE OF YEAR RUNS
    
    for (i in (2+(1/dt)*(k-year1)):((1/dt)*(k-year1)+1/dt)){
      start <- 0 # Not the start of the year
      
      lambda[i-1,1:Mnage] <- t(neta * (1-exp(colSums(-(myneta[1:16,1:Mnage]) * z * ((Imatrix[i-1,1:16])/(psizematrix[i-1,1:16]))))))
      
      ##  Scaling of population size in 1950
      if (k == 1950) {demscale<-psize[i-1,1]/psz1900
      S[i-1,]<-S[i-1,]/demscale 
      L[i-1,]<-L[i-1,]/demscale 
      I[i-1,]<-I[i-1,]/demscale 
      NI[i-1,]<-NI[i-1,]/demscale 
      R[i-1,]<-R[i-1,]/demscale 
      }
      
      #### TB model ####
      
      S[i,1:Mnage] = S[i-1,1:Mnage] - (u[1:Mnage]+lambda[i-1,1:Mnage])*S[i-1,1:Mnage]*dt
      L[i,1:Mnage] = L[i-1,1:Mnage] + lambda[i-1,1:Mnage]*(1 - p[1:Mnage])*(S[i-1,1:Mnage] + g*R[i-1,1:Mnage])*dt - (v[1:Mnage] + lambda[i-1,1:Mnage]*p[1:Mnage]*x + u[1:Mnage])*L[i-1,1:Mnage]*dt 
      
      new_infect[i,1:Mnage] = lambda[i-1,1:(Mnage)]*S[i-1,1:(Mnage)]*dt + lambda[i-1,1:(Mnage)]*(x*L[i-1,1:(Mnage)] + g*R[i-1,1:(Mnage)])*dt  
      new_I_react[i,1:Mnage] = v[1:Mnage]*f[1:(Mnage)]*(L[i-1,1:(Mnage)])*dt + r[1:Mnage]*h[1:(Mnage)]*R[i-1,1:(Mnage)]*dt 
      new_NI_react[i,1:Mnage] =  v[1:Mnage]*(1 - f[1:(Mnage)])*L[i-1,1:(Mnage)]*dt + r[1:Mnage]*(1 - h[1:(Mnage)])*R[i-1,1:(Mnage)]*dt  
      new_actv_react[i,1:Mnage] = v[1:Mnage]*(L[i-1,1:Mnage])*dt + r[1:Mnage]*R[i-1,1:Mnage]*dt 
      new_actv_inf[i,1:Mnage] = lambda[i-1,1:Mnage]*p[1:Mnage]*S[i-1,1:Mnage]*dt + lambda[i-1,1:Mnage]*p[1:Mnage]*x*(L[i-1,1:Mnage])*dt + lambda[i-1,1:Mnage]*p[1:Mnage]*g*R[i-1,1:Mnage]*dt 
      new_I[i,1:Mnage] = lambda[i-1,1:Mnage]*p[1:Mnage]*f[1:Mnage]*(S[i-1,1:Mnage] + g*R[i-1,1:(Mnage)])*dt + (v[1:Mnage] + lambda[i-1,1:Mnage]*p[1:(Mnage)]*x)*f[1:(Mnage)]*L[i-1,1:(Mnage)]*dt + r[1:Mnage]*h[1:(Mnage)]*R[i-1,1:(Mnage)]*dt + w*NI[i-1,1:(Mnage)]*dt
      new_I_noconv[i,1:Mnage] = lambda[i-1,1:Mnage]*p[1:Mnage]*f[1:Mnage]*(S[i-1,1:Mnage] + g*R[i-1,1:Mnage])*dt + (v[1:Mnage] + lambda[i-1,1:Mnage]*p[1:Mnage]*x)*f[1:Mnage]*L[i-1,1:Mnage]*dt + r[1:Mnage]*h[1:Mnage]*R[i-1,1:Mnage]*dt
      new_NI[i,1:Mnage] = lambda[i-1,1:Mnage]*p[1:(Mnage)]*(1 - f[1:(Mnage)])*(S[i-1,1:(Mnage)] + g*R[i-1,1:(Mnage)])*dt + (v[1:Mnage] + lambda[i-1,1:Mnage]*p[1:(Mnage)]*x)*(1 - f[1:(Mnage)])*L[i-1,1:(Mnage)]*dt + r[1:Mnage]*(1 - h[1:(Mnage)])*R[i-1,1:(Mnage)]*dt  
      new_actv[i,1:Mnage] = lambda[i-1,1:Mnage]*p[1:Mnage]*S[i-1,1:Mnage]*dt + (v[1:Mnage] + lambda[i-1,1:Mnage]*p[1:Mnage]*x)*L[i-1,1:Mnage]*dt + (r[1:Mnage] + lambda[i-1,1:Mnage]*p[1:Mnage]*g)*R[i-1,1:Mnage]*dt
      new_notif[i,1:Mnage] = CDR[1:(Mnage)]*(new_I[i,1:Mnage] + e*new_NI[i,1:Mnage])
      new_actvinf[i,2:Mnage] = lambda[i-1,1:(Mnage-1)]*p[1:(Mnage-1)]*S[i-1,1:(Mnage-1)]*dt
      
      R[i,1:Mnage] = R[i-1,1:Mnage] + n[1:Mnage]*(I[i-1,1:(Mnage)] + NI[i-1,1:(Mnage)])*dt + CDR[1:Mnage]*CoT*(new_I[i,1:Mnage] + e*new_NI[i,1:Mnage]) - (r[1:Mnage] + g*lambda[i-1,1:Mnage] + u[1:(Mnage)])*R[i-1,1:(Mnage)]*dt 
      I[i,1:Mnage] = I[i-1,1:Mnage] + (1 - CDR[1:Mnage]*CoT)*(new_I[i,1:Mnage]) - (n[1:Mnage] + u[1:(Mnage)] + ui[1:Mnage])*I[i-1,1:Mnage]*dt
      NI[i,1:Mnage] = NI[i-1,1:Mnage] + ((1 - CDR[1:Mnage]*CoT*e)*new_NI[i,1:Mnage]) - (n[1:Mnage] + u[1:Mnage] + uni[1:Mnage] + w)*NI[i-1,1:(Mnage)]*dt                    
      
      #############
      
      Sv[i,1:Mnage] = Sv[i-1,1:Mnage] - (u[1:Mnage]+(lambda[i-1,1:Mnage]))*Sv[i-1,1:Mnage]*dt 
      Lv[i,1:Mnage] = Lv[i-1,1:Mnage] + lambda[i-1,1:Mnage]*Sv[i-1,1:Mnage]*dt-u[1:Mnage]*Lv[i-1,1:Mnage]*dt+lambda[i-1,1:Mnage]*x*Rv[i-1,1:Mnage]*dt
      Rv[i,1:Mnage] = Rv[i-1,1:Mnage] -lambda[i-1,1:Mnage]*x*Rv[i-1,1:Mnage]*dt- u[1:Mnage]*Rv[i-1,1:Mnage]*dt
      
      #### Vaccine coverage and duration ####
      
      ## Vaccination and Waning ##
      S2 = S[i,] + Sv[i,]*(d[i,]*(1-thetaS[i,])) - thetaS[i,]*S[i,]  
      L2 = L[i,] + Lv[i,]*(d[i,]*(1-thetaL[i,])) - thetaL[i,]*L[i,]
      R2 = R[i,] + Rv[i,]*(d[i,]*(1-thetaR[i,])) - thetaR[i,]*R[i,]
      
      Sv2 = Sv[i,] - Sv[i,]*(d[i,]*(1-thetaS[i,])) + thetaS[i,]*S[i,]
      Lv2 = Lv[i,] - Lv[i,]*(d[i,]*(1-thetaL[i,])) + thetaL[i,]*L[i,]
      Rv2 = Rv[i,] - Rv[i,]*(d[i,]*(1-thetaR[i,])) + thetaR[i,]*R[i,]
      
      ##calculate the number just vaccinated before the S gets replaced by S2 etc. Everyone except active cases included.
      if (vaccine==1){num_vac[i,] = (S[i,] + L[i,] + R[i,] + Sv[i,]*d[i,] + Lv[i,]*d[i,] + Rv[i,]*d[i,])*thetaV2[i,]/effD}
      if (vaccine==2){num_vac[i,] = (L[i,] + R[i,] + Lv[i,]*d[i,] + Rv[i,]*d[i,])*thetaV2[i,]/effD}
      if (vaccine==3){num_vac[i,] = (S[i,] + L[i,] + R[i,] + Sv[i,] + Lv[i,] + Rv[i,])*thetaV2[i,]/effD}
      
      S[i,]<-S2;L[i,]<-L2;R[i,]<-R2;   
      Sv[i,]<-Sv2;Lv[i,]<-Lv2;Rv[i,]<-Rv2;
      
      
      ## POPULATION SIZE ##
      psize[i]<-sum(S[i,],L[i,],R[i,],I[i,],NI[i,],Sv[i,],Lv[i,],Rv[i,])
      
      #ages needed to fit to incidence and population size
      psize014[i]<-sum(S[i,1:15],L[i,1:15],R[i,1:15],I[i,1:15],NI[i,1:15],Sv[i,1:15],Lv[i,1:15],Rv[i,1:15])
      psize1554[i]<-sum(S[i,16:55],L[i,16:55],R[i,16:55],I[i,16:55],NI[i,16:55],Sv[i,16:55],Lv[i,16:55],Rv[i,16:55])
      psize5564[i]<-sum(S[i,56:65],L[i,56:65],R[i,56:65],I[i,56:65],NI[i,56:65],Sv[i,56:65],Lv[i,56:65],Rv[i,56:65])
      psize65plus[i]<-sum(S[i,66:Mnage],L[i,66:Mnage],R[i,66:Mnage],I[i,66:Mnage],NI[i,66:Mnage],Sv[i,66:Mnage],Lv[i,66:Mnage],Rv[i,66:Mnage])
      psize55plus[i]<-sum(S[i,56:Mnage],L[i,56:Mnage],R[i,56:Mnage],I[i,56:Mnage],NI[i,56:Mnage],Sv[i,56:Mnage],Lv[i,56:Mnage],Rv[i,56:Mnage])
      psize55minus[i]<-sum(S[i,1:55],L[i,1:55],R[i,1:55],I[i,1:55],NI[i,1:55],Sv[i,1:55],Lv[i,1:55],Rv[i,1:55])
      
      #ages needed to fit to mort  and prevalence as have different groupings
      psize1559[i]<-sum(S[i,16:60],L[i,16:60],R[i,16:60],I[i,16:60],NI[i,16:60],Sv[i,16:60],Lv[i,16:60],Rv[i,16:60])
      psize1529[i]<-sum(S[i,16:30],L[i,16:30],R[i,16:30],I[i,16:30],NI[i,16:30],Sv[i,16:30],Lv[i,16:30],Rv[i,16:30])
      psize3044[i]<-sum(S[i,31:45],L[i,31:45],R[i,31:45],I[i,31:45],NI[i,31:45],Sv[i,31:45],Lv[i,31:45],Rv[i,31:45])
      psize4559[i]<-sum(S[i,46:60],L[i,46:60],R[i,46:60],I[i,46:60],NI[i,46:60],Sv[i,46:60],Lv[i,46:60],Rv[i,46:60])
      psize60plus[i]<-sum(S[i,61:Mnage],L[i,61:Mnage],R[i,61:Mnage],I[i,61:Mnage],NI[i,61:Mnage],Sv[i,61:Mnage],Lv[i,61:Mnage],Rv[i,61:Mnage])
      
      #other ages
      psize0509[i]<-sum(S[i,6:10],L[i,6:10],R[i,6:10],I[i,6:10],NI[i,6:10],Sv[i,6:10],Lv[i,6:10],Rv[i,6:10])
      psize1019[i]<-sum(S[i,11:20],L[i,11:20],R[i,11:20],I[i,11:20],NI[i,11:20],Sv[i,11:20],Lv[i,11:20],Rv[i,11:20])
      psize2029[i]<-sum(S[i,21:30],L[i,21:30],R[i,21:30],I[i,21:30],NI[i,21:30],Sv[i,21:30],Lv[i,21:30],Rv[i,21:30])
      psize3039[i]<-sum(S[i,31:40],L[i,31:40],R[i,31:40],I[i,31:40],NI[i,31:40],Sv[i,31:40],Lv[i,31:40],Rv[i,31:40])
      psize4049[i]<-sum(S[i,41:50],L[i,41:50],R[i,41:50],I[i,41:50],NI[i,41:50],Sv[i,41:50],Lv[i,41:50],Rv[i,41:50])
      psize5059[i]<-sum(S[i,51:60],L[i,51:60],R[i,51:60],I[i,51:60],NI[i,51:60],Sv[i,51:60],Lv[i,51:60],Rv[i,51:60])
      psize6069[i]<-sum(S[i,61:70],L[i,61:70],R[i,61:70],I[i,61:70],NI[i,61:70],Sv[i,61:70],Lv[i,61:70],Rv[i,61:70])
      psize70plus[i]<-sum(S[i,71:Mnage],L[i,71:Mnage],R[i,71:Mnage],I[i,71:Mnage],NI[i,71:Mnage],Sv[i,71:Mnage],Lv[i,71:Mnage],Rv[i,71:Mnage])
      psize5574[i]<-sum(S[i,56:75],L[i,56:75],R[i,56:75],I[i,56:75],NI[i,56:75],Sv[i,56:75],Lv[i,56:75],Rv[i,56:75])
      psize75plus[i]<-sum(S[i,76:Mnage],L[i,76:Mnage],R[i,76:Mnage],I[i,76:Mnage],NI[i,76:Mnage],Sv[i,76:Mnage],Lv[i,76:Mnage],Rv[i,76:Mnage])
      psize1524[i]<-sum(S[i,16:25],L[i,16:25],R[i,16:25],I[i,16:25],NI[i,16:25],Sv[i,76:Mnage],Lv[i,16:25],Rv[i,16:25])
      psize2554[i]<-sum(S[i,26:55],L[i,26:55],R[i,26:55],I[i,26:55],NI[i,26:55],Sv[i,26:55],Lv[i,26:55],Rv[i,26:55])
      psize15plus[i]<-sum(S[i,16:Mnage],L[i,16:Mnage],R[i,16:Mnage],I[i,16:Mnage],NI[i,16:Mnage],Sv[i,16:Mnage],Lv[i,16:Mnage],Rv[i,16:Mnage])
      psize60[i]<-sum(S[i,61],L[i,61],R[i,61],I[i,61],NI[i,61],Sv[i,61],Lv[i,61],Rv[i,61])
      psize15[i]<-sum(S[i,16],L[i,16],R[i,16],I[i,16],NI[i,16],Sv[i,16],Lv[i,16],Rv[i,16])
      psize6174[i]<-sum(S[i,62:75],L[i,62:75],R[i,62:75],I[i,62:75],NI[i,62:75],Lv[i,62:75],Rv[i,62:75],Iv[i,62:75],NIv[i,62:75])
      
      #Population size and incidence for contact matrices
      psizematrix[i,1]<-sum(S[i,1:5],L[i,1:5],R[i,1:5],I[i,1:5],NI[i,1:5],Sv[i,1:5],Lv[i,1:5],Rv[i,1:5])
      psizematrix[i,2]<-sum(S[i,6:10],L[i,6:10],R[i,6:10],I[i,6:10],NI[i,6:10],Sv[i,6:10],Lv[i,6:10],Rv[i,6:10])
      psizematrix[i,3]<-sum(S[i,11:15],L[i,11:15],R[i,11:15],I[i,11:15],NI[i,11:15],Sv[i,11:15],Lv[i,11:15],Rv[i,11:15])
      psizematrix[i,4]<-sum(S[i,16:20],L[i,16:20],R[i,16:20],I[i,16:20],NI[i,16:20],Sv[i,16:20],Lv[i,16:20],Rv[i,16:20])
      psizematrix[i,5]<-sum(S[i,21:25],L[i,21:25],R[i,21:25],I[i,21:25],NI[i,21:25],Sv[i,21:25],Lv[i,21:25],Rv[i,21:25])
      psizematrix[i,6]<-sum(S[i,26:30],L[i,26:30],R[i,26:30],I[i,26:30],NI[i,26:30],Sv[i,26:30],Lv[i,26:30],Rv[i,26:30])
      psizematrix[i,7]<-sum(S[i,31:35],L[i,21:65],R[i,21:65],I[i,21:65],NI[i,21:65],Sv[i,21:65],Lv[i,21:65],Rv[i,21:65])
      psizematrix[i,8]<-sum(S[i,36:40],L[i,36:40],R[i,36:40],I[i,36:40],NI[i,36:40],Sv[i,36:40],Lv[i,36:40],Rv[i,36:40])
      psizematrix[i,9]<-sum(S[i,41:45],L[i,41:45],R[i,41:45],I[i,41:45],NI[i,41:45],Sv[i,41:45],Lv[i,41:45],Rv[i,41:45])
      psizematrix[i,10]<-sum(S[i,46:50],L[i,46:50],R[i,46:50],I[i,46:50],NI[i,46:50],Sv[i,46:50],Lv[i,46:50],Rv[i,46:50])
      psizematrix[i,11]<-sum(S[i,51:55],L[i,51:55],R[i,51:55],I[i,51:55],NI[i,51:55],Sv[i,51:55],Lv[i,51:55],Rv[i,51:55])
      psizematrix[i,12]<-sum(S[i,56:60],L[i,56:60],R[i,56:60],I[i,56:60],NI[i,56:60],Sv[i,56:60],Lv[i,56:60],Rv[i,56:60])
      psizematrix[i,13]<-sum(S[i,61:65],L[i,61:65],R[i,61:65],I[i,61:65],NI[i,61:65],Sv[i,61:65],Lv[i,61:65],Rv[i,61:65])
      psizematrix[i,14]<-sum(S[i,66:70],L[i,66:70],R[i,66:70],I[i,66:70],NI[i,66:70],Sv[i,66:70],Lv[i,66:70],Rv[i,66:70])
      psizematrix[i,15]<-sum(S[i,71:75],L[i,71:75],R[i,71:75],I[i,71:75],NI[i,71:75],Sv[i,71:75],Lv[i,71:75],Rv[i,71:75])
      psizematrix[i,16]<-sum(S[i,76:Mnage],L[i,76:Mnage],R[i,76:Mnage],I[i,76:Mnage],NI[i,76:Mnage],Sv[i,76:Mnage],Lv[i,76:Mnage],Rv[i,76:Mnage])
      
      
      Imatrix[i,1]<-sum(I[i,1:5])
      Imatrix[i,2]<-sum(I[i,6:10])
      Imatrix[i,3]<-sum(I[i,11:15])
      Imatrix[i,4]<-sum(I[i,16:20])
      Imatrix[i,5]<-sum(I[i,21:25])
      Imatrix[i,6]<-sum(I[i,26:30])
      Imatrix[i,7]<-sum(I[i,31:35])
      Imatrix[i,8]<-sum(I[i,36:40])
      Imatrix[i,9]<-sum(I[i,41:45])
      Imatrix[i,10]<-sum(I[i,46:50])
      Imatrix[i,11]<-sum(I[i,51:55])
      Imatrix[i,12]<-sum(I[i,56:60])
      Imatrix[i,13]<-sum(I[i,61:65])
      Imatrix[i,14]<-sum(I[i,66:70])
      Imatrix[i,15]<-sum(I[i,71:75])
      Imatrix[i,16]<-sum(I[i,76:Mnage])
      
      # TB deaths
      TBDeaths[i,1:Mnage]=dt*((ui[1:Mnage]*I[i-1,1:Mnage])+(uni[1:Mnage]*NI[i-1,1:Mnage]));
      # All deaths
      ADeaths[i,]=dt*(u*S[i-1,]+u*L[i-1,]+(u+ui)*I[i-1,]+(u+uni)*NI[i-1,]+u*R[i-1,]+u*Sv[i-1,]+u*Lv[i-1,]+u*Rv[i-1,])
      
      # If last timestep of the year
      if(i == ((1/dt)*(k-year1)+1/dt)){
        
        i1<-((1/dt)*(k-year1)+1); i2<-((1/dt)*(k-year1)+1/dt); i3<-((1/dt)*(k-year1)+2) #i1 is number of first time step of the yr, i2 is number of last timestep of the yr, i3 is second timestep of the year
        
        # TB INCIDENCE AND MORTALITY etc for MODEL FITTING and RESEARCH OUTCOMES
        
        ## (1) population size 
        PSIZEy[(k-year1+1),1]<-mean(psize[i1:i2]); PSIZEy[(k-year1+1),2]<-mean(psize014[i1:i2]); PSIZEy[(k-year1+1),3]<-mean(psize1554[i1:i2]); PSIZEy[(k-year1+1),4]<-mean(psize5564[i1:i2])
        PSIZEy[(k-year1+1),5]<-mean(psize65plus[i1:i2]); PSIZEy[(k-year1+1),6]<-mean(psize1559[i1:i2]); PSIZEy[(k-year1+1),7]<-mean(psize1529[i1:i2]); PSIZEy[(k-year1+1),8]<-mean(psize3044[i1:i2]); PSIZEy[(k-year1+1),9]<-mean(psize4559[i1:i2])
        PSIZEy[(k-year1+1),10]<-mean(psize60plus[i1:i2]); PSIZEy[(k-year1+1),11]<-mean(psize55plus[i1:i2]); PSIZEy[(k-year1+1),12]<-mean(psize0509[i1:i2])
        PSIZEy[(k-year1+1),13]<-mean(psize1019[i1:i2]); PSIZEy[(k-year1+1),14]<-mean(psize2029[i1:i2]); PSIZEy[(k-year1+1),15]<-mean(psize3039[i1:i2])
        PSIZEy[(k-year1+1),16]<-mean(psize4049[i1:i2]); PSIZEy[(k-year1+1),17]<-mean(psize5059[i1:i2]); PSIZEy[(k-year1+1),18]<-mean(psize6069[i1:i2])
        PSIZEy[(k-year1+1),19]<-mean(psize70plus[i1:i2]); PSIZEy[(k-year1+1),20]<-mean(psize5574[i1:i2]); PSIZEy[(k-year1+1),21]<-mean(psize75plus[i1:i2])
        PSIZEy[(k-year1+1),22]<-mean(psize1524[i1:i2]); PSIZEy[(k-year1+1),23]<-mean(psize2554[i1:i2]); PSIZEy[(k-year1+1),24]<-mean(psize15plus[i1:i2])
        PSIZEy[(k-year1+1),25]<-mean(psize15[i1:i2]); PSIZEy[(k-year1+1),26]<-mean(psize60[i1:i2]); PSIZEy[(k-year1+1),27]<-mean(psize6174[i1:i2])
        
        ## (2) TB incidence rate
        TBI[(k-year1+1),1]<-100000*sum(new_I_noconv[i1:i2,],new_NI[i1:i2,],new_Iv_noconv[i1:i2,],new_NIv[i1:i2,])/mean(psize[i1:i2])
        TBI[(k-year1+1),2]<-100000*sum(new_I_noconv[i1:i2,1:15],new_NI[i1:i2,1:15],new_Iv_noconv[i1:i2,1:15],new_NIv[i1:i2,1:15])/mean(psize014[i1:i2])
        TBI[(k-year1+1),3]<-100000*sum(new_I_noconv[i1:i2,16:55],new_NI[i1:i2,16:55],new_Iv_noconv[i1:i2,16:55],new_NIv[i1:i2,16:55])/mean(psize1554[i1:i2])
        TBI[(k-year1+1),4]<-100000*sum(new_I_noconv[i1:i2,56:65],new_NI[i1:i2,56:65],new_Iv_noconv[i1:i2,56:65],new_NIv[i1:i2,56:65])/mean(psize5564[i1:i2])
        TBI[(k-year1+1),5]<-100000*sum(new_I_noconv[i1:i2,66:Mnage],new_NI[i1:i2,66:Mnage],new_Iv_noconv[i1:i2,66:Mnage],new_NIv[i1:i2,66:Mnage])/mean(psize65plus[i1:i2])
        TBI[(k-year1+1),6]<-100000*sum(new_I_noconv[i1:i2,56:Mnage],new_NI[i1:i2,56:Mnage],new_Iv_noconv[i1:i2,56:Mnage],new_NIv[i1:i2,56:Mnage])/mean(psize55plus[i1:i2])
        TBI[(k-year1+1),7]<-100000*sum(new_I_noconv[i1:i2,1:55],new_NI[i1:i2,1:55],new_Iv_noconv[i1:i2,1:55],new_NIv[i1:i2,1:55])/mean(psize55minus[i1:i2])
        TBI[(k-year1+1),8]<-100000*sum(new_I_noconv[i1:i2,16:25],new_NI[i1:i2,16:25],new_Iv_noconv[i1:i2,16:25],new_NIv[i1:i2,16:25])/mean(psize1524[i1:i2])
        TBI[(k-year1+1),9]<-100000*sum(new_I_noconv[i1:i2,26:55],new_NI[i1:i2,26:55],new_Iv_noconv[i1:i2,26:55],new_NIv[i1:i2,26:55])/mean(psize2554[i1:i2])
        TBI[(k-year1+1),10]<-100000*sum(new_I_noconv[i1:i2,16:Mnage],new_NI[i1:i2,16:Mnage],new_Iv_noconv[i1:i2,16:Mnage],new_NIv[i1:i2,16:Mnage])/mean(psize15plus[i1:i2])
        
        
        ## (2b) TB treatment initiation rate
        TBN[(k-year1+1),1]<-100000*sum(new_notif[i1:i2,],new_notifv[i1:i2,])/mean(psize[i1:i2])
        TBN[(k-year1+1),2]<-100000*sum(new_notif[i1:i2,1:15],new_notifv[i1:i2,1:15])/mean(psize014[i1:i2])
        TBN[(k-year1+1),3]<-100000*sum(new_notif[i1:i2,16:55],new_notifv[i1:i2,16:55])/mean(psize1554[i1:i2])
        TBN[(k-year1+1),4]<-100000*sum(new_notif[i1:i2,56:65],new_notifv[i1:i2,56:65])/mean(psize5564[i1:i2])
        TBN[(k-year1+1),5]<-100000*sum(new_notif[i1:i2,66:Mnage],new_notifv[i1:i2,66:Mnage])/mean(psize65plus[i1:i2])  
        TBN[(k-year1+1),6]<-100000*sum(new_notif[i1:i2,56:Mnage],new_notifv[i1:i2,56:Mnage])/mean(psize55plus[i1:i2])  
        TBN[(k-year1+1),7]<-100000*sum(new_notif[i1:i2,1:55],new_notifv[i1:i2,1:55])/mean(psize55minus[i1:i2])
        TBN[(k-year1+1),8]<-100000*sum(new_notif[i1:i2,16:Mnage],new_notifv[i1:i2,16:Mnage])/mean(psize15plus[i1:i2])  
        TBN[(k-year1+1),9]<-100000*sum(new_notif[i1:i2,16:65],new_notifv[i1:i2,16:65])/mean(psize1554[i1:i2]+psize5564[i1:i2])
        
        ## (3) TB prevalence rate 
        
        TBP[(k-year1+1),1]<-100000*(sum(I[i1:i2,],NI[i1:i2,],Iv[i1:i2,],NIv[i1:i2,])/(1/dt))/mean(psize[i1:i2])
        TBP[(k-year1+1),2]<-100000*(sum(I[i1:i2,1:15],NI[i1:i2,1:15],Iv[i1:i2,1:15],NIv[i1:i2,1:15])/(1/dt))/mean(psize014[i1:i2])
        TBP[(k-year1+1),3]<-100000*(sum(I[i1:i2,16:30],NI[i1:i2,16:30],Iv[i1:i2,16:30],NIv[i1:i2,16:30])/(1/dt))/mean(psize1529[i1:i2])
        TBP[(k-year1+1),4]<-100000*(sum(I[i1:i2,31:45],NI[i1:i2,31:45],Iv[i1:i2,31:45],NIv[i1:i2,31:45])/(1/dt))/mean(psize3044[i1:i2])
        TBP[(k-year1+1),5]<-100000*(sum(I[i1:i2,46:60],NI[i1:i2,46:60],Iv[i1:i2,46:60],NIv[i1:i2,46:60])/(1/dt))/mean(psize4559[i1:i2])
        TBP[(k-year1+1),6]<-100000*(sum(I[i1:i2,61:Mnage],NI[i1:i2,61:Mnage],Iv[i1:i2,61:Mnage],NIv[i1:i2,61:Mnage])/(1/dt))/mean(psize60plus[i1:i2])
        TBP[(k-year1+1),7]<-100000*(sum(I[i1:i2,56:Mnage],NI[i1:i2,56:Mnage],Iv[i1:i2,56:Mnage],NIv[i1:i2,56:Mnage])/(1/dt))/mean(psize55plus[i1:i2])
        TBP[(k-year1+1),8]<-100000*(sum(I[i1:i2,31:60],NI[i1:i2,31:60],Iv[i1:i2,31:60],NIv[i1:i2,31:60])/(1/dt))/mean(psize3044[i1:i2]+psize4559[i1:i2])
        
        ## (4) TB bacteriologically positive prevalence rate - comparable to calibration data
        
        TBPb[(k-year1+1),1]<-100000*((sum(I[i1:i2,],Iv[i1:i2,]))/(1/dt))/mean(psize[i1:i2])
        TBPb[(k-year1+1),2]<-100000*((sum(I[i1:i2,1:15],Iv[i1:i2,1:15]))/(1/dt))/mean(psize014[i1:i2])
        TBPb[(k-year1+1),3]<-100000*((sum(I[i1:i2,16:30],Iv[i1:i2,16:30]))/(1/dt))/mean(psize1529[i1:i2])
        TBPb[(k-year1+1),4]<-100000*((sum(I[i1:i2,31:45],Iv[i1:i2,31:45]))/(1/dt))/mean(psize3044[i1:i2])
        TBPb[(k-year1+1),5]<-100000*((sum(I[i1:i2,46:60],Iv[i1:i2,46:60]))/(1/dt))/mean(psize4559[i1:i2])
        TBPb[(k-year1+1),6]<-100000*((sum(I[i1:i2,61:Mnage],Iv[i1:i2,61:Mnage]))/(1/dt))/mean(psize60plus[i1:i2])
        TBPb[(k-year1+1),7]<-100000*((sum(I[i1:i2,56:Mnage],Iv[i1:i2,56:Mnage]))/(1/dt))/mean(psize55plus[i1:i2])
        TBPb[(k-year1+1),8]<-100000*((sum(I[i1:i2,16:Mnage],Iv[i1:i2,16:Mnage]))/(1/dt))/mean(psize15plus[i1:i2])
        TBPb[(k-year1+1),9]<-100000*((sum(I[i1:i2,16:45],Iv[i1:i2,16:45]))/(1/dt))/((mean(psize1529[i1:i2]))+(mean(psize3044[i1:i2])))
        
        ## (5) TB mortality rate
        TBM[(k-year1+1),1]<-100000*sum(TBDeaths[i1:i2,])/mean(psize[i1:i2])
        TBM[(k-year1+1),2]<-100000*sum(TBDeaths[i1:i2,1:15])/mean(psize014[i1:i2])
        TBM[(k-year1+1),3]<-100000*sum(TBDeaths[i1:i2,16:55])/mean(psize1554[i1:i2])
        TBM[(k-year1+1),4]<-100000*sum(TBDeaths[i1:i2,56:65])/mean(psize5564[i1:i2])
        TBM[(k-year1+1),5]<-100000*sum(TBDeaths[i1:i2,66:Mnage])/mean(psize65plus[i1:i2])
        TBM[(k-year1+1),6]<-100000*sum(TBDeaths[i1:i2,16:60])/mean(psize1559[i1:i2])
        TBM[(k-year1+1),7]<-100000*sum(TBDeaths[i1:i2,61:Mnage])/mean(psize60plus[i1:i2])
        TBM[(k-year1+1),8]<-100000*sum(TBDeaths[i1:i2,56:Mnage])/mean(psize55plus[i1:i2])
        TBM[(k-year1+1),9]<-100000*sum(TBDeaths[i1:i2,1:55])/mean(psize55minus[i1:i2])
        TBM[(k-year1+1),10]<-100000*sum(TBDeaths[i1:i2,16:25])/mean(psize1524[i1:i2])
        TBM[(k-year1+1),11]<-100000*sum(TBDeaths[i1:i2,26:55])/mean(psize2554[i1:i2])
        TBM[(k-year1+1),12]<-100000*sum(TBDeaths[i1:i2,16:30])/mean(psize1529[i1:i2])
        TBM[(k-year1+1),13]<-100000*sum(TBDeaths[i1:i2,31:45])/mean(psize3044[i1:i2])
        TBM[(k-year1+1),14]<-100000*sum(TBDeaths[i1:i2,46:60])/mean(psize4559[i1:i2])
        
        
        ## (6) Prevalence of latent infection 
        TBPI[(k-year1+1),1]<-100*(((sum(L[i1,])/psize[i1])+(sum(L[i2,])/psize[i2]))/2)
        TBPI[(k-year1+1),2]<-100*(((sum(L[i1,1:15])/psize014[i1])+(sum(L[i2,1:15])/psize014[i2]))/2)
        TBPI[(k-year1+1),3]<-100*(((sum(L[i1,16:55])/psize1554[i1])+(sum(L[i2,16:55])/psize1554[i2]))/2)
        TBPI[(k-year1+1),4]<-100*(((sum(L[i1,56:65])/psize5564[i1])+(sum(L[i2,56:65])/psize5564[i2]))/2)
        TBPI[(k-year1+1),5]<-100*(((sum(L[i1,66:Mnage])/psize65plus[i1])+(sum(L[i2,66:Mnage])/psize65plus[i2]))/2)
        TBPI[(k-year1+1),6]<-100*(((sum(L[i1,56:Mnage])/psize55plus[i1])+(sum(L[i2,56:Mnage])/psize55plus[i2]))/2)
        TBPI[(k-year1+1),7]<-100*(((sum(L[i1,6:10])/psize0509[i1])+(sum(L[i2,6:10])/psize0509[i2]))/2)
        TBPI[(k-year1+1),8]<-100*(((sum(L[i1,11:20])/psize1019[i1])+(sum(L[i2,11:20])/psize1019[i2]))/2)
        TBPI[(k-year1+1),9]<-100*(((sum(L[i1,21:30])/psize2029[i1])+(sum(L[i2,21:30])/psize2029[i2]))/2)
        TBPI[(k-year1+1),10]<-100*(((sum(L[i1,31:40])/psize3039[i1])+(sum(L[i2,31:40])/psize3039[i2]))/2)
        TBPI[(k-year1+1),11]<-100*(((sum(L[i1,41:50])/psize4049[i1])+(sum(L[i2,41:50])/psize4049[i2]))/2)
        TBPI[(k-year1+1),12]<-100*(((sum(L[i1,51:60])/psize5059[i1])+(sum(L[i2,51:60])/psize5059[i2]))/2)
        TBPI[(k-year1+1),13]<-100*(((sum(L[i1,61:70])/psize6069[i1])+(sum(L[i2,61:70])/psize6069[i2]))/2)
        TBPI[(k-year1+1),14]<-100*(((sum(L[i1,71:Mnage])/psize70plus[i1])+(sum(L[i2,71:Mnage])/psize70plus[i2]))/2)
        TBPI[(k-year1+1),15]<-100*(((sum(L[i1,66:75])/(psize5574[i1]-psize5564[i1]))+(sum(L[i2,66:75])/(psize5574[i1]-psize5564[i1])))/2)
        TBPI[(k-year1+1),16]<-100*(((sum(L[i1,76:Mnage])/psize75plus[i1])+(sum(L[i2,76:Mnage])/psize75plus[i2]))/2)
        TBPI[(k-year1+1),17]<-100*(((sum(L[i1,16:25])/psize1524[i1])+(sum(L[i2,16:25])/psize1524[i2]))/2)
        
        ## (7) % reactivation
        TBRa[(k-year1+1),1]<-sum(new_actv_react[i1:i2,])
        TBRa[(k-year1+1),2]<-sum(new_actv_react[i1:i2,1:15])
        TBRa[(k-year1+1),3]<-sum(new_actv_react[i1:i2,16:30])
        TBRa[(k-year1+1),4]<-sum(new_actv_react[i1:i2,31:45])
        TBRa[(k-year1+1),5]<-sum(new_actv_react[i1:i2,46:60])
        TBRa[(k-year1+1),6]<-sum(new_actv_react[i1:i2,61:Mnage])
        
        
        ## (8) % new/reinfection
        TBRi[(k-year1+1),1]<-(sum(new_actv[i1:i2,])-sum(new_actv_react[i1:i2,])-sum(new_actvinf[i1:i2,]))
        TBRi[(k-year1+1),2]<-(sum(new_actv[i1:i2,1:15])-sum(new_actv_react[i1:i2,1:15])-sum(new_actvinf[i1:i2,1:15]))
        TBRi[(k-year1+1),3]<-(sum(new_actv[i1:i2,16:30])-sum(new_actv_react[i1:i2,16:30])-sum(new_actvinf[i1:i2,16:30]))
        TBRi[(k-year1+1),4]<-(sum(new_actv[i1:i2,31:45])-sum(new_actv_react[i1:i2,31:45])-sum(new_actvinf[i1:i2,31:45]))
        TBRi[(k-year1+1),5]<-(sum(new_actv[i1:i2,46:60])-sum(new_actv_react[i1:i2,46:60])-sum(new_actvinf[i1:i2,46:60]))
        TBRi[(k-year1+1),6]<-(sum(new_actv[i1:i2,61:Mnage])-sum(new_actv_react[i1:i2,61:Mnage])-sum(new_actvinf[i1:i2,61:Mnage]))
        
        TBRinew[(k-year1+1),1]<-sum(new_actvinf[i1:i2,])
        TBRinew[(k-year1+1),2]<-sum(new_actvinf[i1:i2,1:15])
        TBRinew[(k-year1+1),3]<-sum(new_actvinf[i1:i2,16:30])
        TBRinew[(k-year1+1),4]<-sum(new_actvinf[i1:i2,31:45])
        TBRinew[(k-year1+1),5]<-sum(new_actvinf[i1:i2,46:60])
        TBRinew[(k-year1+1),6]<-sum(new_actvinf[i1:i2,61:Mnage])
        
        
        ###(9) New active cases and TB deaths per year
        TBAc[(k-year1+1),1]<-sum(new_actv[i1:i2,],new_actvv[i1:i2,])
        
        TBAc_age[(k-year1+1),1]<-sum(new_actv[i1:i2,1:15],new_actvv[i1:i2,1:15])
        TBAc_age[(k-year1+1),2]<-sum(new_actv[i1:i2,16:25],new_actvv[i1:i2,16:25])
        TBAc_age[(k-year1+1),3]<-sum(new_actv[i1:i2,26:55],new_actvv[i1:i2,26:55])
        TBAc_age[(k-year1+1),4]<-sum(new_actv[i1:i2,56:65],new_actvv[i1:i2,56:65])
        TBAc_age[(k-year1+1),5]<-sum(new_actv[i1:i2,66:Mnage],new_actvv[i1:i2,66:Mnage])
        
        TBMo[(k-year1+1),1]<-sum(TBDeaths[i1:i2,])
        
        ##(10) number receiving vaccine
        NV[(k-year1+1),1]<-sum(num_vac[i1:i2,])
        
        TBcases[(k-year1+1),1]<-sum(I[i1:i2,1:20],NI[i1:i2,1:20],Iv[i1:i2,1:20],NIv[i1:i2,1:20])/(1/dt)
        TBcases[(k-year1+1),2]<-sum(I[i1:i2,21:60],NI[i1:i2,21:60],Iv[i1:i2,21:60],NIv[i1:i2,21:60])/(1/dt)
        TBcases[(k-year1+1),3]<-sum(I[i1:i2,61:Mnage],NI[i1:i2,61:Mnage],Iv[i1:i2,61:Mnage],NIv[i1:i2,61:Mnage])/(1/dt)
        TBcases[(k-year1+1),4]<-sum(I[i1:i2,],NI[i1:i2,],Iv[i1:i2,],NIv[i1:i2,])/(1/dt)
        TBcases[(k-year1+1),5]<-sum(I[i1:i2,1:15],NI[i1:i2,1:15],Iv[i1:i2,1:15],NIv[i1:i2,1:15])/(1/dt)
        TBcases[(k-year1+1),6]<-sum(I[i1:i2,16],NI[i1:i2,16],Iv[i1:i2,16],NIv[i1:i2,16])/(1/dt)
        TBcases[(k-year1+1),7]<-sum(I[i1:i2,17:30],NI[i1:i2,17:30],Iv[i1:i2,17:30],NIv[i1:i2,17:30])/(1/dt)
        TBcases[(k-year1+1),8]<-sum(I[i1:i2,31],NI[i1:i2,31],Iv[i1:i2,31],NIv[i1:i2,31])/(1/dt)
        TBcases[(k-year1+1),9]<-sum(I[i1:i2,32:45],NI[i1:i2,32:45],Iv[i1:i2,32:45],NIv[i1:i2,32:45])/(1/dt)
        TBcases[(k-year1+1),10]<-sum(I[i1:i2,46],NI[i1:i2,46],Iv[i1:i2,46],NIv[i1:i2,46])/(1/dt)
        TBcases[(k-year1+1),11]<-sum(I[i1:i2,47:60],NI[i1:i2,47:60],Iv[i1:i2,47:60],NIv[i1:i2,47:60])/(1/dt)
        TBcases[(k-year1+1),12]<-sum(I[i1:i2,61],NI[i1:i2,61],Iv[i1:i2,61],NIv[i1:i2,61])/(1/dt)
        TBcases[(k-year1+1),13]<-sum(I[i1:i2,62:75],NI[i1:i2,62:75],Iv[i1:i2,62:75],NIv[i1:i2,62:75])/(1/dt)
        TBcases[(k-year1+1),14]<-sum(I[i1:i2,16:75],NI[i1:i2,16:75],Iv[i1:i2,16:75],NIv[i1:i2,16:75])/(1/dt)
        
        TBcases_age[(k-year1+1),1:Mnage]<-(((I[i1, 1:Mnage] + I[i2,1:Mnage]) / (1 / dt)) + ((Iv[i1, 1:Mnage] + Iv[i2, 1:Mnage]) / (1 / dt)) +
                                             ((NI[i1, 1:Mnage] + NI[i2, 1:Mnage]) / (1 / dt)) + ((NIv[i1, 1:Mnage] +NIv[i2, 1:Mnage]) / (1 / dt)))
        TBdead[(k-year1+1),]<-colSums(TBDeaths[i1:i2,])
        
      }
      #### END OF MIDDLE YEAR RUNS
    }
  } 
  
  ## Outputs 
  assign('Sv',Sv,envir = .GlobalEnv);
  assign('S',S,envir = .GlobalEnv);assign('L',L,envir = .GlobalEnv);assign('I',I,envir = .GlobalEnv);assign('NI',NI,envir = .GlobalEnv);assign('R',R,envir = .GlobalEnv);assign('new_I',new_I,envir = .GlobalEnv);assign('new_I_noconv',new_I_noconv,envir = .GlobalEnv);assign('new_NI',new_NI,envir = .GlobalEnv);assign('new_notif',new_notif,envir = .GlobalEnv);
  assign('NBirths',BIRTHS,envir=.GlobalEnv);  assign('brate',brate,envir=.GlobalEnv)  
  assign('Lv',Lv,envir = .GlobalEnv);assign('Rv',Rv,envir = .GlobalEnv);assign('Iv',Iv,envir = .GlobalEnv);assign('NIv',NIv,envir = .GlobalEnv); assign('new_Iv',new_Iv,envir = .GlobalEnv);assign('new_Iv_noconv',new_Iv_noconv,envir = .GlobalEnv);assign('new_NIv',new_NIv,envir = .GlobalEnv);assign('new_notifv',new_notifv,envir = .GlobalEnv);
  assign('lambda',lambda,envir=.GlobalEnv);
  assign('thetaL',thetaL,envir=.GlobalEnv);
  assign('d',d,envir=.GlobalEnv);
  assign('TBDeaths',TBDeaths,envir=.GlobalEnv);
  assign('ADeaths',ADeaths,envir=.GlobalEnv);
  assign('u',u,envir=.GlobalEnv);
  assign('psize',psize,envir=.GlobalEnv);assign('bv',bv,envir=.GlobalEnv);
  assign('psize014',psize014,envir=.GlobalEnv)
  assign('psize1529',psize1529,envir=.GlobalEnv)
  assign('psize1554',psize1554,envir=.GlobalEnv)
  assign('psize1559',psize1559,envir=.GlobalEnv)
  assign('psize3044',psize3044,envir=.GlobalEnv)
  assign('psize4559',psize4559,envir=.GlobalEnv)
  assign('psize5564',psize5564,envir=.GlobalEnv)
  assign('psize60plus',psize60plus,envir=.GlobalEnv)
  assign('psize65plus',psize65plus,envir=.GlobalEnv)
  assign('psize55plus',psize55plus,envir=.GlobalEnv)
  assign('psize55minus',psize55minus,envir=.GlobalEnv)
  assign('psize0509',psize0509,envir=.GlobalEnv)
  assign('psize1019',psize1019,envir=.GlobalEnv)
  assign('psize2029',psize2029,envir=.GlobalEnv)
  assign('psize3039',psize3039,envir=.GlobalEnv)
  assign('psize4049',psize4049,envir=.GlobalEnv)
  assign('psize5059',psize5059,envir=.GlobalEnv)
  assign('psize6069',psize6069,envir=.GlobalEnv)
  assign('psize70plus',psize70plus,envir=.GlobalEnv)
  assign('psize5574',psize5574,envir=.GlobalEnv)
  assign('psize75plus',psize75plus,envir=.GlobalEnv)
  assign('psize1524',psize1524,envir=.GlobalEnv)
  assign('psize2554',psize2554,envir=.GlobalEnv)
  assign('psize15plus',psize15plus,envir=.GlobalEnv)
  assign('psizeVX',psizeVX,envir=.GlobalEnv)
  
  assign('TBI',TBI,envir=.GlobalEnv);assign('TBN',TBN,envir=.GlobalEnv);assign('TBM',TBM,envir=.GlobalEnv);
  assign('VX',VX,envir=.GlobalEnv);
  assign('TBP',TBP,envir=.GlobalEnv);assign('TBPb',TBPb,envir=.GlobalEnv)
  assign('TBPI',TBPI,envir=.GlobalEnv);assign('PSIZEy',PSIZEy,envir=.GlobalEnv);
  assign('TBRa',TBRa,envir=.GlobalEnv);assign('TBRi',TBRi,envir=.GlobalEnv); 
  assign('TBRinew',TBRinew,envir=.GlobalEnv);
  assign('TBInew',TBInew,envir=.GlobalEnv);assign('TBAc',TBAc,envir=.GlobalEnv);assign('TBAc_age',TBAc_age,envir=.GlobalEnv); assign('TBMo',TBMo,envir=.GlobalEnv);assign('num_vac',num_vac,envir=.GlobalEnv); assign('NV',NV,envir=.GlobalEnv);
  assign('CDR',CDR,envir=.GlobalEnv);
  assign('new_actv',new_actv, envir=.GlobalEnv);assign('new_actvv',new_actvv, envir=.GlobalEnv);
  assign('Imatrix',Imatrix,envir=.GlobalEnv);
  assign('psizematrix',psizematrix,envir=.GlobalEnv);
  assign('TBcases',TBcases,envir=.GlobalEnv);
  assign('TBdead',TBdead,envir=.GlobalEnv);
  assign('TBcases_age',TBcases_age,envir=.GlobalEnv);
  
  colnames(TBI)<-c("All ages","0-14", "15-54","55-64","65+","55+","<55","15-24","25-54","15+")
  colnames(TBN)<-c("All ages","0-14", "15-54", "55-64", "65+", "55+","<55","15+","15-64")
  colnames(TBM)<-c("All ages", "0-14", "15-54", "55-64", "65+", "15-59", "60+", "55+", "<55","15-24","25-54", "15-29", "30-44", "45-59")
  colnames(TBP)<-c("All ages","0-14", "15-29", "30-44", "45-59", "60+", "55+",'30-59')
  colnames(TBPb)<-c("All ages","0-14", "15-29", "30-44", "45-59", "60+", "55+","15+","15-44")
  colnames(TBPI)<-c("All ages","0-14", "15-54", "55-64", "65+", "55+", "5-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70+","6574", "75+","1524")
  colnames(PSIZEy)<-c("All ages", "0-14", "15-54", "55-64", "65+", "15-59", "15-29", "30-44", "45-59", "60+", "55+","5-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70+", "55-74", "75+","15-24", "25-54","15+",'15','60','61-74')
  
  ##  Output required 
  
  X<-cbind(psize,rowSums(S),S[,1],rowSums(I),rowSums(NI),rowSums(L),rowSums(R),rowSums(new_I),rowSums(new_NI),
           rowSums(new_I_react),rowSums(new_NI_react), 
           rowSums(Sv),rowSums(Lv),rowSums(Rv),
           TBI[,1],TBI[,2], TBI[,3], TBI[,4], TBI[,5], TBI[,6], TBI[,7], 
           TBN[,1],TBN[,2], TBN[,3], TBN[,4], TBN[,5], TBN[,6], TBN[,7],
           TBM[,1],TBM[,2],TBM[,3], TBM[,4], TBM[,5], TBM[,6], TBM[,7],TBM[,8],TBM[,9], 
           TBP[,1],TBP[,2],TBP[,3], TBP[,4], TBP[,5], TBP[,6],TBP[,7], 
           TBPb[,1],TBPb[,2],TBPb[,3], TBPb[,4], TBPb[,5], TBPb[,6],TBPb[,7],
           TBPI[,1],TBPI[,2],TBPI[,3], TBPI[,4], TBPI[,5],  TBPI[,6], TBPI[,7],TBPI[,8],TBPI[,9], TBPI[,10], TBPI[,11],  TBPI[,12], TBPI[,13], TBPI[,14],
           PSIZEy[,1],PSIZEy[,2], PSIZEy[,3], PSIZEy[,4], PSIZEy[,5], PSIZEy[,6], PSIZEy[,7], PSIZEy[,8], PSIZEy[,9], PSIZEy[,10], PSIZEy[,11], PSIZEy[,12], PSIZEy[,13], PSIZEy[,14], PSIZEy[,15], PSIZEy[,16], PSIZEy[,17], PSIZEy[,18], PSIZEy[,19],PSIZEy[,20], PSIZEy[,21], PSIZEy[,22], PSIZEy[,23],
           TBI[,8], TBI[,9],TBM[,10],TBM[,11],TBM[,12],TBM[,13],TBM[,14],TBPb[,8],PSIZEy[,24],TBRa[,1],TBRa[,2],TBRa[,3],TBRa[,4],TBRa[,5],TBRi[,1],TBRi[,2],TBRi[,3],TBRi[,4],TBRi[,5], TBPb[,9],rowSums(Iv),rowSums(NIv),rowSums(new_Iv),rowSums(new_NIv), rowSums(new_Iv_noconv), rowSums(new_notifv),TBI[,10],TBN[,8],
           TBP[,8],TBN[,9],
           TBcases[,1],TBcases[,2],TBcases[,3],TBcases[,4],
           TBcases[,5],TBcases[,6],TBcases[,7],TBcases[,8],TBcases[,9],
           TBcases[,10],TBcases[,11],TBcases[,12],TBcases[,13],TBcases[,14],
           PSIZEy[,25],PSIZEy[,26],PSIZEy[,27],
           TBRinew[,1],TBRinew[,2],TBRinew[,3],TBRinew[,4],TBRinew[,5],TBRinew[,6],
           TBRa[,1],TBRa[,2],TBRa[,3],TBRa[,4],TBRa[,5],TBRa[,6],
           TBRi[,1],TBRi[,2],TBRi[,3],TBRi[,4],TBRi[,5],TBRi[,6])
  colnames(X)<-c("PSIZE","S","Births","I","NI","L","R","new_I","new_NI","new_I_react","new_NI_react","Sv","Lv","Rv",
                 "TBItot","TBI0-14","TBI15-54","TBI55-64","TBI65+","TBI55+","TBI<55",
                 "TBNtot","TBN0-14","TBN15-54","TBN55-64","TBN65+","TBN55+","TBN<55",
                 "TBMtot","TBM0-14", "TBM15-54", "TBM55-64", "TBM65+", "TBM15-59", "TBM60+","TBM55+","TBM<55",
                 "TBPtot","TBP0-14", "TBP15-29", "TBP30-44", "TBP45-59", "TBP60+", "TBP55+",
                 "TBPbtot","TBPb0-14", "TBPb15-29", "TBPb30-44", "TBPb45-59", "TBPb60+","TBPb55+", 
                 "TBPItot","TBPI0-14", "TBPI15-54", "TBPI55-64", "TBPI65+", "TBPI55+", "TBPI5-9", "TBPI10-19", "TBPI20-29", "TBPI30-39", "TBPI40-49", "TBPI50-59", "TBPI60-69", "TBPI70+",
                 "YearPsizetot", "YearPsize0-14", "YearPsize15-54", "YearPsize55-64", "YearPsize65+", "YearPsize15-59", "YearPsize15-29", "YearPsize30-44", "YearPsize45-59", "YearPsize60+", "YearPsize55+", "YearPsize5-9", "YearPsize10-19", "YearPsize20-29", "YearPsize30-39", "YearPsize40-49", "YearPsize50-59", "YearPsize60-69", "YearPsize70+","YearPsize55-74", "YearPsize75+", "YearPsize15-24","YearPsize25-54",
                 "TBI15-24","TBI25-54","TBM15-24","TBM25-54", "TBM15-29", "TBM30-44", "TBM45-59","TBPb15+","YearPsize15plus","TBRatot","TBRa0-14","TBRa15-64","TBRa55-64","TBRa65+","TBRitot","TBRi0-14","TBRi15-64","TBRi55-64","TBRi65+", "TBPb15-44","Iv","NIv","new_Iv","new_NIv","new_Iv_noconv","new_notifv","TBI15+","TBN15+",
                 'TBP30-59',"TBN15-64",'TBcase0-19','TBcase20-59','TBcase60+','TBcasetot',
                 'TBcase0-14','TBcase15','TBcase16-29','TBcase30','TBcase31-44',
                 'TBcase45','TBcase46-59','TBcase60','TBcase61-74','TBcase15-74',
                 "YearPsize15","YearPsize60","YearPsize61-74",
                 "TBRinewtot","TBRinew0-14","TBRinew15-29","TBRinew30-44","TBRinew45-59","TBRinew60+",
                 "TBRatot","TBRa0-14","TBRa15-29","TBRa30-44","TBRa45-59","TBRa60+",
                 "TBRitot","TBRi0-14","TBRi15-29","TBRi30-44","TBRi45-59","TBRi60+")
  
  return(X)
}