text: # -*- coding: utf-8 -*-
  """
  Created on Tue Jan 21 16:43:56 2020
  
  @author: ssterl
  """
  
  ##########################################
  ######### REVUB plotting results #########
  ##########################################
  
  # REVUB model © 2019 CIREG project
  # Author: Sebastian Sterl, Vrije Universiteit Brussel
  # This code accompanies the paper "Turbines of the Caribbean: Decarbonising Suriname's electricity mix through hydro-supported integration of wind power" by Sterl et al.
  # All equation, section &c. numbers refer to the official REVUB manual (see corresponding GitHub page, https://github.com/VUB-HYDR/REVUB).
  
  import numpy as np
  import pandas as pd
  import numbers as nb
  import matplotlib.pyplot as plt
  import numpy.matlib
  
  # [set by user] select hydropower plant (starting count at zero) and year (starting count at zero) for which to display results
  plot_HPP_multiple = np.array([0])
  plot_year_multiple = 0
  
  # [set by user] select month of year (1 = Jan, 2 = Feb, &c.) and day of month, and number of days to display results
  plot_month_multiple = 1
  plot_day_month_multiple = 14
  plot_num_days_multiple = 3
  
  # [set by user] total electricity demand to be met (MW)
  P_total_av = 146.84                         # MW
  E_total_av = (1e-3)*P_total_av*hrs_day*365  # GWh/year
  P_total_hourly = P_total_av*L_norm[:,:,0]   # MW
  
  # [calculate] non-hydro-solar-wind (thermal) power contribution (difference between total and hydro-solar-wind)
  P_BAL_thermal_hourly = P_total_hourly - np.nansum(P_BAL_hydro_stable_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_flexible_hourly[:,:,plot_HPP_multiple] + P_BAL_wind_hourly[:,:,plot_HPP_multiple] + P_BAL_solar_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_RoR_hourly[:,:,plot_HPP_multiple], axis = 2)
  P_STOR_thermal_hourly = P_total_hourly - np.nansum(P_STOR_hydro_stable_hourly[:,:,plot_HPP_multiple] + P_STOR_hydro_flexible_hourly[:,:,plot_HPP_multiple] + P_STOR_wind_hourly[:,:,plot_HPP_multiple] + P_STOR_solar_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_RoR_hourly[:,:,plot_HPP_multiple] - P_STOR_pump_hourly[:,:,plot_HPP_multiple], axis = 2)
  P_BAL_thermal_hourly[P_BAL_thermal_hourly < 0] = 0
  P_STOR_thermal_hourly[P_STOR_thermal_hourly < 0] = 0
  
  # [calculate] excess (to-be-curtailed) power
  P_BAL_curtailed_hourly = np.nansum(P_BAL_hydro_stable_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_flexible_hourly[:,:,plot_HPP_multiple] + P_BAL_wind_hourly[:,:,plot_HPP_multiple] + P_BAL_solar_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_RoR_hourly[:,:,plot_HPP_multiple], axis = 2) + P_BAL_thermal_hourly - P_total_hourly
  P_STOR_curtailed_hourly = np.nansum(P_STOR_hydro_stable_hourly[:,:,plot_HPP_multiple] + P_STOR_hydro_flexible_hourly[:,:,plot_HPP_multiple] + P_STOR_wind_hourly[:,:,plot_HPP_multiple] + P_STOR_solar_hourly[:,:,plot_HPP_multiple] + P_BAL_hydro_RoR_hourly[:,:,plot_HPP_multiple] - P_STOR_pump_hourly[:,:,plot_HPP_multiple], axis = 2) + P_STOR_thermal_hourly - P_total_hourly
  
  # [preallocate] extra variables for thermal power generation assessment
  E_total_bymonth = np.zeros(shape = (months_yr,len(simulation_years)))
  E_thermal_BAL_bymonth = np.zeros(shape = (months_yr,len(simulation_years)))
  E_thermal_STOR_bymonth = np.zeros(shape = (months_yr,len(simulation_years)))
  E_curtailed_BAL_bymonth = np.zeros(shape = (months_yr,len(simulation_years)))
  E_curtailed_STOR_bymonth = np.zeros(shape = (months_yr,len(simulation_years)))
  
  # [loop] across all years in the simulation
  for y in range(len(simulation_years)):
      # [loop] across all months of the year, converting hourly values (MW or MWh/h) to GWh/month (see eq. S24, S25)
      for m in range(months_yr):
          E_total_bymonth[m,y] = 10**(-3)*np.sum(P_total_hourly[int(positions[m,y]):int(positions[m+1,y]),y])
          E_thermal_BAL_bymonth[m,y] = 10**(-3)*np.sum(P_BAL_thermal_hourly[int(positions[m,y]):int(positions[m+1,y]),y])
          E_thermal_STOR_bymonth[m,y] = 10**(-3)*np.sum(P_STOR_thermal_hourly[int(positions[m,y]):int(positions[m+1,y]),y])
          E_curtailed_BAL_bymonth[m,y] = 10**(-3)*np.sum(P_BAL_curtailed_hourly[int(positions[m,y]):int(positions[m+1,y]),y])
          E_curtailed_STOR_bymonth[m,y] = 10**(-3)*np.sum(P_STOR_curtailed_hourly[int(positions[m,y]):int(positions[m+1,y]),y])
          
  
  # [read] vector with hours in each year
  hrs_year = range(int(hrs_byyear[plot_year_multiple]))
  
  # [identify] index of day of month to plot
  plot_day_load = np.sum(days_year[range(plot_month_multiple - 1),plot_year_multiple]) + plot_day_month_multiple - 1
  
  # [strings] string arrays containing the names and abbreviations of the different months
  months_names_full = np.array(["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
  months_names_short = np.array(["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"])
  months_byyear = np.empty(shape = (months_yr,len(simulation_years)), dtype = 'object')
  
  # [arrange] create string for each month-year combination in the time series
  for y in range(len(simulation_years)):
      for m in range(months_yr):
          months_byyear[m,y] = months_names_full[m] + np.str(simulation_years[y])
  
  # [arrange] create string for each day-month-year combination in the time series
  days_bymonth_byyear = np.empty(shape = (int(np.max(days_year)), months_yr,len(simulation_years)), dtype = 'object')
  for y in range(len(simulation_years)):
      for m in range(months_yr):
          for d in range(int(days_year[m,y])):
              days_bymonth_byyear[d,m,y] = np.str(d+1) + months_names_full[m] + 'Yr' + np.str(y+1)
  
  days_bymonth_byyear_axis = (np.transpose(days_bymonth_byyear[:,:,plot_year_multiple])).ravel()
  days_bymonth_byyear_axis = list(filter(None, days_bymonth_byyear_axis))
  
  # [colours] for plotting
  colour_hydro_stable = np.array([55, 126, 184]) / 255
  colour_hydro_flexible = np.array([106, 226, 207]) / 255
  colour_solar = np.array([255, 255, 51]) / 255
  colour_wind = np.array([77, 175, 74]) / 255
  colour_hydro_RoR = np.array([100, 100, 100]) / 255
  colour_hydro_pumped = np.array([77, 191, 237]) / 255
  colour_thermal = np.array([75, 75, 75]) / 255
  colour_curtailed = np.array([200, 200, 200]) / 255
  
  
  # [figure] (cf. Fig. S4a, S9a)
  # [plot] average monthly power mix in user-selected year
  fig = plt.figure()
  area_mix_BAL_bymonth = [np.nansum(E_hydro_BAL_stable_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_hydro_BAL_flexible_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_wind_BAL_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1) - E_curtailed_BAL_bymonth[:,plot_year_multiple], np.nansum(E_solar_BAL_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_hydro_BAL_RoR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), E_thermal_BAL_bymonth[:,plot_year_multiple], E_curtailed_BAL_bymonth[:,plot_year_multiple]]/days_year[:,plot_year_multiple]*10**3/hrs_day
  labels_generation_BAL = ['Hydropower (stable)', 'Hydropower (flexible)', 'Wind power', 'Solar power', 'Hydropower (RoR)', 'Thermal', 'Curtailed VRE']
  plt.stackplot(np.array(range(months_yr)), area_mix_BAL_bymonth, labels = labels_generation_BAL, colors = [colour_hydro_stable, colour_hydro_flexible, colour_wind, colour_solar, colour_hydro_RoR, colour_thermal, colour_curtailed])
  plt.plot(np.array(range(months_yr)), E_total_bymonth[:,plot_year_multiple]/days_year[:,plot_year_multiple]*10**3/hrs_day, label = 'Total load', color = 'black', linewidth = 3)
  plt.plot(np.array(range(months_yr)), np.nansum(ELCC_BAL_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
  plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
  plt.xticks(np.array(range(months_yr)), months_names_full, rotation = 'vertical')
  plt.ylabel('Power generation (MWh/h)')
  plt.title('monthly power generation (selected year #' + str(plot_year_multiple + 1) + ', BAL)')
  plt.savefig("Total_Fig1.png", dpi = 300, bbox_inches = 'tight')
  
  
  # [figure] (cf. Fig. S4b, S9b)
  # [plot] power mix by year
  fig = plt.figure()
  E_generated_BAL_bymonth_sum = [np.nansum(np.sum(E_hydro_BAL_stable_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_hydro_BAL_flexible_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_wind_BAL_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1) - np.sum(E_curtailed_BAL_bymonth, axis = 0), np.nansum(np.sum(E_solar_BAL_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_hydro_BAL_RoR_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.sum(E_thermal_BAL_bymonth, axis = 0), np.sum(E_curtailed_BAL_bymonth, axis = 0)]
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[0], bottom = np.sum(E_generated_BAL_bymonth_sum[0:0], axis = 0), label = 'Hydropower (stable)', color = colour_hydro_stable)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[1], bottom = np.sum(E_generated_BAL_bymonth_sum[0:1], axis = 0), label = 'Hydropower (flexible)', color = colour_hydro_flexible)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[2], bottom = np.sum(E_generated_BAL_bymonth_sum[0:2], axis = 0), label = 'Wind power', color = colour_wind)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[3], bottom = np.sum(E_generated_BAL_bymonth_sum[0:3], axis = 0), label = 'Solar power', color = colour_solar)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[4], bottom = np.sum(E_generated_BAL_bymonth_sum[0:4], axis = 0), label = 'Hydropower (RoR)', color = colour_hydro_RoR)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[5], bottom = np.sum(E_generated_BAL_bymonth_sum[0:5], axis = 0), label = 'Thermal', color = colour_thermal)
  plt.bar(np.array(range(len(simulation_years))), E_generated_BAL_bymonth_sum[6], bottom = np.sum(E_generated_BAL_bymonth_sum[0:6], axis = 0), label = 'Curtailed VRE', color = colour_curtailed)
  plt.plot(np.array(range(len(simulation_years))), np.sum(E_total_bymonth, axis = 0), label = 'Total load', color = 'black', linewidth = 3)
  plt.plot(np.array(range(len(simulation_years))), np.sum(ELCC_BAL_yearly[:,plot_HPP_multiple], axis = 1)/10**3, label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
  plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
  plt.xticks(np.array(range(len(simulation_years))), np.array(range(len(simulation_years))) + 1)
  plt.xlabel('year')
  plt.ylabel('Power generation (GWh/year)')
  plt.ylim([0, np.nanmax(np.sum(E_generated_BAL_bymonth_sum, axis = 0))*1.1])
  plt.title('Multiannual generation (BAL)')
  plt.savefig("Total_Fig2.png", dpi = 300, bbox_inches = 'tight')
  
  
  # [figure] (cf. Fig. 2 main paper, Fig. S5)
  # [plot] power mix for selected days of selected month
  fig = plt.figure()
  area_mix_full = [np.nansum(P_BAL_hydro_stable_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_BAL_hydro_flexible_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_BAL_wind_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_BAL_solar_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_BAL_hydro_RoR_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), P_BAL_thermal_hourly[hrs_year,plot_year_multiple], -1*P_BAL_curtailed_hourly[hrs_year,plot_year_multiple]]
  plt.stackplot(np.array(hrs_year), area_mix_full, labels = labels_generation_BAL, colors = [colour_hydro_stable, colour_hydro_flexible, colour_wind, colour_solar, colour_hydro_RoR, colour_thermal, colour_curtailed])
  plt.plot(np.array(hrs_year), P_total_hourly[hrs_year,plot_year_multiple], label = 'Total load', color = 'black', linewidth = 3)
  plt.plot(np.array(hrs_year), np.nansum(L_followed_BAL_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
  plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
  plt.xticks(np.array(np.arange(hrs_year[0],hrs_year[-1] + hrs_day,hrs_day)), days_bymonth_byyear_axis)
  plt.xlim([hrs_day*plot_day_load, hrs_day*(plot_day_load + plot_num_days_multiple)])
  plt.ylim([0, np.nanmax(np.sum(area_mix_full, axis = 0)*1.1)])
  plt.xlabel('Day of the year')
  plt.ylabel('Power generation (MWh/h)')
  plt.title('Daily generation & load profiles (BAL)')
  plt.savefig("Total_Fig3.png", dpi = 300, bbox_inches = 'tight')
  
  
  # [check] if STOR scenario available
  if option_storage == 1 and np.min(STOR_break[plot_HPP_multiple]) == 0:
      
      # [figure] (cf. Fig. S4a, S9a)
      # [plot] average monthly power mix in user-selected year
      fig = plt.figure()
      area_mix_STOR_bymonth = [np.nansum(E_hydro_STOR_stable_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_hydro_STOR_flexible_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_wind_STOR_bymonth[:,plot_year_multiple,plot_HPP_multiple] - E_hydro_pump_STOR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1) - E_curtailed_STOR_bymonth[:,plot_year_multiple], np.nansum(E_solar_STOR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), np.nansum(E_hydro_BAL_RoR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), E_thermal_STOR_bymonth[:,plot_year_multiple], E_curtailed_STOR_bymonth[:,plot_year_multiple]]/days_year[:,plot_year_multiple]*10**3/hrs_day
      labels_generation_STOR = ['Hydropower (stable)', 'Hydropower (flexible)', 'Wind power', 'Solar power', 'Hydropower (RoR)', 'Thermal', 'Curtailed VRE']
      plt.stackplot(np.array(range(months_yr)), area_mix_STOR_bymonth, labels = labels_generation_STOR, colors = [colour_hydro_stable, colour_hydro_flexible, colour_wind, colour_solar, colour_hydro_RoR, colour_thermal, colour_curtailed])
      plt.fill_between(np.array(range(months_yr)), -1*np.nansum(E_hydro_pump_STOR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), label = 'Stored VRE', facecolor = colour_hydro_pumped)
      plt.plot(np.array(range(months_yr)), E_total_bymonth[:,plot_year_multiple]/days_year[:,plot_year_multiple]*10**3/hrs_day, label = 'Total load', color = 'black', linewidth = 3)
      plt.plot(np.array(range(months_yr)), np.nansum(ELCC_STOR_bymonth[:,plot_year_multiple,plot_HPP_multiple], axis = 1), label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
      plt.plot(np.array(range(months_yr)), np.zeros(months_yr), color = 'black', linewidth = 1)
      plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
      plt.xticks(np.array(range(months_yr)),months_names_full, rotation = 'vertical')
      plt.ylabel('Power generation (MWh/h)')
      plt.title('monthly power generation (selected year #' + str(plot_year_multiple + 1) + ', STOR)')
      plt.savefig("Total_Fig1_b.png", dpi = 300, bbox_inches = 'tight')
      
      
      # [figure] (cf. Fig. S4b, S9b)
      # [plot] power mix by year
      fig = plt.figure()
      E_generated_STOR_bymonth_sum = [np.nansum(np.sum(E_hydro_STOR_stable_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_hydro_STOR_flexible_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_wind_STOR_bymonth[:,:,plot_HPP_multiple] - E_hydro_pump_STOR_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1)  - np.sum(E_curtailed_STOR_bymonth, axis = 0), np.nansum(np.sum(E_solar_STOR_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.nansum(np.sum(E_hydro_BAL_RoR_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), np.sum(E_thermal_STOR_bymonth, axis = 0), np.sum(E_curtailed_STOR_bymonth, axis = 0)]
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[0], bottom = np.sum(E_generated_STOR_bymonth_sum[0:0], axis = 0), label = 'Hydropower (stable)', color = colour_hydro_stable)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[1], bottom = np.sum(E_generated_STOR_bymonth_sum[0:1], axis = 0), label = 'Hydropower (flexible)', color = colour_hydro_flexible)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[2], bottom = np.sum(E_generated_STOR_bymonth_sum[0:2], axis = 0), label = 'Wind power', color = colour_wind)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[3], bottom = np.sum(E_generated_STOR_bymonth_sum[0:3], axis = 0), label = 'Solar power', color = colour_solar)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[4], bottom = np.sum(E_generated_STOR_bymonth_sum[0:4], axis = 0), label = 'Hydropower (RoR)', color = colour_hydro_RoR)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[5], bottom = np.sum(E_generated_STOR_bymonth_sum[0:5], axis = 0), label = 'Thermal', color = colour_thermal)
      plt.bar(np.array(range(len(simulation_years))), E_generated_STOR_bymonth_sum[6], bottom = np.sum(E_generated_STOR_bymonth_sum[0:6], axis = 0), label = 'Curtailed VRE', color = colour_curtailed)
      plt.bar(np.array(range(len(simulation_years))), -1*np.nansum(np.sum(E_hydro_pump_STOR_bymonth[:,:,plot_HPP_multiple], axis = 0), axis = 1), label = 'Stored VRE', color = colour_hydro_pumped)
      plt.plot(np.array(range(len(simulation_years))), np.sum(E_total_bymonth, axis = 0), label = 'Total load', color = 'black', linewidth = 3)
      plt.plot(np.array(range(len(simulation_years))), np.sum(ELCC_STOR_yearly[:,plot_HPP_multiple], axis = 1)/10**3, label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
      plt.plot(np.array(range(len(simulation_years))), np.zeros(len(simulation_years)), color = 'black', linewidth = 1)
      plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
      plt.xticks(np.array(range(len(simulation_years))), np.array(range(len(simulation_years))) + 1)
      plt.xlabel('year')
      plt.ylabel('Power generation (GWh/year)')
      plt.ylim([np.nanmin(-1*np.sum(E_hydro_pump_STOR_bymonth[:,:,plot_HPP_multiple], axis = 0))*1.1, np.nanmax(np.sum(E_generated_STOR_bymonth_sum, axis = 0))*1.1])
      plt.title('Multiannual generation (STOR)')
      plt.savefig("Total_Fig2_b.png", dpi = 300, bbox_inches = 'tight')
      
      
      # [figure] (cf. Fig. 2 main paper, Fig. S5)
      # [plot] power mix for selected days of selected month
      fig = plt.figure()
      area_mix_full = [np.nansum(P_STOR_hydro_stable_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_STOR_hydro_flexible_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_STOR_wind_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]] - P_STOR_pump_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_STOR_solar_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), np.nansum(P_BAL_hydro_RoR_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), P_STOR_thermal_hourly[hrs_year,plot_year_multiple], -1*P_STOR_curtailed_hourly[hrs_year,plot_year_multiple]]
      plt.stackplot(np.array(hrs_year), area_mix_full, labels = labels_generation_STOR, colors = [colour_hydro_stable, colour_hydro_flexible, colour_wind, colour_solar, colour_hydro_RoR, colour_thermal, colour_curtailed])
      plt.fill_between(np.array(hrs_year), -1*np.nansum(P_STOR_pump_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), label = 'Stored VRE', color = colour_hydro_pumped)
      plt.plot(np.array(hrs_year), P_total_hourly[hrs_year,plot_year_multiple], label = 'Total load', color = 'black', linewidth = 3)
      plt.plot(np.array(hrs_year), np.nansum(L_followed_STOR_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0), label = 'ELCC$_{tot}$', color = 'black', linestyle = '--', linewidth = 3)
      plt.plot(np.array(hrs_year), np.zeros(len(hrs_year)), color = 'black', linewidth = 1)
      plt.legend(loc = 'center left', bbox_to_anchor = (1, 0.5))
      plt.xticks(np.array(np.arange(hrs_year[0],hrs_year[-1] + hrs_day,hrs_day)), days_bymonth_byyear_axis)
      plt.xlim([hrs_day*plot_day_load, hrs_day*(plot_day_load + plot_num_days_multiple)])
      plt.ylim([np.nanmin(-1*np.nansum(P_STOR_pump_hourly[hrs_year,plot_year_multiple,plot_HPP_multiple[:,np.newaxis]], axis = 0))*1.1, np.nanmax(np.sum(area_mix_full, axis = 0)*1.1)])
      plt.xlabel('Day of the year')
      plt.ylabel('Power generation (MWh/h)')
      plt.title('Daily generation & load profiles (STOR)')
      plt.savefig("Total_Fig3_b.png", dpi = 300, bbox_inches = 'tight')
  