# Daily dry matter feed intake (kg/day) of a given bovine 
FeedIntake <- function(Dig, # Feed digestibility (%)
                       LW, 
                       catnb = "")
  {
  
  if(catnb == "LactCow"){
    DMI <- ((5.4 * LW)/500) / ((100 - Dig)/100)
    
    return(DMI)
  }
  
  NEma <- 7.5
  DMI  <- LW^0.75 * (0.2444 * NEma - 0.0111 * NEma^2 - 0.472)/NEma # kg DMI/day/animal

  return(DMI)
}

# IPCC Tier 2 annual methane (CH4) emissions (kg/yr)  dry matter feed intake (kg/day) of a given bovine 
MethaneEmission <- function(FI, allInp)  # food intake (MJ/day)
  {
  
  GEI <- FI*allInp$GEFeed
  CH4 <- (GEI * allInp$Ym * allInp$DIYr)/allInp$CH4Energy # IPCC
  CH4 <- CH4/allInp$DIYr # kg CH4/day/animal
  
  return(CH4)
}

# IPCC Tier 2 annual nitrous oxide (N2O) emissions (kg/yr) of a given bovine 
N2OManure <- function(DMI, # Dry matter intake (kg/day)
                      LW,  # Animal liveweight (kg)
                      CP,  # Crude protein concentration (% of mass)
                      allInp)
  {
  
  MetLW   <- LW^0.75
  NIntake <- DMI * (CP/100) * allInp$NinCP
  FecalN  <- NIntake * 0.7 * (1-allInp$CPDig) + 0.0038 * MetLW * allInp$NinCP
  UrineN  <- NIntake * 0.3
  TotalN  <- FecalN + UrineN
  NExcrYr <- TotalN * allInp$DIYr
  N2OIPCC <- NExcrYr * allInp$N2O_N * allInp$EfManure # IPCC
  
  N2O <- N2OIPCC/allInp$DIYr # kg N2O/day/animal
  
  return(N2O)
}

# Calculate emissions of the mature categories of a herd with a given herd structure  
calcSteadyCatEmissions <- function(HerdStruct, # Numeric vector with number of animals in each category
                                   catNames,
                                   allInp)
  {

  Diet     <- allInp$Diet
  Category <- allInp$Category
  
  # get vectors of Digestibility and Crude Protein concentration of each of category diet in each period
  lDiets        <- lapply(catNames, function(catName) GetCatDiet(aCateg = catName, allInp = allInp))
  names(lDiets) <- catNames
  
  dig <- sapply(catNames, function(catName) getDietAttrMths(diets = lDiets[[catName]], 
                                                            attribute = "Dig", mthnames = allInp$mthnames))
  
  CP <- sapply(catNames, function(catName) getDietAttrMths(diets = lDiets[[catName]], 
                                                           attribute = "CP", mthnames = allInp$mthnames))
  
  FICH4N2O <- lapply(catNames, function(catnb){
    
    CatFICH4N2O <- extractFICH4N2O(Dig      = as.numeric(dig[,catnb]),   # DF with categ as columns and mth as rows 
                                   CP       = as.numeric(CP[,catnb]),    # DF with categ as columns and mth as rows
                                   LW       = allInp$averageLW[[catnb]], # vector with values corresponding to categ (columns)
                                   allInp   = allInp,
                                   catnb    = catnb)        
    
    CatFICH4N2O_mean <- sapply(CatFICH4N2O, mean)
    
    return(CatFICH4N2O_mean)
    }
  )
  
  names(FICH4N2O) <- catNames
  
  res <- as.data.frame(do.call(rbind, FICH4N2O))
  
  return(res)
}

# Calculate emissions of the growing animals of the herd 
calcCohortEmissions <- function(initPop  = 100,
                                birthMth = 1:12,
                                maxAge   = 100,
                                allInp,
                                second   = F){
  
  cohData <- getCohortMtx(initPop  = initPop, 
                          birthMth = birthMth,
                          allInp   = allInp)
  
  if(second){
    LW <- unname(sapply(cohData$CatStr, function(nm) allInp$averageLW[[nm]]))
    catnb <- paste0(cohData$CatStr, "FI")
  }else{
    LW <- cohData$LW 
    catnb <- cohData$CatStr
  }
  
  FICH4N2O <- extractFICH4N2O(Dig      = as.numeric(cohData$Dig),
                              LW       = LW,
                              CP       = cohData$CP,
                              allInp   = allInp,
                              catnb    = catnb)
  
  cohData <- cbind.data.frame(cohData, FI = FICH4N2O$FI, CH4 = FICH4N2O$CH4, N2O = FICH4N2O$N2O)  
  
  cohData[(cohData$CatStr == "SucklingCalves") & (cohData$age < 3), "CH4"] <- 0
  
  return(cohData)
}

extractFICH4N2O <- function(Dig,
                            LW,
                            CP,
                            allInp,
                            catnb,
                            ...){
  
  FI <- multiApply(FUN = FeedIntake, Dig = Dig, LW = LW, catnb = catnb)
 
  CH4 <- multiApply(FUN = MethaneEmission, FI = FI, MoreArgs = list(allInp = allInp))
  
  N2O <- multiApply(FUN = N2OManure, DMI = FI, LW = LW, CP = CP, MoreArgs = list(allInp = allInp))
  
  return(list(FI = FI, CH4 = CH4, N2O = N2O))
}

