#########################################################################
#												                                                #
#				                     Cost Analysis                              #
#												                                                #
#########################################################################

# Author: George Gardner
# Date Created: 06/06/23
# Last Updated: 11/20/23
#   Replaced 1.7% discount rate (proposed Circular A-4 Revisions) with 2% (final Circular A-4 Revisions)

# The purpose of this code is to calculate the costs to WWTP associated with
# achieving a higher WQS in the fishery maintenance area (FMA) of the Delaware River.
# Changes to original analysis:
# - period of analysis is 30 years (2024 - 2053)
# - capital costs occur only in the first year (2024).
# - no costs occur for the next 4 years (2025-2028).
# - O&M costs occur over a 25-year period (2029-2053)
# These alternative costs assume all capital costs occur within the first year of the
# period of analysis (changed to 30-years), no costs occur for a 4-year period, and
# O&M costs occur over a 25-year period.

# Preliminaries
rm(list=ls())
gc()

options(stringsAsFactors=FALSE)
pacman::p_load(dplyr, tidyr, purrr, magrittr, readxl)

Scen_RunDate <-"11-20-23"

# Set file paths:
# input path:
input_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/Costs/Inputs/"
# R-generated output path:
output_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/Costs/Output/"
# functions path:
functions_path <- "C:/Users/40179/ICF/CWA Anniversary & Econ Support - General/06_Technical Support for Economic Analysis/TD3_Delaware River WQS/Costs/Functions/"

##################################
# Step 1: Read in cost data      #
##################################

# Kleinfelder cost estimates
Kleinfelder <- read_excel(paste0(input_path, "Lower_Delaware_River_Costs_5-24-2023.xlsx"))

######################################################################
# Step 2: Filter cost data to original and revised HADO scenarios    #
######################################################################
## Original HADO:
# Scenario AA15 only includes Class A' and A facilities (filter out Class B facilities).
# Scenario AA15 only includes ammonia limits of 1.5 mg/L for Class A' and 5 mg/L for Class A facilities.

Cost <- Kleinfelder %>%
  select(-(c(3,4,5))) %>%
  filter((Fac_Class == "A'" & NH3 == 1.5) | (Fac_Class == "A" & NH3 == 5))

## Revised HADO:
# Similar to Original HADO but also includes costs associated with achieving 6 mg/L DO

Cost_Rev_a <- Kleinfelder %>%
  select(-(c(3,4,5))) %>%
  filter((Fac_Class == "A'" & NH3 == 1.5) | (Fac_Class == "A" & NH3 == 5))

Cost_Rev_b <- Kleinfelder[,1:5] %>%
  filter(DO == 6)

Cost_Rev <- left_join(Cost_Rev_a, Cost_Rev_b, by='Plant') %>%
  select(-c(6))
Cost_Rev <- Cost_Rev %>%
  mutate(DO = ifelse(is.na(DO),"NA",DO)) %>%
  mutate(Capital_Cost_DO_19 = ifelse(is.na(Capital_Cost_DO_19), 0, Capital_Cost_DO_19)) %>%
  mutate(OM_Cost_DO_19 = ifelse(is.na(OM_Cost_DO_19), 0, OM_Cost_DO_19)) %>%
  mutate(Tot_Capital_Cost_19 = Capital_Cost_NH3_19 + Capital_Cost_DO_19) %>%
  mutate(Tot_OM_Cost_19 = OM_Cost_NH3_19 + OM_Cost_DO_19)

Cost_Rev <- Cost_Rev[,c(1:3, 6, 4:5, 7:10)]
colnames(Cost_Rev)[colnames(Cost_Rev) == 'Fac_Class.x'] <- 'Fac_Class'

rm(Cost_Rev_a, Cost_Rev_b)

######################################################################
# Step 3: Convert cost from 2019$ to 2022$ using the CCI             #
######################################################################
CCI <- read_excel(paste0(input_path, "Construction_Cost_Index_2023.xlsx"))

# NOTE: REPLACE WITH 2022$
Conv <- CCI$AVG.[2]/CCI$AVG.[5]

# For original HADO
Cost$Capital_Cost_NH3_22 <- Cost$Capital_Cost_NH3_19*Conv
Cost$OM_Cost_NH3_22 <- Cost$OM_Cost_NH3_19*Conv
Cost <- Cost[,-c(4,5)]

