# Copyright (c) 2018 Kelvin Say
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.


# The coming disruption: the movement towards the customer renewable energy transition
# Energy Policy, 2018
# ------------------------------------------------------------------------------------
# Figure 1: (600x400): plot_dayprofile(technical_data, 1, 42, 2000, 5000)
# Figure 2: (530x400): plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
#                      plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
#                      plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
#                      plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
# Figure 3: (600x300): casestudy_plotfig_a(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"))
#                      casestudy_plotfig_pv_range(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
#                      casestudy_plotfig_battery_range(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# Figure 4: (600x300): casestudy_plotfig_b(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
#                      casestudy_plotfig_c(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# Figure 5: (600x300): casestudy_plotfig_a(list(economic_2_full, economic_2_partial, economic_2_none), c("Full\nRebate", "Partial\nRebate", "No\nRebate"))
#                      casestudy_plotfig_pv_range(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
#                      casestudy_plotfig_battery_range(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
# Figure 6: (600x300): casestudy_plotfig_b(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
#                      casestudy_plotfig_c(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))

# plot_dayprofile(technical_data, 1, 42, 2000, 5000)
plot_dayprofile <- function(technical_data, year, day_of_year, pv_capacity, battery_capacity) {

  # Search for a valid PV and Battery Index
  pv_index <- 0
  battery_index <- 0
  for (i in 1:nrow(technical_data$Operational)) {
    if (pv_capacity == technical_data$Operational[i,1][[1]]$'PV Capacity (Wp)') {
      pv_index <- i
      break
    }
  }
  for (j in 1:ncol(technical_data$Operational)) {
    if (battery_capacity == technical_data$Operational[1,j][[1]]$'Battery Capacity (Wh)') {
      battery_index <- j
      break
    }
  }

  # Normal Operation
  if ((pv_index > 0) & (battery_index > 0) & (day_of_year > 0) & (day_of_year <= 365) & (year <= technical_data$Inputs$Timeline$Info$'Operational Years')) {

    time_step <- technical_data$Inputs$Timeline$Info$'Time Step (s)'
    start_index <- (day_of_year-1)*24*3600/time_step +1
    end_index <- day_of_year*24*3600/time_step
    n_samples <- end_index - start_index + 1

    timeline <- technical_data$Inputs$Timeline
    timeline$Data$'Time' <- timeline$Data$'Time'[start_index:end_index]

    baseload <- technical_data$Inputs$Baseload
    baseload$Data$'Load (W)' <- baseload$Data$'Load (W)'[start_index:end_index]

    solar <- technical_data$Inputs$Solar
    solar$Info$'Nominal Capacity (Wp)' <- pv_capacity
    solar$Data$'Insolation Reference (W/Wp)' <- solar$Data$'Insolation Reference (W/Wp)'[start_index:end_index]
    if (pv_capacity > 0) {
      gradient <- (solar$Info$'25 Year Capacity (%)' - 1) / (25*365)
      solar_degradation_start <- ((year - 1) * 365 + (day_of_year-1)) * gradient + 1
      solar_degradation_end <- ((year - 1) * 365 + day_of_year) * gradient + 1
      solar_degradation <- solar_degradation_start + (solar_degradation_end - solar_degradation_start)/n_samples*0:(n_samples-1)
      solar$Data$'Generation (W)' <- pv_capacity * solar$Data$'Insolation Reference (W/Wp)' * solar_degradation
    } else {
      solar$Data$'Generation (W)' <- NULL
    }

    battery <- technical_data$Inputs$Battery
    battery$Info$'Rated Capacity (Wh)' <- battery_capacity
    if (battery_capacity > 0) {
      gradient <- (battery$Info$'10 Year Capacity (%)' - 1) / (10*365)
      capacity_start_Wh <- battery_capacity * (((year - 1) * 365 + (day_of_year-1)) * gradient + 1)
      capacity_end_Wh <- battery_capacity * (((year - 1) * 365 + day_of_year) * gradient + 1)
      v_capacity_Wh <- capacity_start_Wh + (capacity_end_Wh - capacity_start_Wh)/n_samples*0:(n_samples-1)
      battery$Data$'Capacity (Wh)' = v_capacity_Wh
      battery$Initial$'SoC (Wh)' = technical_data$Operational[pv_index, battery_index][[1]]$Daily$'B.Soc.Initial (Wh)'[(year-1)*365 + day_of_year]
    } else {
      battery$Data$'Capacity (Wh)' = NULL
      battery$Initial$'SoC (Wh)' = 0
    }

    pv_battery_operation <- process_PV_Battery_Load(timeline, baseload, solar, battery)
    meter_operation <- data.table('Time' = timeline$Data$Time,
                                  'Load (W)' = baseload$Data$'Load (W)',
                                  'PV (W)' = solar$Data$'Generation (W)',
                                  'NetLoad (W)' = pv_battery_operation$'NetLoad (W)',
                                  'PVUsed (W)' = pv_battery_operation$'PVUsed (W)',
                                  'B.Charging (W)' = pv_battery_operation$'B.Charging (W)',
                                  'B.Discharging (W)' = pv_battery_operation$'B.Discharging (W)',
                                  'B.InvLosses (W)' = pv_battery_operation$'B.InvLosses (W)',
                                  'B.CapLosses (Wh)' = pv_battery_operation$'B.CapLosses (Wh)',
                                  'B.SoC (Wh)' = pv_battery_operation$'B.SoC (Wh)',
                                  'Imported (W)' = pv_battery_operation$'Imported (W)',
                                  'Exported (W)' = pv_battery_operation$'Exported (W)')

    max_y <- max(meter_operation$'PV (W)', meter_operation$'Load (W)')/1000
    min_y <- min(meter_operation$'NetLoad (W)', 0)/1000
    max_x <- last(meter_operation$'Time')
    min_x <- first(meter_operation$'Time')

    df <- data.table(Time = meter_operation$'Time',
                     'Net Demand' = meter_operation$'NetLoad (W)'/1000,
                     'Demand' = meter_operation$'Load (W)'/1000,
                     'PV Output' = meter_operation$'PV (W)'/1000,
                     'Battery SoC' = meter_operation$'B.SoC (Wh)'/battery_capacity * (max_y-min_y) + min_y)
    melt_df <- melt(df, id="Time")

    results_linetype <- rep("", 4)
    names(results_linetype) <- levels(melt_df$variable)
    results_linetype['Net Demand'] <- 'solid'
    results_linetype['Demand'] <- 'dashed'
    results_linetype['PV Output'] <- 'longdash'
    results_linetype['Battery SoC'] <- 'twodash'

    y_break_upper <- floor((max_y*10)/5)*0.5
    y_break_lower <- ceiling((min_y*10)/5)*0.5
    y_breaks <- 0.5*((y_break_lower/0.5):(y_break_upper/0.5))

    plot1 <- ggplot(melt_df, aes(x=Time, y=value, group=variable)) +
      theme_bw() +
      theme(legend.title = element_blank(),
            legend.position = c(0.003,0.999),
            legend.justification = c(0.01,1),
            legend.background = element_rect(linetype="solid", colour = "grey"),
            legend.key.width = unit(1,"line"),
            panel.grid.major.x = element_blank(),
            panel.grid.minor.x = element_blank(),
            plot.title = element_text(hjust = 0.5),
            axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 12),
            axis.text.y = element_text(size = 12),
            axis.line.y = element_line(colour = "grey"),
            axis.line.x = element_line(colour = "grey"),
            panel.grid.major.y = element_line(colour = "grey"),
            panel.border = element_rect(colour="grey")) +
      scale_y_continuous(minor_breaks = NULL,
                         breaks = y_breaks,
                         sec.axis = sec_axis(trans = ~.*100/(max_y-min_y)-min_y*100/(max_y-min_y),
                                             name = "State of Charge (%)"),
                         expand = c(0.002,0.002)) +
      scale_x_datetime(breaks = date_breaks("2 hours"), minor_breaks = date_breaks("1 hour"), labels = date_format("%H:%M"), expand = c(0.001,0.001)) +
      geom_line(aes(linetype=variable)) +
      scale_linetype_manual(values=results_linetype) +
      #      geom_point() +
      coord_cartesian(ylim = c(min_y, max_y), xlim = c(min_x,max_x)) +
      xlab(paste0("Time of Day")) +
      ylab("Demand & PV Output (kW)")

    # Error Case: Return an empty plot
  } else {
    df <- data.frame()
    plot1 <- ggplot(df)
  }

  return(plot1)
}

