# Daily dry matter feed intake (kg/day) of a given bovine 
FeedIntake <- function(Dig, # Feed digestibility (%)
                       LW, 
                       RefCowLW, 
                       LWG, 
                       MilkProd, # Feed digestibility (%)
                       allInp,
                       MilkFatProp = NULL,
                       sexLWMult   = 1, 
                       isDairy     = T,  # Is it a dairy animal?
                       BreedMaint  = 1,  # Multiplier maintenance factor for the animal breed
                       Activity    = 1,  # Multiplier activity factor for the animal (often used for grazing)
                       catnb       = catnb)
  {
  # if(catnb == "Heifer") browser()
  
  if(is.null(MilkFatProp)) MilkFatProp <- allInp$MilkFatProp
  
  MetLW         <- LW^0.75
  MaintCoef     <- ifelse(isDairy, allInp$CMaintDairy, allInp$CMaintBeef) 
  CorrMaintCoef <- MaintCoef * Activity * BreedMaint
  Maint         <- MetLW * CorrMaintCoef
  corrMatureLW  <- (0.8 + sexLWMult * 0.2) * RefCowLW
  LWGEnergy     <- 22.02 * (LW/corrMatureLW)^0.75
  LWGFactor     <- abs(LWG)^1.097
  LactReq       <- MilkProd * (allInp$FFM_Req + allInp$MFat_Req * MilkFatProp)
  PregGT3Prop   <- 0.6
  PregReq       <- 0.1 * 3 * PregGT3Prop * Maint
  LWGReq        <- sapply(LWG, function(aLWG) ifelse (aLWG > 0, LWGEnergy*LWGFactor, -0.8 * LWGEnergy * LWGFactor))
  ED            <- Dig/100
  REM           <- 1.123 - (0.4092*ED) + 1.126 * 0.1 * ED^2 - 0.254/ED
  REG           <- 1.164 - (0.516*ED) + 1.308 * 0.1 * ED^2 - 0.374/ED
  GE_Demand     <- ((Maint + LactReq + PregReq)/REM + LWGReq/REG)/ED
  DMI           <- GE_Demand/allInp$GEFeed
  
  return(DMI)
}

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

# 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
  
  return(N2OIPCC)
}

# 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
  
  ECRefLW <- allInp$ECRefLW # Mature empty cow reference liveweight
  
  # get vector with number of animals
  n.cat <- sapply(catNames, function(nm) unname(HerdStruct[nm]))
  
  # 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))
  
  # get vectors of reference LW and actual liveweights for each category
  vRefLWs <- c("LactCow" = ECRefLW, "DryCow" = ECRefLW, "Bull" = ECRefLW * 1.2)
  vLWs    <- c("LactCow" = ECRefLW, "DryCow" = ECRefLW, "Bull" = ECRefLW * 1.2)
  
  # get vector of milk production
  LLDays <- (allInp$CowsCalvIntv - allInp$PreCalvingDryMths)*allInp$meanDiM
  vMP    <- c("LactCow" = allInp$MilkProd/LLDays, "DryCow" = 0, "Bull" = 0)
  
  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       = vLWs[[catnb]],            # vector with values corresponding to categ (columns)
                                   RefCowLW = vRefLWs[[catnb]],         # vector with values corresponding to categ (columns)
                                   LWG      = 0,                        # vector with values corresponding to categ (columns)
                                   MilkProd = vMP[[catnb]],             # vector with values corresponding to categ (columns)
                                   GEFeed   = allInp$GEFeed,            # scalar 
                                   Ym       = allInp$Ym,                # scalar
                                   allInp   = allInp,
                                   catnb    = catnb)        
    
    CatFICH4N2O_mean <- sapply(CatFICH4N2O, mean)
    return(CatFICH4N2O_mean)
    }
  )
  
  names(FICH4N2O) <- catNames
  
  FICH4N2O <- as.data.frame(do.call(rbind, FICH4N2O))
  
  propMthYr <- allInp$meanDiM/allInp$DIYr
  
  res <- cbind.data.frame(FI = FICH4N2O$FI*n.cat*allInp$DIYr, CH4 = FICH4N2O$CH4*n.cat, N2O = FICH4N2O$N2O*n.cat) 
  
  return(res)
}