# For revised HADO
Cost_Rev$Capital_Cost_NH3_22 <- Cost_Rev$Capital_Cost_NH3_19*Conv
Cost_Rev$OM_Cost_NH3_22 <- Cost_Rev$OM_Cost_NH3_19*Conv
Cost_Rev$Capital_Cost_DO_22 <- Cost_Rev$Capital_Cost_DO_19*Conv
Cost_Rev$OM_Cost_DO_22 <- Cost_Rev$OM_Cost_DO_19*Conv
Cost_Rev$Tot_Capital_Cost_22 <- Cost_Rev$Tot_Capital_Cost_19*Conv
Cost_Rev$Tot_OM_Cost_22 <- Cost_Rev$Tot_OM_Cost_19*Conv
Cost_Rev <- Cost_Rev[c(3:4, 1:2, 5:9),]

######################################################################
# Step 4: Create a matrix containing year-to-year costs incurred for #
# each facility. Calculate the TPV and annualized value.             #
######################################################################
# Assumptions:
# Capital Costs are spread evenly over 7-year period from 2024 to 2030
# O&M Costs occur over a 25-year period (to 2055).
# Discount cash flow at 3% and 7% back to the year 2024

promulg_yr <- 2024
analysis_pd <- 30
discount_1_5 <- 0.015
discount_2 <- 0.02
discount_3 <- 0.03
discount_7 <- 0.07

# prepare a matrix of annual costs
# Original HADO
cost_tbl <- matrix(0,nrow = nrow(Cost), ncol = analysis_pd)
colnames(cost_tbl) <- c("2024","2025","2026","2027","2028","2029","2030","2031","2032","2033",
                     "2034","2035","2036","2037","2038","2039","2040","2041","2042","2043",
                     "2044","2045","2046","2047","2048","2049", "2050", "2051", "2052", "2053")
# capital cost component
cost_tbl_capital <- cost_tbl
cost_tbl_capital[,c(1)] <- Cost$Capital_Cost_NH3_22
# O&M cost component
cost_tbl_om <- cost_tbl
cost_tbl_om[,6:30] <- Cost$OM_Cost_NH3_22
# capital and O&M
cost_tbl[,c(1)] <- Cost$Capital_Cost_NH3_22
cost_tbl[,6:30] <- Cost$OM_Cost_NH3_22

# Revised HADO
cost_tbl_rev <- matrix(0,nrow = nrow(Cost), ncol = analysis_pd)
colnames(cost_tbl_rev) <- c("2024","2025","2026","2027","2028","2029","2030","2031","2032","2033",
                        "2034","2035","2036","2037","2038","2039","2040","2041","2042","2043",
                        "2044","2045","2046","2047","2048","2049", "2050", "2051", "2052", "2053")
# capital cost component
cost_tbl_rev_capital <- cost_tbl_rev
cost_tbl_rev_capital[,c(1)] <- Cost_Rev$Tot_Capital_Cost_22
# O&M cost component
cost_tbl_rev_om <- cost_tbl_rev
cost_tbl_rev_om[,6:30] <- Cost_Rev$Tot_OM_Cost_22
# capital and O&M
cost_tbl_rev[,c(1)] <- Cost_Rev$Tot_Capital_Cost_22
cost_tbl_rev[,6:30] <- Cost_Rev$Tot_OM_Cost_22

Cost_Rev_Subset = subset(Cost_Rev, select = c(Plant, Fac_Class, NH3, DO, Tot_Capital_Cost_22, Tot_OM_Cost_22))

source(paste0(functions_path, "TPV_func.R"))  
source(paste0(functions_path, "annualization.R"))

# create vector of years to apply to the TPV function
Year_vec = c(2024:2053)

# Original HADO ################################################
## Capital costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl_capital, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl_capital, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl_capital, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl_capital, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# capital costs
write.csv(plant_costs, paste0(output_path, "Alternate_Orig_HADO_Plant_Capital_Costs_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Orig_HADO_Total_Capital_Costs_",Scen_RunDate,".csv"), row.names = TRUE)

## O&M Costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl_om, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl_om, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl_om, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl_om, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# O&M costs
write.csv(plant_costs, paste0(output_path, "Alternate_Orig_HADO_Plant_OM_Costs_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Orig_HADO_Total_OM_Costs_",Scen_RunDate,".csv"), row.names = TRUE)