# ---------------------------------------------------------------------------------------------------------------
# Display the NPV spectrum of every PV and Battery combination and highlight the maximum
# ---------------------------------------------------------------------------------------------------------------
# 530 x 400 pixels
# plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
# plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
plot_contour <- function(simulation_info, technical_summary, financial_data) {

  max_found <- FALSE
  pv_sizes <- technical_summary$'PV Capacity (Wp)' / 1000
  battery_sizes <- technical_summary$'Battery Capacity (Wh)' / 1000
  x_diff <- simulation_info$'PV Delta (Wp)' / 1000
  y_diff <- simulation_info$'Battery Delta (Wh)' / 1000

  mat_npv <- financial_data$'System.NPV'
  colnames(mat_npv) <- paste("batt=", battery_sizes, sep="")
  rownames(mat_npv) <- paste("pv=", pv_sizes, sep="")

  # Determine the level of self-sufficiency for the most profitable system
  v_max <- which((mat_npv==max(mat_npv)) & (mat_npv > 0))
  if (length(v_max) > 0) {
    max_found <- TRUE
    max_battery_index <- ceiling(v_max[1]/length(pv_sizes))
    max_pv_index <- (v_max[1]-1) %% length(pv_sizes) + 1
    self_sufficiency <- round(technical_summary$Total$'Self-Sufficiency (%)'[max_pv_index, max_battery_index] * 100, 2)

    melt_npv <- melt(mat_npv)
    names(melt_npv) <- c("pv","batt","npv")
    melt_npv$pv <- as.numeric(str_sub(melt_npv$pv, str_locate(melt_npv$pv, "=")[1,1] + 1))
    melt_npv$batt <- as.numeric(str_sub(melt_npv$batt, str_locate(melt_npv$batt, "=")[1,1] + 1))
    melt_npv$npv <- sapply(melt_npv$npv, function(x) if (x < 0) -10 else x)
  }

  # Error Condition: return an empty plot
  if (!max_found) {
    df <- data.frame()
    plot1 <- ggplot(df) +
      geom_point() +
      theme_bw() +
      xlim(0,max(pv_sizes)) + ylim(0,max(battery_sizes)) +
      xlab(paste("Solar PV Capacity (kWp)   $", round(financial_data$Inputs$Costs$'PV Cost ($/kWp)',0), "/kWp", sep="")) +
      ylab(paste("Battery Capacity (kWh)   $", round(financial_data$Inputs$Costs$'Battery Cost ($/kWh)',0), "/kWh", sep="")) +
      ggtitle(paste("System NPV ($) Over ", financial_data$Inputs$Finance$'Analysis Years', " Years", sep=""),
              subtitle=paste("Consumption = ", round(technical_summary$Total$'Load (kWh)'/financial_data$Inputs$Finance$'Analysis Years'/1000,1), " MWh/year,  "
                             , "Inverter Limit = ", financial_data$Inputs$Tariff$'FiT System Limit (Wp)' / 1000, "kW\n"
                             , "Retail Price = ", round(financial_data$Inputs$Tariff$'Import ($/kWh)' * 100,3), "c, "
                             , "Feed-in Tariff = ", round(financial_data$Inputs$Tariff$'Export ($/kWh)' * 100,3), "c, "
                             , "Discount Rate = ", financial_data$Inputs$Finance$'Discount Rate (%)' * 100, "%\n"
                             , "Import Inflation = ", financial_data$Inputs$Tariff$'Import Inflation (%)' * 100, "%, "
                             , "Export Inflation = ", financial_data$Inputs$Tariff$'Export Inflation (%)' * 100, "%\n"
                             , "Maximum Profit = none and 0% self-sufficient"
                             , sep=""))


    # Normal Condition
  } else {

    plot1 <- ggplot(melt_npv) +
      aes(x = pv, y = batt, z = npv, fill = npv) +
      geom_tile() +
      geom_contour(color = "white", alpha = 0.5) +
      coord_equal(0.5) +
      scale_fill_gradient(low = "lightgrey",
                          high = "black",
                          na.value = "transparent",
                          breaks = c(0:4*1/4*max(melt_npv$npv)),
                          labels = dollar_format(),
                          limits = c(0, max(melt_npv$npv))) +
      theme_bw() +
      theme(legend.key.width=unit(0.03, "npc"),
            legend.key.height=unit(0.145, "npc"),
            legend.spacing = unit(0.04, "npc"),
            legend.text = element_text(size = 9),
            axis.text.x = element_text(size = 12),
            axis.text.y = element_text(size = 12)) +
      xlab(paste0("Rated PV Capacity (kWp)")) +
      ylab(paste0("Battery Capacity (kWh)")) +
      guides(fill=guide_colourbar(title="NPV\n")) +
      scale_x_continuous(expand = c(0.001,0.001)) +
      scale_y_continuous(expand = c(0.001,0.001))
  }

  return(plot1)
}