# Calculate emissions of the growing animals of the herd 
calcCohortEmissions <- function(initPop     = 100,
                                birthMth    = 1:12,
                                RefMatCowLW = 500,
                                maxAge      = 100,
                                MilkProd    = 0,
                                allInp){
  
  cohData <- getCohortMtx(initPop     = initPop, 
                          birthMth    = birthMth,
                          allInp      = allInp)
  
  FICH4N2O <- extractFICH4N2O(Dig      = as.numeric(cohData$Dig),
                              LW       = cohData$LW,
                              RefCowLW = rep(RefMatCowLW, nrow(cohData)), 
                              LWG      = cohData$LWG,
                              MilkProd = rep(MilkProd, nrow(cohData)),
                              CP       = cohData$CP,
                              GEFeed   = allInp$GEFeed,
                              Ym       = allInp$Ym,
                              allInp   = allInp,
                              catnb    = cohData$CatStr)
    
  propMthYr <- cohData$mthdays/allInp$DIYr
  
  cohData <- cbind.data.frame(cohData, FI = FICH4N2O$FI, CH4 = FICH4N2O$CH4,
                              N2O = FICH4N2O$N2O, FIMth = FICH4N2O$FI*cohData$mthdays,
                              CH4Mth = FICH4N2O$CH4*propMthYr, N2OMth = FICH4N2O$N2O*propMthYr)
  
  cohData[(cohData$CatI == 1) & (cohData$age < 3), "CH4"] <- 0
  
  return(cohData)
}

extractFICH4N2O <- function(Dig,
                            LW,
                            RefCowLW,
                            LWG,
                            MilkProd,
                            CP,
                            GEFeed,
                            Ym,
                            allInp,
                            catnb,
                            ...){
  
  FI <- multiApply(FUN = FeedIntake, Dig = Dig, LW = LW, RefCowLW = RefCowLW, LWG = LWG,
                   MilkProd = MilkProd, catnb = catnb, MoreArgs = list(allInp = allInp))
  
  CH4 <- multiApply(FUN = MethaneEmission, GEI = FI*GEFeed, Ym = Ym, DIYr = allInp$DIYr,
                    CH4Energy = allInp$CH4Energy)
  
  N2O <- multiApply(FUN = N2OManure, DMI = FI, LW = LW, CP = CP, MoreArgs = list(allInp = allInp))
  
  return(list(FI = FI, CH4 = CH4, N2O = N2O))
}

getTotalHerdEmissions <- function(CvSsnalDist  = NULL, 
                                  retFemCalves = NULL, 
                                  CohFun       = calcCohortEmissions,
                                  HerdStruct,
                                  HerdStats,
                                  allInp)
  {
  
  if(is.null(CvSsnalDist)) CvSsnalDist <- allInp$CvSsnalDist
  if(is.null(retFemCalves)) retFemCalves <- HerdStats[["retFemCalves"]]
  
  # calc emissions of growing animals
  GrwFcohGEEData <- multiApply(FUN = CohFun, birthMth = 1:12, 
                               initPop = CvSsnalDist * retFemCalves, SIMPLIFY = F,
                               MoreArgs = list(allInp = allInp))
  
  names(GrwFcohGEEData) <- allInp$mthnames
  # browser()
  # get aggregate GEE emissions of growing animals
  totalsCohMth    <- GrwFcohMth2(GrwFcohGEEData, allInp = allInp)
  totalsCohMth2   <- lapply(totalsCohMth, function(x) do.call(cbind, x))
  totalsGrwCat    <- lapply(totalsCohMth2, function(x) rowSums(x))
  totalsGrwCatFmt <- do.call(cbind.data.frame, totalsGrwCat)
  
  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)
  
  names(totalsGrwCatFmt) <- names(totalMatCat)

  cats <- data.frame(Category = c(rownames(totalMatCat), rownames(totalsGrwCatFmt)))
  res  <- rbind.data.frame(totalMatCat, totalsGrwCatFmt)
  res  <- cbind.data.frame(cats, 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)
  {
  FeedProdGHGHa <- list(
    CO2LimingPerHa = c(allInp$LimePast, allInp$LimeSuppl) *  
      allInp$EfCO2Liming * allInp$CO2_C * 1e+3,                        # Mg/ha/year * 1000 = kg/ha/year
    CO2UreaPerHa = c(allInp$NFertPast, allInp$NFertSuppl) * 
      c(allInp$PropUreaPast, allInp$PropUreaSuppl) * allInp$EfCO2Urea, # kg/ha/year
    N2OFertilizer =  c(allInp$EfFertPast, allInp$EfFertOther) * 
      c(allInp$NFertPast, allInp$NFertSuppl) / 1e+3                    # g/ha/year / 1000 = kg/ha/year
  )
  
  FeedProdGHGTot <- lapply(FeedProdGHGHa, function(x) x*c(allInp$pastArea, allInp$supplArea))
  names(FeedProdGHGTot) <- c("CO2Limestone", "CO2Urea", "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]][-1]))
    
    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]){
    eq <- HerdGHGTotal[[h]]/multSimpl[h]
    yrPerHead  <- HerdGHGTotal[[h]]/allHerdStructure$Equilibrium
    dayPerHead <- yrPerHead/allInp$DIYr
    cur <- (yrPerHead * allHerdStructure$Current)/multSimpl[h]
    
    dfList[[h]] <- data.frame(Category = HerdGHGTotal$Category,
                              YearlyPerHead = yrPerHead,
                              DailyPerHead  = dayPerHead,
                              TotalEquilibriumHerd = eq,
                              TotalCurrentHerd     = cur)
  }
  
  return(dfList)
}