getTotalHerdEmissions <- function(CohFun = calcCohortEmissions,
                                  HerdStruct,
                                  allInp)
  {
  
  # calc emissions of growing animals 
  GrwFcohGEEData <- multiApply(FUN = CohFun, birthMth = 1:12,
                               initPop = allInp$CvSsnalDist * HerdStruct[["WeanedCalves"]], SIMPLIFY = F,
                               MoreArgs = list(allInp = allInp, second = T))
  
  names(GrwFcohGEEData) <- allInp$mthnames
  
  # get aggregate GEE emissions of growing animals
  totalsGrwCatFmt <- GrwFcohMth(cohGEEData = GrwFcohGEEData, allInp = allInp)
  
  catnames = c("LactCow", "DryCow")
  
  if(allInp$BullsPerCow > 0){
    catnames = c(catnames, "Bull")
  }
  
  # get aggregate GEE emissions of mature animals
  totalMatCat <- calcSteadyCatEmissions(HerdStruct,
                                        catNames = catnames,
                                        allInp   = allInp)
  
  res  <- rbind(totalMatCat, totalsGrwCatFmt)
  res  <- cbind(data.frame(Category = rownames(res)), res)
  rownames(res) <- NULL
  
  if(allInp$BullsPerCow == 0){
    res <- rbind(res, data.frame(Category = "Bull", FI = NA, CH4 = NA, N2O = NA))
  }
  
  res <- res[order(res$Category),]
  
  return(res)
}

FeedProdGHG <- function(allInp){
  # Mg/ha/year C * 1000 * 44/12 = kg/ha/year CO2     # IPCC 11.12
  CO2Limestone <- c("past"  = allInp$LimePast * allInp$pastArea * allInp$EfCO2Liming,      
                    "suppl" = allInp$LimeSuppl * allInp$supplArea * allInp$EfCO2Liming)
  CO2Limestone <- CO2Limestone * allInp$CO2_C * 1e+3
  
  # kg/ha/year C * 44/12 = kg/ha/year CO2     # IPCC 11.13
  CO2Urea <- c("past"  = allInp$NFertPast * allInp$pastArea * allInp$PropUreaPast/100 * allInp$EfCO2Urea,   
               "suppl" = allInp$NFertSuppl * allInp$supplArea * allInp$PropUreaSuppl/100 * allInp$EfCO2Urea)
  CO2Urea <- CO2Urea * allInp$CO2_C
  
  # g/ha/year / 1000 = kg/ha/year
  N2OFertilizer <- c("past"  = allInp$NFertPast * allInp$pastArea * allInp$EfFertPast, 
                     "suppl" = allInp$NFertSuppl * allInp$supplArea * allInp$EfFertOther)  
  N2OFertilizer <- N2OFertilizer/1e+3 
  
  FeedProdGHGTot <- list(CO2Limestone  = CO2Limestone,
                         CO2Urea       = CO2Urea,
                         N2OFertilizer = N2OFertilizer)

  return(FeedProdGHGTot)
}