# ---------------------------------------------------------------------------------------------------------------
# Display the grid-dependency and annual energy exports for systems above the 5th NPV percentile (NPV Skim Top)
# ---------------------------------------------------------------------------------------------------------------
# plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
# plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
plot_systems <- function(simulation_info, technical_summary, financial_data) {
  maxnpv <- max(financial_data$'System.NPV')
  minnpv <- Inf
  npv_skim <- simulation_info$'NPV Skim Top (%)'
  region_indexes <- which( (financial_data$'System.NPV' >= (maxnpv * (1-npv_skim))) &  (financial_data$'System.NPV' > 0) )
  num_pv <- length(technical_summary$'PV Capacity (Wp)')
  num_systems <- length(region_indexes)
  legend_type <- 0 #0 = None, #1 = PV Only, #2 = PV-Batteries Only, #3 = PV and PV-Batteries

  # Error condition: return an empty plot
  if (num_systems <= 0) {
    df <- data.frame()
    plot1 <- ggplot(df) +
      geom_point() +
      theme_bw() +
      xlim(0,10) + ylim(0,100) +
      xlab("Annual Energy Exported (MWh)") +
      ylab("Remaining Grid Consumption (%)") +
      ggtitle("Near-Optimal DER Configurations",
              subtitle = paste("Top ", simulation_info$'NPV Skim Top (%)'*100, "% of NPV\n", "Year ", year, "\n\n", sep=""))

    # Normal condition
  } else {
    maximum_index <- which(financial_data$'System.NPV'==max(financial_data$'System.NPV'))
    max_dependence <- NaN
    max_exported <- NaN

    df <- data.frame(exported = numeric(num_systems),
                     dependence = numeric(num_systems),
                     npv = numeric(num_systems),
                     cfg_type = character(num_systems),
                     stringsAsFactors = FALSE)


    for (i in 1:num_systems) {
      battery_index <- ceiling(region_indexes[i]/num_pv)
      pv_index <- (region_indexes[i]-1) %% num_pv + 1

      # Figure out if all systems in the region are either PV-only, PV-Battery or a mix of both
      if (technical_summary$'Battery Capacity (Wh)'[battery_index] == 0) {
        df$cfg_type[i] = "PV"
        if (i == 1) {
          legend_type <- 1 # PV only
        } else {
          if (legend_type == 2) {
            legend_type <- 3 # PV and PV-Batteries
          }
        }
      } else {
        df$cfg_type[i] = "PV-Battery"
        if (i == 1) {
          legend_type <- 2 # PV-Batteries Only
        } else {
          if (legend_type == 1) {
            legend_type <- 3 # PV and PV-Batteries
          }
        }
      }

      # Record the data for the plot
      df$exported[i] <- technical_summary$Total$'Exported (kWh)'[pv_index, battery_index]/1000/technical_summary$Inputs$Timeline$Info$'Operational Years'/technical_summary$Inputs$Baseload$Info$'Annual Consumption (MWh)'*100
      df$dependence[i] <- technical_summary$Total$'Dependence (%)'[pv_index, battery_index]*100
      df$npv[i] <- financial_data$'System.NPV'[pv_index, battery_index]
      if (region_indexes[i] == maximum_index[1]) {
        max_dependence <- technical_summary$Total$'Dependence (%)'[pv_index, battery_index]*100
        max_exported <- df$exported[i]
      }
      if (df$npv[i] < minnpv) {
        minnpv <- df$npv[i]
      }
    }

    # Figure out the length shapes
    if (legend_type == 1) { # All systems in the region are PV-only
      manual_shapes = c(17)
    }
    if (legend_type == 2) { # All systems in the region are PV-Battery
      manual_shapes = c(15)
    }
    if (legend_type == 3) { # All systems in the region are a mix of PV-only and PV-Battery
      manual_shapes = c(17, 15)
    }

    df[order(df$cfg_type, df$npv),]

    plot1 <- ggplot(df) +
      aes(x=exported, y=dependence, z=npv) +
      geom_point(size = 4,
                 aes(colour = npv,
                     shape = factor(cfg_type))) +
      scale_colour_gradient(low = "lightgrey",
                            high = "black",
                            breaks = c(0:4*1/4*(maxnpv-minnpv)+minnpv),
                            labels = dollar_format()) +
      scale_shape_manual(values = manual_shapes) +
      theme_bw() +
      theme(legend.key.width=unit(0.03, "npc"),
            legend.key.height=unit(0.145, "npc"),
            legend.spacing = unit(0.04, "npc"),
            legend.text = element_text(size = 9),
            axis.text  = element_text(size = 12)) +
      coord_cartesian(xlim=c(0,200), ylim=c(0,100)) +
      xlab("Energy Exported (%)") +
      ylab("Remaining Grid Consumption (%)") +
      guides(colour=guide_colourbar(title="NPV (95th)\n", order=2), shape= "none") +
      scale_x_continuous(expand = c(0.001,0.001)) +
      scale_y_continuous(expand = c(0.001,0.001))

    # plot1 <- ggplot(melt_npv) +
    #   aes(x = pv, y = batt, z = npv, fill = npv) +
    #   geom_tile() +
    #   coord_equal(0.5) +
    #   scale_fill_gradient(low = "lightgrey",
    #                       high = "black",
    #                       na.value = "transparent",
    #                       breaks = c(0:4*1/4*max(melt_npv$npv)),
    #                       labels = dollar_format(),
    #                       limits = c(0, max(melt_npv$npv))) +
    #   theme_bw() +
    #   theme(legend.key.width=unit(0.03, "npc"),
    #         legend.key.height=unit(0.145, "npc"),
    #         legend.spacing = unit(0.04, "npc")) +
    #   xlab(paste0("Rated PV Capacity (kW)")) +
    #   ylab(paste0("Battery Capacity (kWh)")) +
    #   guides(fill=guide_colourbar(title="NPV\n")) +

  }

  return(plot1)
}