## Capital and O&M costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# total costs
write.csv(plant_costs, paste0(output_path, "Alternate_Orig_HADO_Plant_Costs_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Orig_HADO_Total_Costs_",Scen_RunDate,".csv"), row.names = TRUE)

# Revised HADO ################################################

## Capital costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl_rev_capital, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl_rev_capital, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl_rev_capital, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl_rev_capital, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost_Rev_Subset,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "DO Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# capital costs
write.csv(plant_costs, paste0(output_path, "Alternate_Revised_HADO_Plant_Capital_Costs_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Revised_HADO_Total_Capital_Costs_",Scen_RunDate,".csv"), row.names = TRUE)

## O&M Costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl_rev_om, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl_rev_om, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl_rev_om, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl_rev_om, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost_Rev_Subset,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "DO Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# O&M costs
write.csv(plant_costs, paste0(output_path, "Alternate_Revised_HADO_Plant_OM_Costs_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Revised_HADO_Total_OM_Costs_",Scen_RunDate,".csv"), row.names = TRUE)

## Capital and O&M costs
# TPV using 1.5% discount rate. 
TPV_1_5pct <- tot_pv(cost_tbl_rev, discount_1_5, Year_vec, promulg_yr)
# TPV using 2% discount rate. 
TPV_2pct <- tot_pv(cost_tbl_rev, discount_2, Year_vec, promulg_yr)
# TPV using 3% discount rate. 
TPV_3pct <- tot_pv(cost_tbl_rev, discount_3, Year_vec, promulg_yr)
# TPV using 7% discount rate. 
TPV_7pct <- tot_pv(cost_tbl_rev, discount_7, Year_vec, promulg_yr)
# annualization of costs using a 1.5% discount rate
AN_1_5pct <- annualize(TPV_1_5pct, discount_1_5, analysis_pd, 1) 
# annualization of costs using a 2% discount rate
AN_2pct <- annualize(TPV_2pct, discount_2, analysis_pd, 1) 
# annualization of costs using a 3% discount rate
AN_3pct <- annualize(TPV_3pct, discount_3, analysis_pd, 1) 
# annualization of costs using a 7% discount rate
AN_7pct <- annualize(TPV_7pct, discount_7, analysis_pd, 1) 

# convert to million $
TPV_1_5pct <- TPV_1_5pct/1000000
TPV_2pct <- TPV_2pct/1000000
TPV_3pct <- TPV_3pct/1000000 
TPV_7pct <- TPV_7pct/1000000
AN_1_5pct <- AN_1_5pct/1000000
AN_2pct <- AN_2pct/1000000
AN_3pct <- AN_3pct/1000000
AN_7pct <- AN_7pct/1000000

plant_costs <- cbind(Cost_Rev_Subset,TPV_1_5pct, AN_1_5pct, TPV_2pct, AN_2pct, TPV_3pct, AN_3pct, TPV_7pct, AN_7pct)
tot_costs <- c(sum(TPV_1_5pct), sum(AN_1_5pct), sum(TPV_2pct), sum(AN_2pct), sum(TPV_3pct), sum(AN_3pct), sum(TPV_7pct), sum(AN_7pct))
tot_costs <- as.data.frame(tot_costs)
names(tot_costs) <- c("Cost")
row.names(tot_costs) <- c("TPV (1.5%, million 2022$)", "total annualized value (1.5%, million 2022$)",
                          "TPV (2%, million 2022$)", "total annualized value (2%, million 2022$)",
                          "TPV (3%, million 2022$)", "total annualized value (3%, million 2022$)", 
                          "TPV (7%, million 2022$)", "total annualized value (7%, million 2022$)")

names(plant_costs) <- c("Plant", "Facility Class", "Ammonia Limit", "DO Limit", "Capital Cost (2022$)", "O&M Cost (2022$)",
                        "TPV (1.5%, million 2022$)", "annualized (1.5%, million 2022$)",
                        "TPV (2%, million 2022$)", "annualized (2%, million 2022$)",
                        "TPV (3%, million 2022$)", "annualized (3%, million 2022$)", 
                        "TPV (7%, million 2022$)", "annualized (7%, million 2022$)")

# total costs
write.csv(plant_costs, paste0(output_path, "Alternate_Revised_HADO_Plant_Costs_Rev_",Scen_RunDate,".csv"))
write.csv(tot_costs, paste0(output_path, "Alternate_Revised_HADO_Total_Costs_Rev_",Scen_RunDate,".csv"), row.names = TRUE)