GHGSynth <- function(allHerdStructure,
                     emissionList,
                     wasteManage,
                     soilEmissions,
                     allInp)
  {
 
  GHGlist <- list()
  
  FeedProdGHGTot <- FeedProdGHG(allInp = allInp)
  v_feedGHG      <- sapply(FeedProdGHGTot, function(x) sum(x)) 
  v_feedGHG2     <- c(CO2 = v_feedGHG[["CO2Limestone"]] + v_feedGHG[["CO2Urea"]], CH4 = 0, N2O = v_feedGHG[["N2OFertilizer"]]) 
  
  ### Waste Management ###
  for(hd in c("Equilibrium", "Current")){
    colHerd <- paste0("Total", hd, "Herd")
    
    wasteManagement <- as.list(colSums(wasteManage[[hd]][-c(1,2)]))
    
    v_herdGHG <- c(CO2 = 0, CH4 = sum(emissionList$CH4[[colHerd]], na.rm = T), N2O = sum(emissionList$N2O[[colHerd]], na.rm = T))
    farmGHG   <- v_herdGHG + v_feedGHG2 
    farmCO2e  <- farmGHG * c(1, allInp$GWPCH4, allInp$GWPN2O)
    
    ### Total emissions ###
    WM_CH4  <- wasteManagement$CH4emission
    WM_DN2O <- wasteManagement$directN2OEmissions
    WM_IN2O <- wasteManagement$indirectN2OEmissions
    
    fertilizeHerdEmissions_CH4 <- v_feedGHG2[["CH4"]] + sum(emissionList$CH4[[colHerd]], na.rm = T) + WM_CH4
    fertilizeHerdEmissions_N2O <- v_feedGHG2[["N2O"]] + WM_DN2O + WM_IN2O
    fertilizeHerdEmissions_CO2 <- v_feedGHG2[["CO2"]] + soilEmissions[soilEmissions$Statistic == "annualCO2stockChange", "Total"]
    
    totalEmissions <- ((fertilizeHerdEmissions_CH4 * allInp$GWPCH4) + (fertilizeHerdEmissions_N2O * allInp$GWPN2O) + fertilizeHerdEmissions_CO2) / allInp$FPCM
    entericFermentationEmission <- (fertilizeHerdEmissions_CH4 * allInp$GWPCH4) / allInp$FPCM
    
    GHG <- list(
      feedProdApplicationEmissionCO2Limestone  = v_feedGHG[["CO2Limestone"]],
      feedProdApplicationEmissionCO2Urea       = v_feedGHG[["CO2Urea"]],
      feedProdApplicationEmissionN2OFertilizer = v_feedGHG[["N2OFertilizer"]],
      feedProdEmissionCH4  = v_feedGHG2[["CH4"]],
      feedProdEmissionN2O  = v_feedGHG2[["N2O"]],
      feedProdEmissionCO2  = v_feedGHG2[["CO2"]],
      herdWasteEmissionCH4 = v_herdGHG[["CH4"]],
      herdWasteEmissionN2O = v_herdGHG[["N2O"]],
      herdWasteEmissionCO2 = v_herdGHG[["CO2"]],
      prodSystemTotalEmissionCH4 = farmGHG[["CH4"]],
      prodSystemTotalEmissionN2O = farmGHG[["N2O"]],
      prodSystemTotalEmissionCO2 = farmGHG[["CO2"]],
      prodCO2eqTotalEmissionCH4  = farmCO2e[["CH4"]],
      prodCO2eqTotalEmissionN2O  = farmCO2e[["N2O"]],
      prodCO2eqTotalEmissionCO2  = farmCO2e[["CO2"]],
      fertilizeHerdEmissionsCH4  = fertilizeHerdEmissions_CH4,
      fertilizeHerdEmissionsN2O  = fertilizeHerdEmissions_N2O,
      fertilizeHerdEmissionsCO2  = fertilizeHerdEmissions_CO2,
      totalEmissions = totalEmissions,
      entericFermentationEmission = entericFermentationEmission
    )
    
    GHGlist[[hd]] <- c(wasteManagement, GHG)
  }
  
  resDF <- data.frame(Statistic = names(GHGlist[[1]]), Description = "", Unit = "", TotalEquilibriumHerd = as.numeric(GHGlist[["Equilibrium"]]), TotalCurrentHerd = as.numeric(GHGlist[["Current"]]))
  
  return(resDF)
}

calculateHerdEmissions <- function(HerdGHGTotal,
                                   allHerdStructure, 
                                   allInp,
                                   multSimpl){
  
  dfList <- list()
  
  for(h in colnames(HerdGHGTotal)[-1]){
    yrPerHead  <- HerdGHGTotal[[h]] * allInp$DIYr
    eq  <- (yrPerHead * allHerdStructure$Equilibrium)/multSimpl[h]
    cur <- (yrPerHead * allHerdStructure$Current)/multSimpl[h]
    
    dfList[[h]] <- data.frame(Category = HerdGHGTotal$Category,
                              DailyPerHead  = HerdGHGTotal[[h]],
                              YearlyPerHead = yrPerHead,
                              TotalEquilibriumHerd = eq,
                              TotalCurrentHerd     = cur)
  }
  
  return(dfList)
}

GECalc <- function(LW, Dig, LWG, catnb, N, allInp){
  
  coefLWG  <- ifelse(allInp$Category[[catnb]][["Sex"]] %in% "F", 0.8, 1.2)
  coefAct  <- ifelse(allInp$Category[[catnb]][["isStabled"]] %in% "yes", 0, 0.17)
  MilkProd <- ifelse(catnb %in% "LactCow", (allInp$annualMilkProdL/allInp$DIYr)/N, 0)
  
  Maint   <- LW^0.75 * 0.386 # IPCC 10.3
  LWGReq  <- 22.02 * (LW/(coefLWG * allInp$ECRefLW))^0.75 * LWG^1.097 # IPCC 10.6
  LactReq <- MilkProd * (allInp$FFM_Req + allInp$MFat_Req * allInp$MilkFatProp) # IPCC 10.8
  PregReq <- ifelse(catnb %in% c("SucklingCalves", "WeanedCalves"), 0, 0.1 * Maint) # IPCC 10.8
  WorkReq <- 0.1 * Maint * allInp$workHours # IPCC 10.11 
  ActReq  <- coefAct * Maint # IPCC 10.4
  REM     <- 1.123 - (4.092 * 10^(-3) * Dig) + (1.126 * 10^(-5) * Dig^2) - (25.4/Dig) # IPCC 10.14
  REG     <- 1.164 - (5.16 * 10^(-3) * Dig) + (1.308 * 10^(-5) * Dig^2) - (37.4/Dig)  # IPCC 10.15
  GE      <- ((Maint + LactReq + PregReq + WorkReq + ActReq)/REM + LWGReq/REG)/(Dig/100) # IPCC 10.16
  
  return(GE)
}