# casestudy_plotfig_a(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"))
# casestudy_plotfig_a(list(economic_2_full, economic_2_partial, economic_2_none), c("Full\nRebate", "Partial\nRebate", "No\nRebate"))
casestudy_plotfig_a <- function(scenarios_econ_data, scenarios_label) {

  num_scenarios <- length(scenarios_econ_data)
  forecast_years <- NROW(scenarios_econ_data[[1]]$'Detailed Financials')
  npv_skim <- scenarios_econ_data[[1]]$Scenario$Simulation$'NPV Skim Top (%)'
  num_pv <- length(scenarios_econ_data[[1]]$'Detailed Financials'[[1]]$'PV.Costs')

  # Create each scenario's data tables
  dt_default <- data.table(year = numeric(forecast_years),
                           cfg = character(forecast_years),
                           scenario = character(forecast_years))
  dt_scenario_list <- rep(list(dt_default), num_scenarios)

  # Evaluate each scenario within this case study
  for (n in 1:num_scenarios) {

    for (i in 1:forecast_years) {
      financial_data <- scenarios_econ_data[[n]]$'Detailed Financials'[[i]]

      maxnpv <- max(financial_data$'System.NPV')
      minnpv <- Inf
      region_indexes <- which( (financial_data$'System.NPV' >= (maxnpv * (1-npv_skim))) &  (financial_data$'System.NPV' > 0) )
      num_systems <- length(region_indexes)
      legend_type <- 0 #0 = None, #1 = PV, #2 = PV-Battery, #3 = Equivalent

      if (num_systems > 0) {
        for (j in 1:num_systems) {
          battery_index <- ceiling(region_indexes[j]/num_pv)
          pv_index <- (region_indexes[j]-1) %% num_pv + 1

          # Figure out if all systems in the region are either PV-only, PV-Battery or a mix of both
          if (scenarios_econ_data[[n]]$`Detailed Financials`[[1]]$Battery.Costs[battery_index] == 0) {
            if (j == 1) {
              legend_type <- 1 # PV
            } else {
              if (legend_type == 2) {
                legend_type <- 3 # Equivalent
              }
            }
          } else {
            if (j == 1) {
              legend_type <- 2 # PV-Battery
            } else {
              if (legend_type == 1) {
                legend_type <- 3 # Equivalent
              }
            }
          }
        }

        switch(paste0(legend_type),
               '0' = {dt_scenario_list[[n]][i]$cfg <- "None"},
               '1' = {dt_scenario_list[[n]][i]$cfg <- "PV-only"},
               '2' = {dt_scenario_list[[n]][i]$cfg <- "PV-battery"},
               '3' = {dt_scenario_list[[n]][i]$cfg <- "Comparable"})
        dt_scenario_list[[n]][i]$year <- i
        dt_scenario_list[[n]][i]$scenario <- scenarios_label[n]
      }
    }
  }

  # Setup the combined results in 'dt' for ggplot
  num_rows <- forecast_years * num_scenarios
  dt <- data.table(year = numeric(num_rows),
                   cfg = character(num_rows),
                   scenario = character(num_rows))
  index_start <- 1
  index_end <- index_start + forecast_years - 1
  for (n in 1:num_scenarios) {
    dt[index_start:index_end] <- dt_scenario_list[[n]]
    index_start <- index_end + 1
    index_end <- index_start + forecast_years - 1
  }
  dt$cfg <- factor(dt$cfg, levels = c("PV-only", "Comparable", "PV-battery"))
  dt$scenario <- factor(dt$scenario, levels = rev(scenarios_label))

  plot1 <- ggplot(dt, aes(scenario)) +
    geom_bar(aes(fill = cfg),
             position = position_stack(reverse = TRUE)) +
    scale_fill_grey(start = 0.8, end = 0.37) +
    coord_flip() +
    theme_bw() +
    theme(legend.position = "top",
          legend.key.height = unit(0.04, "npc"),
          legend.box.spacing = unit(0, "npc"),
          legend.text = element_text(size = 12, colour="black"),
          axis.text.y = element_text(size = 12),
          axis.text.x = element_text(angle = 90, hjust = 0.5, vjust=-0.5, size = 12),
          plot.margin = unit(c(0,0.2,0,0),"cm")) +
    scale_x_discrete(name = "") +
    scale_y_continuous(name = "",
                       breaks = c(1:(forecast_years)),
                       labels = names(scenarios_econ_data[[1]]$'Detailed Financials'),
                       minor_breaks = NULL,
                       expand = c(0.01,0)) +
    guides(fill = guide_legend(title = ""))

  return(plot1)
}

# casestudy_plotfig_c(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# casestudy_plotfig_c(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
casestudy_plotfig_c <- function(scenarios_econ_data, scenarios_label, scenarios_line_colours, scenarios_fill_colours) {

  # Configure the plot output
  plotcfg_y_step <- 20
  plotcfg_textsize <- 12

  # Wrangle the economic data into a combined melt format
  # -----------------------------------------------------
  num_scenarios <- length(scenarios_econ_data)
  years <- scenarios_econ_data[[1]]$'Economic Outcomes'$year

  # Main line
  y <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y <- cbind(y, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$gen_nondispatch))
  }
  names(y) <- c("year", scenarios_label)
  y_melt <- melt(y, id.vars = c("year"), variable.name = "group", value.name = "y")

  # Lower range
  y_lower <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_lower <- cbind(y_lower, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$gen_nondispatch.min))
  }
  names(y_lower) <- c("year", scenarios_label)
  ymelt_lower <- melt(y_lower, id.vars = c("year"), variable.name = "group", value.name = "y_lower")

  # Upper range
  y_upper <- data.table(year = years)
  y_max <- 100
  for (n in 1:num_scenarios) {
    y_upper <- cbind(y_upper, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$gen_nondispatch.max))
    y_max <- max(y_max, scenarios_econ_data[[n]]$'Economic Outcomes'$gen_nondispatch.max)
  }
  names(y_upper) <- c("year", scenarios_label)
  ymelt_upper <- melt(y_upper, id.vars = c("year"), variable.name = "group", value.name = "y_upper")

  combined <- merge(y_melt, ymelt_lower, by = c('year', 'group'))
  combined <- merge(combined, ymelt_upper, by = c('year', 'group'))

  # Setup ggplot parameters
  # -----------------------
  # Determine breaks
  initial_year <- min(years)
  final_year <- max(years)
  y_lim_max <- ceiling(y_max/plotcfg_y_step)*plotcfg_y_step
  y_breaks <- 0:(y_lim_max/plotcfg_y_step)*plotcfg_y_step

  # Setup the colours
  colour_fill <- scenarios_fill_colours
  names(colour_fill) <- levels(combined$group)
  colour_line <- scenarios_line_colours
  names(colour_line) <- levels(combined$group)

  # Setup the line types
  results_linetype <- c('solid', 'dashed', 'twodash')
  names(results_linetype) <- levels(combined$group)

  # Generate the plot
  plot1 <- ggplot(combined,
                  aes(x = year, y = y, group = group, colour = group)) +
    scale_colour_manual(name = "",
                        values = colour_line) +
    scale_fill_manual(name = "",
                      values = colour_fill) +
    geom_ribbon(aes(ymin = y_lower, ymax = y_upper, fill = group),
                colour = NA,
                alpha = 0.3) +
    geom_line(aes(linetype=group), size = 1) +
    scale_linetype_manual(name = "",
                          values = results_linetype) +
    scale_y_continuous(name = "Energy Exported (%)",
                       breaks = y_breaks,
                       expand = c(0, 0)) +
    scale_x_continuous(name = "",
                       breaks = years,
                       minor_breaks = NULL,
                       expand = c(0.01,0)) +
    coord_cartesian(ylim = c(0,y_lim_max)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust=0.5,
                                     size = plotcfg_textsize),
          axis.text.y = element_text(size = plotcfg_textsize),
          legend.text = element_text(size = plotcfg_textsize),
          legend.position = c(0.99, 0.99),
          legend.justification = c("right","top"),
          legend.direction = "horizontal",
          legend.key.width=unit(1.5,"line"),
          legend.background = element_rect(fill = alpha('white', 0.05)),
          plot.margin = unit(c(0.2,0.2,0,0.2),"cm"))

  return(plot1)
}


# casestudy_plotfig_b(list(economic_1_high, economic_1_low, economic_1_flat), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# casestudy_plotfig_b(list(economic_2_full, economic_2_partial, economic_2_none), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
casestudy_plotfig_b <- function(scenarios_econ_data, scenarios_label, scenarios_line_colours, scenarios_fill_colours) {

  # Configure the plot output
  plotcfg_y_step <- 20
  plotcfg_textsize <- 12

  # Wrangle the economic data into a combined melt format
  # -----------------------------------------------------
  num_scenarios <- length(scenarios_econ_data)
  years <- scenarios_econ_data[[1]]$'Economic Outcomes'$year

  # Main line
  y <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y <- cbind(y, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$dependence))
  }
  names(y) <- c("year", scenarios_label)
  y_melt <- melt(y, id.vars = c("year"), variable.name = "group", value.name = "y")

  # Lower range
  y_lower <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_lower <- cbind(y_lower, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$dependence.min))
  }
  names(y_lower) <- c("year", scenarios_label)
  ymelt_lower <- melt(y_lower, id.vars = c("year"), variable.name = "group", value.name = "y_lower")

  # Upper range
  y_upper <- data.table(year = years)
  y_max <- 100
  for (n in 1:num_scenarios) {
    y_upper <- cbind(y_upper, data.table(scenarios_econ_data[[n]]$'Economic Outcomes'$dependence.max))
    y_max <- max(y_max, scenarios_econ_data[[n]]$'Economic Outcomes'$dependence.max)
  }
  names(y_upper) <- c("year", scenarios_label)
  ymelt_upper <- melt(y_upper, id.vars = c("year"), variable.name = "group", value.name = "y_upper")

  combined <- merge(y_melt, ymelt_lower, by = c('year', 'group'))
  combined <- merge(combined, ymelt_upper, by = c('year', 'group'))

  # Setup ggplot parameters
  # -----------------------
  # Determine breaks
  initial_year <- min(years)
  final_year <- max(years)
  y_lim_max <- ceiling(y_max/plotcfg_y_step)*plotcfg_y_step
  y_breaks <- 0:(y_lim_max/plotcfg_y_step)*plotcfg_y_step

  # Setup the colours
  colour_fill <- scenarios_fill_colours
  names(colour_fill) <- levels(combined$group)
  colour_line <- scenarios_line_colours
  names(colour_line) <- levels(combined$group)

  # Setup the line types
  results_linetype <- c('solid', 'dashed', 'twodash')
  names(results_linetype) <- levels(combined$group)

  # Generate the plot
  plot1 <- ggplot(combined,
                  aes(x = year, y = y, group = group, colour = group)) +
    scale_colour_manual(name = "",
                        values = colour_line) +
    scale_fill_manual(name = "",
                      values = colour_fill) +
    geom_ribbon(aes(ymin = y_lower, ymax = y_upper, fill = group),
                colour = NA,
                alpha = 0.3) +
    geom_line(aes(linetype=group), size = 1) +
    scale_linetype_manual(name = "",
                          values = results_linetype) +
    scale_y_continuous(name = "Remaining Grid Consumption (%)",
                       breaks = y_breaks,
                       expand = c(0, 0)) +
    scale_x_continuous(name = "",
                       breaks = years,
                       minor_breaks = NULL,
                       expand = c(0.01,0)) +
    coord_cartesian(ylim = c(0,y_lim_max)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust=0.5,
                                     size = plotcfg_textsize),
          axis.text.y = element_text(size = plotcfg_textsize),
          legend.text = element_text(size = plotcfg_textsize),
          legend.position = c(0.99, 0.99),
          legend.justification = c("right","top"),
          legend.direction = "horizontal",
          legend.key.width=unit(1.5,"line"),
          legend.background = element_rect(fill = alpha('white', 0.05)),
          plot.margin = unit(c(0.2,0.2,0,0.2),"cm"))

  return(plot1)
}

# ------------------------------------------------------------------------------------------------------------
obtain_pv_battery_ranges <- function(economic_summary) {
  starting_year <- economic_summary$Scenario$Simulation$'Starting Year'
  num_years <- economic_summary$Scenario$Simulation$'Number of Years'
  pv_range <- (economic_summary$Scenario$Simulation$'PV Start (Wp)' + 0:((economic_summary$Scenario$Simulation$'PV End (Wp)' - economic_summary$Scenario$Simulation$'PV Start (Wp)') / economic_summary$Scenario$Simulation$'PV Delta (Wp)')) * economic_summary$Scenario$Simulation$'PV Delta (Wp)'
  num_pv <- length(pv_range)
  battery_range <- (economic_summary$Scenario$Simulation$'Battery Start (Wh)' + 0:((economic_summary$Scenario$Simulation$'Battery End (Wh)' - economic_summary$Scenario$Simulation$'Battery Start (Wh)') / economic_summary$Scenario$Simulation$'Battery Delta (Wh)')) * economic_summary$Scenario$Simulation$'Battery Delta (Wh)'
  npv_skim <- economic_summary$Scenario$Simulation$'NPV Skim Top (%)'

  output <- data.table(year = starting_year + 0:(num_years - 1),
                       pv_capacity.min = rep(NaN, num_years),
                       pv_capacity = rep(NaN, num_years),
                       pv_capacity.max = rep(NaN, num_years),
                       battery_capacity.min = rep(NaN, num_years),
                       battery_capacity = rep(NaN, num_years),
                       battery_capacity.max = rep(NaN, num_years),
                       npv_max = rep(NaN, num_years))

  for (i in 1:num_years) {

    financial_data <- economic_summary$'Detailed Financials'[[i]]

    # Determine the PV + battery combination that has the highest NPV
    maxnpv <- max(financial_data$'System.NPV')
    if (maxnpv > 0) {
      v_max <- which(financial_data$'System.NPV' == maxnpv)
      max_battery_index <- ceiling(v_max[1]/num_pv)
      max_pv_index <- (v_max[1]-1) %% num_pv + 1

      pv_capacity <- pv_range[max_pv_index]
      battery_capacity <- battery_range[max_battery_index]
      npv_max <- maxnpv

      output$pv_capacity[i] <- pv_capacity
      output$battery_capacity[i] <- battery_capacity
      output$npv_max[i] <- npv_max
    }

    # Record the values of the region around the maximum NPV system (defined by the "skim %" parameter)
    v_region <- which( (financial_data$'System.NPV' >= (maxnpv * (1-npv_skim))) &  (financial_data$'System.NPV' > 0) )
    if (length(v_region) > 0) {
      output$pv_capacity.min[i] <- Inf
      output$battery_capacity.min[i] <- Inf
      output$pv_capacity.max[i] <- -Inf
      output$battery_capacity.max[i] <- -Inf

      for (j in 1:length(v_region)) {
        battery_index <- ceiling(v_region[j]/num_pv)
        pv_index <- (v_region[j]-1) %% num_pv + 1

        pv_capacity <- pv_range[pv_index]
        output$pv_capacity.min[i] <- min(output$pv_capacity.min[i], pv_capacity)
        output$pv_capacity.max[i] <- max(output$pv_capacity.max[i], pv_capacity)

        battery_capacity <- battery_range[battery_index]
        output$battery_capacity.min[i] <- min(output$battery_capacity.min[i], battery_capacity)
        output$battery_capacity.max[i] <- max(output$battery_capacity.max[i], battery_capacity)
      }
    }
  }
  return(output)
}

# casestudy_plotfig_pv_range(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# casestudy_plotfig_pv_range(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
casestudy_plotfig_pv_range <- function(scenarios_econ_data, scenarios_label, scenarios_line_colours, scenarios_fill_colours) {

  # Configure the plot output
  plotcfg_textsize <- 12

  # Wrangle the economic data into a combined melt format
  # -----------------------------------------------------
  pv_range <- (scenarios_econ_data[[1]]$Scenario$Simulation$'PV Start (Wp)' + 0:((scenarios_econ_data[[1]]$Scenario$Simulation$'PV End (Wp)' - scenarios_econ_data[[1]]$Scenario$Simulation$'PV Start (Wp)') / scenarios_econ_data[[1]]$Scenario$Simulation$'PV Delta (Wp)')) * scenarios_econ_data[[1]]$Scenario$Simulation$'PV Delta (Wp)'
  y_breaks_pv <- (scenarios_econ_data[[1]]$Scenario$Simulation$'PV Start (Wp)' + 0:((scenarios_econ_data[[1]]$Scenario$Simulation$'PV End (Wp)' - scenarios_econ_data[[1]]$Scenario$Simulation$'PV Start (Wp)') / (scenarios_econ_data[[1]]$Scenario$Simulation$'PV Delta (Wp)' * 2))) * scenarios_econ_data[[1]]$Scenario$Simulation$'PV Delta (Wp)' * 2 / 1000

  num_scenarios <- length(scenarios_econ_data)
  pv_battery_ranges <- vector("list", length = num_scenarios)
  for (n in 1:num_scenarios) {
    pv_battery_ranges[[n]] <- obtain_pv_battery_ranges(scenarios_econ_data[[n]])
  }
  years <- scenarios_econ_data[[1]]$'Economic Outcomes'$year

  # Main line
  y <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y <- cbind(y, data.table(pv_battery_ranges[[n]]$'pv_capacity'/1000))
  }
  names(y) <- c("year", scenarios_label)
  y_melt <- melt(y, id.vars = c("year"), variable.name = "group", value.name = "y")

  # Lower range
  y_lower <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_lower <- cbind(y_lower, data.table(pv_battery_ranges[[n]]$'pv_capacity.min'/1000))
  }
  names(y_lower) <- c("year", scenarios_label)
  ymelt_lower <- melt(y_lower, id.vars = c("year"), variable.name = "group", value.name = "y_lower")

  # Upper range
  y_upper <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_upper <- cbind(y_upper, data.table(pv_battery_ranges[[n]]$'pv_capacity.max'/1000))
  }
  names(y_upper) <- c("year", scenarios_label)
  ymelt_upper <- melt(y_upper, id.vars = c("year"), variable.name = "group", value.name = "y_upper")

  combined <- merge(y_melt, ymelt_lower, by = c('year', 'group'))
  combined <- merge(combined, ymelt_upper, by = c('year', 'group'))

  # Setup ggplot parameters
  # -----------------------
  # Determine breaks
  initial_year <- min(years)
  final_year <- max(years)
  y_lim_max <- max(pv_range) / 1000
  y_breaks <- y_breaks_pv

  # Setup the colours
  colour_fill <- scenarios_fill_colours
  names(colour_fill) <- levels(combined$group)
  colour_line <- scenarios_line_colours
  names(colour_line) <- levels(combined$group)

  # Setup the line types
  results_linetype <- c('solid', 'dashed', 'twodash')
  names(results_linetype) <- levels(combined$group)

  # Generate the plot
  plot1 <- ggplot(combined,
                  aes(x = year, y = y, group = group, colour = group)) +
    scale_colour_manual(name = "",
                        values = colour_line) +
    scale_fill_manual(name = "",
                      values = colour_fill) +
    geom_ribbon(aes(ymin = y_lower, ymax = y_upper, fill = group),
                colour = NA,
                alpha = 0.3) +
    geom_line(aes(linetype=group), size = 1) +
    scale_linetype_manual(name = "",
                          values = results_linetype) +
    scale_y_continuous(name = "PV Capacity (kWp)",
                       breaks = y_breaks,
                       expand = c(0, 0)) +
    scale_x_continuous(name = "",
                       breaks = years,
                       minor_breaks = NULL,
                       expand = c(0.01,0)) +
    coord_cartesian(ylim = c(0,y_lim_max)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust=0.5,
                                     size = plotcfg_textsize),
          axis.text.y = element_text(size = plotcfg_textsize),
          legend.text = element_text(size = plotcfg_textsize),
          legend.position = c(0.99, 0.99),
          legend.justification = c("right","top"),
          legend.direction = "horizontal",
          legend.key.width=unit(1.5,"line"),
          legend.background = element_rect(fill = alpha('white', 0.05)),
          plot.margin = unit(c(0.2,0.2,0,0.2),"cm"))

  return(plot1)
}

# casestudy_plotfig_battery_range(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# casestudy_plotfig_battery_range(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
casestudy_plotfig_battery_range <- function(scenarios_econ_data, scenarios_label, scenarios_line_colours, scenarios_fill_colours) {

  # Configure the plot output
  plotcfg_textsize <- 12

  # Wrangle the economic data into a combined melt format
  # -----------------------------------------------------
  battery_range <- (scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Start (Wh)' + 0:((scenarios_econ_data[[1]]$Scenario$Simulation$'Battery End (Wh)' - scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Start (Wh)') / scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Delta (Wh)')) * scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Delta (Wh)'
  y_break_range <- (scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Start (Wh)' + 0:((scenarios_econ_data[[1]]$Scenario$Simulation$'Battery End (Wh)' - scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Start (Wh)') / (scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Delta (Wh)' * 2))) * scenarios_econ_data[[1]]$Scenario$Simulation$'Battery Delta (Wh)' * 2 / 1000

  num_scenarios <- length(scenarios_econ_data)
  pv_battery_ranges <- vector("list", length = num_scenarios)
  for (n in 1:num_scenarios) {
    pv_battery_ranges[[n]] <- obtain_pv_battery_ranges(scenarios_econ_data[[n]])
  }
  years <- scenarios_econ_data[[1]]$'Economic Outcomes'$year

  # Main line
  y <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y <- cbind(y, data.table(pv_battery_ranges[[n]]$'battery_capacity'/1000))
  }
  names(y) <- c("year", scenarios_label)
  y_melt <- melt(y, id.vars = c("year"), variable.name = "group", value.name = "y")

  # Lower range
  y_lower <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_lower <- cbind(y_lower, data.table(pv_battery_ranges[[n]]$'battery_capacity.min'/1000))
  }
  names(y_lower) <- c("year", scenarios_label)
  ymelt_lower <- melt(y_lower, id.vars = c("year"), variable.name = "group", value.name = "y_lower")

  # Upper range
  y_upper <- data.table(year = years)
  for (n in 1:num_scenarios) {
    y_upper <- cbind(y_upper, data.table(pv_battery_ranges[[n]]$'battery_capacity.max'/1000))
  }
  names(y_upper) <- c("year", scenarios_label)
  ymelt_upper <- melt(y_upper, id.vars = c("year"), variable.name = "group", value.name = "y_upper")

  combined <- merge(y_melt, ymelt_lower, by = c('year', 'group'))
  combined <- merge(combined, ymelt_upper, by = c('year', 'group'))

  # Setup ggplot parameters
  # -----------------------
  # Determine breaks
  initial_year <- min(years)
  final_year <- max(years)
  y_lim_max <- max(battery_range) / 1000
  y_breaks <- y_break_range

  # Setup the colours
  colour_fill <- scenarios_fill_colours
  names(colour_fill) <- levels(combined$group)
  colour_line <- scenarios_line_colours
  names(colour_line) <- levels(combined$group)

  # Setup the line types
  results_linetype <- c('solid', 'dashed', 'twodash')
  names(results_linetype) <- levels(combined$group)

  # Generate the plot
  plot1 <- ggplot(combined,
                  aes(x = year, y = y, group = group, colour = group)) +
    scale_colour_manual(name = "",
                        values = colour_line) +
    scale_fill_manual(name = "",
                      values = colour_fill) +
    geom_ribbon(aes(ymin = y_lower, ymax = y_upper, fill = group),
                colour = NA,
                alpha = 0.3) +
    geom_line(aes(linetype=group), size = 1) +
    scale_linetype_manual(name = "",
                          values = results_linetype) +
    scale_y_continuous(name = "Battery Capacity (kWh)",
                       breaks = y_breaks,
                       expand = c(0, 0)) +
    scale_x_continuous(name = "",
                       breaks = years,
                       minor_breaks = NULL,
                       expand = c(0.01,0)) +
    coord_cartesian(ylim = c(0,y_lim_max)) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90,
                                     hjust = 1,
                                     vjust=0.5,
                                     size = plotcfg_textsize),
          axis.text.y = element_text(size = plotcfg_textsize),
          legend.text = element_text(size = plotcfg_textsize),
          legend.position = c(0.99, 0.99),
          legend.justification = c("right","top"),
          legend.direction = "horizontal",
          legend.key.width=unit(1.5,"line"),
          legend.background = element_rect(fill = alpha('white', 0.05)),
          plot.margin = unit(c(0.2,0.2,0,0.2),"cm"))

  return(plot1)
}


# dev.set(dev.next())
#
# pdf("Fig1.pdf", width = 7, height = 7/6*4)
# print(p)
# dev.off()
#
# p <- plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
# pdf("Fig2a.pdf", width =7, height = 5.8)
# print(p)
# dev.off()
#
# p <- plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2018`)
# pdf("Fig2b.pdf", width = 7, height = 5.8)
# print(p)
# dev.off()
#
# p <- plot_contour(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
# pdf("Fig2c.pdf", width = 7, height = 5.8)
# print(p)
# dev.off()
#
# p <- plot_systems(technical_data$Inputs$Scenario$Simulation, technical_summary, economic_summary$`Detailed Financials`$`2023`)
# pdf("Fig2d.pdf", width = 7, height = 5.8)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_a(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"))
# pdf("Fig3a.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_pv_range(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# pdf("Fig3b.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_battery_range(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# pdf("Fig3c.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_b(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# pdf("Fig4a.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_c(list(p1_c1a_high_economic_summary, p1_c1b_low_economic_summary, p1_c1c_flat_economic_summary), c("High", "Low", "Flat"), c( 'red4', 'dodgerblue3', 'forestgreen'), c('lightpink2', 'cadetblue3', 'darkolivegreen3'))
# pdf("Fig4b.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
#
# p <- casestudy_plotfig_a(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full\nRebate", "Partial\nRebate", "No\nRebate"))
# pdf("Fig5a.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_pv_range(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial Rebate", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
# pdf("Fig5b.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_battery_range(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial Rebate", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
# pdf("Fig5c.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_b(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial Rebate", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
# pdf("Fig6a.pdf", width =7, height = 3.5)
# print(p)
# dev.off()
#
# p <- casestudy_plotfig_c(list(p1_c2a_full_economic_summary, p1_c2b_partial_economic_summary, p1_c2c_none_economic_summary), c("Full Rebate", "Partial Rebate", "No Rebate"), c('forestgreen', 'dodgerblue3', 'red4'), c('darkolivegreen3', 'cadetblue3', 'lightpink2'))
# pdf("Fig6b.pdf", width =7, height = 3.5)
# print(p)
# dev.off()