diff --git a/main.py b/main.py
index 5140022c..e69de29b 100644
--- a/main.py
+++ b/main.py
@@ -1,68 +0,0 @@
-from scripts.geojson_creator import process_geojson
-from pathlib import Path
-import subprocess
-from scripts.ep_run_enrich import energy_plus_workflow
-from hub.imports.geometry_factory import GeometryFactory
-from hub.helpers.dictionaries import Dictionaries
-from hub.imports.construction_factory import ConstructionFactory
-from hub.imports.usage_factory import UsageFactory
-from hub.imports.weather_factory import WeatherFactory
-from hub.imports.results_factory import ResultFactory
-from scripts import random_assignation
-from hub.imports.energy_systems_factory import EnergySystemsFactory
-from scripts.energy_system_sizing_and_simulation_factory import EnergySystemsSimulationFactory
-from scripts.costs.cost import Cost
-from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT
-import hub.helpers.constants as cte
-from scripts.solar_angles import CitySolarAngles
-from scripts.pv_sizing_and_simulation import PVSizingSimulation
-from hub.exports.exports_factory import ExportsFactory
-# Specify the GeoJSON file path
-geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.0001)
-file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson')
-# Specify the output path for the PDF file
-output_path = (Path(__file__).parent / 'out_files').resolve()
-# Create city object from GeoJSON file
-city = GeometryFactory('geojson',
-                       path=file_path,
-                       height_field='height',
-                       year_of_construction_field='year_of_construction',
-                       function_field='function',
-                       function_to_hub=Dictionaries().montreal_function_to_hub_function).city
-# Enrich city data
-ConstructionFactory('nrcan', city).enrich()
-
-UsageFactory('nrcan', city).enrich()
-WeatherFactory('epw', city).enrich()
-ExportsFactory('sra', city, output_path).export()
-sra_path = (output_path / f'{city.name}_sra.xml').resolve()
-subprocess.run(['sra', str(sra_path)])
-ResultFactory('sra', city, output_path).enrich()
-solar_angles = CitySolarAngles(city.name,
-                               city.latitude,
-                               city.longitude,
-                               tilt_angle=45,
-                               surface_azimuth_angle=180).calculate
-energy_plus_workflow(city)
-random_assignation.call_random(city.buildings, random_assignation.residential_new_systems_percentage)
-EnergySystemsFactory('montreal_future', city).enrich()
-for building in city.buildings:
-  EnergySystemsSimulationFactory('archetype13', building=building, output_path=output_path).enrich()
-  if 'PV' in building.energy_systems_archetype_name:
-    ghi = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]]
-    pv_sizing_simulation = PVSizingSimulation(building,
-                                              solar_angles,
-                                              tilt_angle=45,
-                                              module_height=1,
-                                              module_width=2,
-                                              ghi=ghi)
-    pv_sizing_simulation.pv_output()
-for building in city.buildings:
-  costs = Cost(building=building, retrofit_scenario=SYSTEM_RETROFIT).life_cycle
-  costs.to_csv(output_path / f'{building.name}_lcc.csv')
-  (costs.loc['global_operational_costs', f'Scenario {SYSTEM_RETROFIT}'].
-   to_csv(output_path / f'{building.name}_op.csv'))
-  costs.loc['global_capital_costs', f'Scenario {SYSTEM_RETROFIT}'].to_csv(
-    output_path / f'{building.name}_cc.csv')
-  costs.loc['global_maintenance_costs', f'Scenario {SYSTEM_RETROFIT}'].to_csv(
-    output_path / f'{building.name}_m.csv')
diff --git a/scripts/energy_system_analysis_report.py b/scripts/energy_system_analysis_report.py
deleted file mode 100644
index 16620fc6..00000000
--- a/scripts/energy_system_analysis_report.py
+++ /dev/null
@@ -1,338 +0,0 @@
-import os
-import hub.helpers.constants as cte
-import matplotlib.pyplot as plt
-import random
-import matplotlib.colors as mcolors
-from matplotlib import cm
-from scripts.report_creation import LatexReport
-
-class EnergySystemAnalysisReport:
-  def __init__(self, city, output_path):
-    self.city = city
-    self.output_path = output_path
-    self.content = []
-    self.report = LatexReport('energy_system_analysis_report.tex')
-
-  def building_energy_info(self):
-
-    table_data = [
-      ["Building Name", "Year of Construction", "function", "Yearly Heating Demand (MWh)",
-       "Yearly Cooling Demand (MWh)", "Yearly DHW Demand (MWh)", "Yearly Electricity Demand (MWh)"]
-    ]
-    intensity_table_data = [["Building Name", "Total Floor Area m2", "Heating Demand Intensity kWh/m2",
-                             "Cooling Demand Intensity kWh/m2", "Electricity Intensity kWh/m2"]]
-
-    for building in self.city.buildings:
-      total_floor_area = 0
-      for zone in building.thermal_zones_from_internal_zones:
-        total_floor_area += zone.total_floor_area
-      building_data = [
-        building.name,
-        str(building.year_of_construction),
-        building.function,
-        str(format(building.heating_demand[cte.YEAR][0] / 3.6e9, '.2f')),
-        str(format(building.cooling_demand[cte.YEAR][0] / 3.6e9, '.2f')),
-        str(format(building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1e6, '.2f')),
-        str(format(
-          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0])
-          / 1e6, '.2f')),
-      ]
-      intensity_data = [
-        building.name,
-        str(format(total_floor_area, '.2f')),
-        str(format(building.heating_demand[cte.YEAR][0] / (3.6e6 * total_floor_area), '.2f')),
-        str(format(building.cooling_demand[cte.YEAR][0] / (3.6e6 * total_floor_area), '.2f')),
-        str(format(
-          (building.lighting_electrical_demand[cte.YEAR][0] + building.appliances_electrical_demand[cte.YEAR][0]) /
-          (1e3 * total_floor_area), '.2f'))
-      ]
-      table_data.append(building_data)
-      intensity_table_data.append(intensity_data)
-
-    self.report.add_table(table_data, caption='City Buildings Energy Demands')
-    self.report.add_table(intensity_table_data, caption='Energy Intensity Information')
-
-  def base_case_charts(self):
-    save_directory = self.output_path
-
-    def autolabel(bars, ax):
-      for bar in bars:
-        height = bar.get_height()
-        ax.annotate('{:.1f}'.format(height),
-                    xy=(bar.get_x() + bar.get_width() / 2, height),
-                    xytext=(0, 3),  # 3 points vertical offset
-                    textcoords="offset points",
-                    ha='center', va='bottom')
-
-    def create_hvac_demand_chart(building_names, yearly_heating_demand, yearly_cooling_demand):
-      fig, ax = plt.subplots()
-      bar_width = 0.35
-      index = range(len(building_names))
-
-      bars1 = ax.bar(index, yearly_heating_demand, bar_width, label='Yearly Heating Demand (MWh)')
-      bars2 = ax.bar([i + bar_width for i in index], yearly_cooling_demand, bar_width,
-                     label='Yearly Cooling Demand (MWh)')
-
-      ax.set_xlabel('Building Name')
-      ax.set_ylabel('Energy Demand (MWh)')
-      ax.set_title('Yearly HVAC Demands')
-      ax.set_xticks([i + bar_width / 2 for i in index])
-      ax.set_xticklabels(building_names, rotation=45, ha='right')
-      ax.legend()
-      autolabel(bars1, ax)
-      autolabel(bars2, ax)
-      fig.tight_layout()
-      plt.savefig(save_directory / 'hvac_demand_chart.jpg')
-      plt.close()
-
-    def create_bar_chart(title, ylabel, data, filename, bar_color=None):
-      fig, ax = plt.subplots()
-      bar_width = 0.35
-      index = range(len(building_names))
-
-      if bar_color is None:
-        # Generate a random color
-        bar_color = random.choice(list(mcolors.CSS4_COLORS.values()))
-
-      bars = ax.bar(index, data, bar_width, label=ylabel, color=bar_color)
-
-      ax.set_xlabel('Building Name')
-      ax.set_ylabel('Energy Demand (MWh)')
-      ax.set_title(title)
-      ax.set_xticks([i + bar_width / 2 for i in index])
-      ax.set_xticklabels(building_names, rotation=45, ha='right')
-      ax.legend()
-      autolabel(bars, ax)
-      fig.tight_layout()
-      plt.savefig(save_directory / filename)
-      plt.close()
-
-    building_names = [building.name for building in self.city.buildings]
-    yearly_heating_demand = [building.heating_demand[cte.YEAR][0] / 3.6e9 for building in self.city.buildings]
-    yearly_cooling_demand = [building.cooling_demand[cte.YEAR][0] / 3.6e9 for building in self.city.buildings]
-    yearly_dhw_demand = [building.domestic_hot_water_heat_demand[cte.YEAR][0] / 1e6 for building in
-                         self.city.buildings]
-    yearly_electricity_demand = [(building.lighting_electrical_demand[cte.YEAR][0] +
-                                  building.appliances_electrical_demand[cte.YEAR][0]) / 1e6 for building in
-                                 self.city.buildings]
-
-    create_hvac_demand_chart(building_names, yearly_heating_demand, yearly_cooling_demand)
-    create_bar_chart('Yearly DHW Demands', 'Energy Demand (MWh)', yearly_dhw_demand, 'dhw_demand_chart.jpg', )
-    create_bar_chart('Yearly Electricity Demands', 'Energy Demand (MWh)', yearly_electricity_demand,
-                     'electricity_demand_chart.jpg')
-
-  def maximum_monthly_hvac_chart(self):
-    save_directory = self.output_path
-    months = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October',
-              'November', 'December']
-    for building in self.city.buildings:
-      maximum_monthly_heating_load = []
-      maximum_monthly_cooling_load = []
-      fig, axs = plt.subplots(1, 2, figsize=(12, 6))  # Create a figure with 2 subplots side by side
-      for demand in building.heating_peak_load[cte.MONTH]:
-        maximum_monthly_heating_load.append(demand / 3.6e6)
-      for demand in building.cooling_peak_load[cte.MONTH]:
-        maximum_monthly_cooling_load.append(demand / 3.6e6)
-
-      # Plot maximum monthly heating load
-      axs[0].bar(months, maximum_monthly_heating_load, color='red')  # Plot on the first subplot
-      axs[0].set_title('Maximum Monthly Heating Load')
-      axs[0].set_xlabel('Month')
-      axs[0].set_ylabel('Load (kWh)')
-      axs[0].tick_params(axis='x', rotation=45)
-
-      # Plot maximum monthly cooling load
-      axs[1].bar(months, maximum_monthly_cooling_load, color='blue')  # Plot on the second subplot
-      axs[1].set_title('Maximum Monthly Cooling Load')
-      axs[1].set_xlabel('Month')
-      axs[1].set_ylabel('Load (kWh)')
-      axs[1].tick_params(axis='x', rotation=45)
-
-      plt.tight_layout()  # Adjust layout to prevent overlapping
-      plt.savefig(save_directory / f'{building.name}_monthly_maximum_hvac_loads.jpg')
-      plt.close()
-
-  def load_duration_curves(self):
-    save_directory = self.output_path
-    for building in self.city.buildings:
-      heating_demand = [demand / 3.6e6 for demand in building.heating_demand[cte.HOUR]]
-      cooling_demand = [demand / 3.6e6 for demand in building.cooling_demand[cte.HOUR]]
-      heating_demand_sorted = sorted(heating_demand, reverse=True)
-      cooling_demand_sorted = sorted(cooling_demand, reverse=True)
-
-      plt.style.use('ggplot')
-
-      # Create figure and axis objects with 1 row and 2 columns
-      fig, axs = plt.subplots(1, 2, figsize=(12, 6))
-
-      # Plot sorted heating demand
-      axs[0].plot(heating_demand_sorted, color='red', linewidth=2, label='Heating Demand')
-      axs[0].set_xlabel('Hour', fontsize=14)
-      axs[0].set_ylabel('Heating Demand (kWh)', fontsize=14)
-      axs[0].set_title('Heating Load Duration Curve', fontsize=16)
-      axs[0].grid(True)
-      axs[0].legend(loc='upper right', fontsize=12)
-
-      # Plot sorted cooling demand
-      axs[1].plot(cooling_demand_sorted, color='blue', linewidth=2, label='Cooling Demand')
-      axs[1].set_xlabel('Hour', fontsize=14)
-      axs[1].set_ylabel('Cooling Demand (kWh)', fontsize=14)
-      axs[1].set_title('Cooling Load Duration Curve', fontsize=16)
-      axs[1].grid(True)
-      axs[1].legend(loc='upper right', fontsize=12)
-
-      # Adjust layout
-      plt.tight_layout()
-      plt.savefig(save_directory / f'{building.name}_load_duration_curve.jpg')
-      plt.close()
-
-  def individual_building_info(self, building):
-    table_data = [
-      ["Maximum Monthly HVAC Demands",
-       f"\\includegraphics[width=1\\linewidth]{{{building.name}_monthly_maximum_hvac_loads.jpg}}"],
-      ["Load Duration Curve", f"\\includegraphics[width=1\\linewidth]{{{building.name}_load_duration_curve.jpg}}"],
-    ]
-
-    self.report.add_table(table_data, caption=f'{building.name} Information', first_column_width=1.5)
-
-  def building_system_retrofit_results(self, building_name, current_system, new_system):
-    current_system_archetype = current_system[f'{building_name}']['Energy System Archetype']
-    current_system_heating = current_system[f'{building_name}']['Heating Equipments']
-    current_system_cooling = current_system[f'{building_name}']['Cooling Equipments']
-    current_system_dhw = current_system[f'{building_name}']['DHW Equipments']
-    current_system_pv = current_system[f'{building_name}']['Photovoltaic System Capacity']
-    current_system_heating_fuel = current_system[f'{building_name}']['Heating Fuel']
-    current_system_hvac_consumption = current_system[f'{building_name}']['Yearly HVAC Energy Consumption (MWh)']
-    current_system_dhw_consumption = current_system[f'{building_name}']['DHW Energy Consumption (MWH)']
-    current_pv_production = current_system[f'{building_name}']['PV Yearly Production (kWh)']
-    current_capital_cost = current_system[f'{building_name}']['Energy System Capital Cost (CAD)']
-    current_operational = current_system[f'{building_name}']['Energy System Average Yearly Operational Cost (CAD)']
-    current_lcc = current_system[f'{building_name}']['Energy System Life Cycle Cost (CAD)']
-    new_system_archetype = new_system[f'{building_name}']['Energy System Archetype']
-    new_system_heating = new_system[f'{building_name}']['Heating Equipments']
-    new_system_cooling = new_system[f'{building_name}']['Cooling Equipments']
-    new_system_dhw = new_system[f'{building_name}']['DHW Equipments']
-    new_system_pv = new_system[f'{building_name}']['Photovoltaic System Capacity']
-    new_system_heating_fuel = new_system[f'{building_name}']['Heating Fuel']
-    new_system_hvac_consumption = new_system[f'{building_name}']['Yearly HVAC Energy Consumption (MWh)']
-    new_system_dhw_consumption = new_system[f'{building_name}']['DHW Energy Consumption (MWH)']
-    new_pv_production = new_system[f'{building_name}']['PV Yearly Production (kWh)']
-    new_capital_cost = new_system[f'{building_name}']['Energy System Capital Cost (CAD)']
-    new_operational = new_system[f'{building_name}']['Energy System Average Yearly Operational Cost (CAD)']
-    new_lcc = new_system[f'{building_name}']['Energy System Life Cycle Cost (CAD)']
-
-    energy_system_table_data = [
-      ["Detail", "Existing System", "Proposed System"],
-      ["Energy System Archetype", current_system_archetype, new_system_archetype],
-      ["Heating Equipments", current_system_heating, new_system_heating],
-      ["Cooling Equipments", current_system_cooling, new_system_cooling],
-      ["DHW Equipments", current_system_dhw, new_system_dhw],
-      ["Photovoltaic System Capacity", current_system_pv, new_system_pv],
-      ["Heating Fuel", current_system_heating_fuel, new_system_heating_fuel],
-      ["Yearly HVAC Energy Consumption (MWh)", current_system_hvac_consumption, new_system_hvac_consumption],
-      ["DHW Energy Consumption (MWH)", current_system_dhw_consumption, new_system_dhw_consumption],
-      ["PV Yearly Production (kWh)",  current_pv_production, new_pv_production],
-      ["Energy System Capital Cost (CAD)", current_capital_cost, new_capital_cost],
-      ["Energy System Average Yearly Operational Cost (CAD)", current_operational, new_operational],
-      ["Energy System Life Cycle Cost (CAD)", current_lcc, new_lcc]
-    ]
-    self.report.add_table(energy_system_table_data, caption=f'Building {building_name} Energy System Characteristics')
-
-  def building_fuel_consumption_breakdown(self, building):
-    save_directory = self.output_path
-    # Initialize variables to store fuel consumption breakdown
-    fuel_breakdown = {
-      "Heating": {"Gas": 0, "Electricity": 0},
-      "Domestic Hot Water": {"Gas": 0, "Electricity": 0},
-      "Cooling": {"Electricity": 0},
-      "Appliance": building.appliances_electrical_demand[cte.YEAR][0] / 1e6,
-      "Lighting": building.lighting_electrical_demand[cte.YEAR][0] / 1e6
-    }
-
-    # Iterate through energy systems of the building
-    for energy_system in building.energy_systems:
-      for demand_type in energy_system.demand_types:
-        if demand_type == cte.HEATING:
-          consumption = building.heating_consumption[cte.YEAR][0] / 3.6e9
-          for generation_system in energy_system.generation_systems:
-            if generation_system.fuel_type == cte.ELECTRICITY:
-              fuel_breakdown[demand_type]["Electricity"] += consumption
-            else:
-              fuel_breakdown[demand_type]["Gas"] += consumption
-        elif demand_type == cte.DOMESTIC_HOT_WATER:
-          consumption = building.domestic_hot_water_consumption[cte.YEAR][0] / 1e6
-          for generation_system in energy_system.generation_systems:
-            if generation_system.fuel_type == cte.ELECTRICITY:
-              fuel_breakdown[demand_type]["Electricity"] += consumption
-            else:
-              fuel_breakdown[demand_type]["Gas"] += consumption
-        elif demand_type == cte.COOLING:
-          consumption = building.cooling_consumption[cte.YEAR][0] / 3.6e9
-          fuel_breakdown[demand_type]["Electricity"] += consumption
-
-    electricity_labels = ['Appliance', 'Lighting']
-    electricity_sizes = [fuel_breakdown['Appliance'], fuel_breakdown['Lighting']]
-    if fuel_breakdown['Heating']['Electricity'] > 0:
-      electricity_labels.append('Heating')
-      electricity_sizes.append(fuel_breakdown['Heating']['Electricity'])
-    if fuel_breakdown['Cooling']['Electricity'] > 0:
-      electricity_labels.append('Cooling')
-      electricity_sizes.append(fuel_breakdown['Cooling']['Electricity'])
-    if fuel_breakdown['Domestic Hot Water']['Electricity'] > 0:
-      electricity_labels.append('Domestic Hot Water')
-      electricity_sizes.append(fuel_breakdown['Domestic Hot Water']['Electricity'])
-
-    # Data for bar chart
-    gas_labels = ['Heating', 'Domestic Hot Water']
-    gas_sizes = [fuel_breakdown['Heating']['Gas'], fuel_breakdown['Domestic Hot Water']['Gas']]
-
-    # Set the style
-    plt.style.use('ggplot')
-
-    # Create plot grid
-    fig, axs = plt.subplots(1, 2, figsize=(12, 6))
-
-    # Plot pie chart for electricity consumption breakdown
-    colors = cm.get_cmap('tab20c', len(electricity_labels))
-    axs[0].pie(electricity_sizes, labels=electricity_labels,
-               autopct=lambda pct: f"{pct:.1f}%\n({pct / 100 * sum(electricity_sizes):.2f})",
-               startangle=90, colors=[colors(i) for i in range(len(electricity_labels))])
-    axs[0].set_title('Electricity Consumption Breakdown')
-
-    # Plot bar chart for natural gas consumption breakdown
-    colors = cm.get_cmap('Paired', len(gas_labels))
-    axs[1].bar(gas_labels, gas_sizes, color=[colors(i) for i in range(len(gas_labels))])
-    axs[1].set_ylabel('Consumption (MWh)')
-    axs[1].set_title('Natural Gas Consumption Breakdown')
-
-    # Add grid to bar chart
-    axs[1].grid(axis='y', linestyle='--', alpha=0.7)
-
-    # Add a title to the entire figure
-    plt.suptitle('Building Energy Consumption Breakdown', fontsize=16, fontweight='bold')
-
-    # Adjust layout
-    plt.tight_layout()
-
-    # Save the plot as a high-quality image
-    plt.savefig(save_directory / f'{building.name}_energy_consumption_breakdown.png', dpi=300)
-    plt.close()
-
-  def create_report(self, current_system, new_system):
-    os.chdir(self.output_path)
-    self.report.add_section('Current Status')
-    self.building_energy_info()
-    self.base_case_charts()
-    self.report.add_image('hvac_demand_chart.jpg', caption='Yearly HVAC Demands')
-    self.report.add_image('dhw_demand_chart.jpg', caption='Yearly DHW Demands')
-    self.report.add_image('electricity_demand_chart.jpg', caption='Yearly Electricity Demands')
-    self.maximum_monthly_hvac_chart()
-    self.load_duration_curves()
-    for building in self.city.buildings:
-      self.individual_building_info(building)
-      self.building_system_retrofit_results(building_name=building.name, current_system=current_system, new_system=new_system)
-      self.building_fuel_consumption_breakdown(building)
-      self.report.add_image(f'{building.name}_energy_consumption_breakdown.png',
-                            caption=f'Building {building.name} Consumption by source and sector breakdown')
-    self.report.save_report()
-    self.report.compile_to_pdf()
diff --git a/scripts/system_simulation.py b/scripts/system_simulation.py
deleted file mode 100644
index a097d0b0..00000000
--- a/scripts/system_simulation.py
+++ /dev/null
@@ -1,135 +0,0 @@
-import csv
-import math
-from typing import List
-
-from pathlib import Path
-
-import hub.helpers.constants as cte
-from hub.helpers.monthly_values import MonthlyValues
-
-
-class SystemSimulation:
-  def __init__(self, building, out_path):
-    self.building = building
-    self.energy_systems = building.energy_systems
-    self.heating_demand = [0] + building.heating_demand[cte.HOUR]
-    self.cooling_demand = building.cooling_demand
-    self.dhw_demand = building.domestic_hot_water_heat_demand
-    self.T_out = building.external_temperature[cte.HOUR]
-    self.maximum_heating_demand = building.heating_peak_load[cte.YEAR][0]
-    self.maximum_cooling_demand = building.cooling_peak_load[cte.YEAR][0]
-    self.name = building.name
-    self.energy_system_archetype = building.energy_systems_archetype_name
-    self.out_path = out_path
-
-  def archetype1(self):
-    out_path = self.out_path
-    T, T_sup, T_ret, m_ch, m_dis, q_hp, q_aux = [0] * len(self.heating_demand), [0] * len(
-      self.heating_demand), [0] * len(self.heating_demand), [0] * len(self.heating_demand), [0] * len(
-      self.heating_demand), [0] * len(self.heating_demand), [0] * len(self.heating_demand)
-    hp_electricity: List[float] = [0.0] * len(self.heating_demand)
-    aux_fuel: List[float] = [0.0] * len(self.heating_demand)
-    heating_consumption: List[float] = [0.0] * len(self.heating_demand)
-    boiler_consumption: List[float] = [0.0] * len(self.heating_demand)
-    T[0], dt = 25, 3600  # Assuming dt is defined somewhere
-    ua, v, hp_cap, hp_efficiency, boiler_efficiency = 0, 0, 0, 0, 0
-    for energy_system in self.energy_systems:
-      if cte.ELECTRICITY not in energy_system.demand_types:
-        generation_systems = energy_system.generation_systems
-        for generation_system in generation_systems:
-          if generation_system.system_type == cte.HEAT_PUMP and cte.HEATING in energy_system.demand_types:
-            hp_cap = generation_system.nominal_heat_output
-            hp_efficiency = float(generation_system.heat_efficiency)
-            for storage in generation_system.energy_storage_systems:
-              if storage.type_energy_stored == 'thermal':
-                v, h = float(storage.volume), float(storage.height)
-                r_tot = sum(float(layer.thickness) / float(layer.material.conductivity) for layer in
-                            storage.layers)
-                u_tot = 1 / r_tot
-                d = math.sqrt((4 * v) / (math.pi * h))
-                a_side = math.pi * d * h
-                a_top = math.pi * d ** 2 / 4
-                ua = u_tot * (2 * a_top + a_side)
-          elif generation_system.system_type == cte.BOILER:
-            boiler_cap = generation_system.nominal_heat_output
-            boiler_efficiency = float(generation_system.heat_efficiency)
-
-    for i in range(len(self.heating_demand) - 1):
-      T[i + 1] = T[i] + ((m_ch[i] * (T_sup[i] - T[i])) + (
-              ua * (self.T_out[i] - T[i])) / cte.WATER_HEAT_CAPACITY - m_dis[i] * (T[i] - T_ret[i])) * (dt / (cte.WATER_DENSITY * v))
-      if T[i + 1] < 35:
-          q_hp[i + 1] = hp_cap * 1000
-          m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * 7)
-          T_sup[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + T[i + 1]
-      elif 35 <= T[i + 1] < 45 and q_hp[i] == 0:
-          q_hp[i + 1] = 0
-          m_ch[i + 1] = 0
-          T_sup[i + 1] = T[i + 1]
-      elif 35 <= T[i + 1] < 45 and q_hp[i] > 0:
-          q_hp[i + 1] = hp_cap * 1000
-          m_ch[i + 1] = q_hp[i + 1] / (cte.WATER_HEAT_CAPACITY * 3)
-          T_sup[i + 1] = (q_hp[i + 1] / (m_ch[i + 1] * cte.WATER_HEAT_CAPACITY)) + T[i + 1]
-      else:
-        q_hp[i + 1], m_ch[i + 1], T_sup[i + 1] = 0, 0, T[i + 1]
-
-      hp_electricity[i + 1] = q_hp[i + 1] / hp_efficiency
-      if self.heating_demand[i + 1] == 0:
-        m_dis[i + 1], t_return, T_ret[i + 1] = 0, T[i + 1], T[i + 1]
-      else:
-        if self.heating_demand[i + 1] > 0.5 * self.maximum_heating_demand:
-          factor = 8
-        else:
-          factor = 4
-        m_dis[i + 1] = self.maximum_heating_demand / (cte.WATER_HEAT_CAPACITY * factor * 3600)
-        t_return = T[i + 1] - self.heating_demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY * 3600)
-        if m_dis[i + 1] == 0 or (m_dis[i + 1] > 0 and t_return < 25):
-          T_ret[i + 1] = max(25, T[i + 1])
-        else:
-          T_ret[i + 1] = T[i + 1] - self.heating_demand[i + 1] / (m_dis[i + 1] * cte.WATER_HEAT_CAPACITY * 3600)
-      tes_output = m_dis[i + 1] * cte.WATER_HEAT_CAPACITY * (T[i + 1] - T_ret[i + 1])
-      if tes_output < (self.heating_demand[i + 1] / 3600):
-        q_aux[i + 1] = (self.heating_demand[i + 1] / 3600) - tes_output
-        aux_fuel[i + 1] = (q_aux[i + 1] * dt) / 35.8e6
-      boiler_consumption[i + 1] = q_aux[i + 1] / boiler_efficiency
-      heating_consumption[i + 1] = boiler_consumption[i + 1] + hp_electricity[i + 1]
-    data = list(zip(T, T_sup, T_ret, m_ch, m_dis, q_hp, hp_electricity, aux_fuel, q_aux, self.heating_demand))
-    file_name = f'simulation_results_{self.name}.csv'
-    with open(out_path / file_name, 'w', newline='') as csvfile:
-      output_file = csv.writer(csvfile)
-      # Write header
-      output_file.writerow(['T', 'T_sup', 'T_ret', 'm_ch', 'm_dis', 'q_hp', 'hp_electricity', 'aux_fuel', 'q_aux', 'heating_demand'])
-      # Write data
-      output_file.writerows(data)
-    return heating_consumption, hp_electricity, boiler_consumption, T_sup
-
-  def enrich(self):
-    if self.energy_system_archetype == 'PV+ASHP+GasBoiler+TES' or 'PV+4Pipe+DHW':
-      building_new_heating_consumption, building_heating_electricity_consumption, building_heating_gas_consumption, supply_temperature = (
-        self.archetype1())
-      self.building.heating_consumption[cte.HOUR] = building_new_heating_consumption
-      self.building.heating_consumption[cte.MONTH] = MonthlyValues.get_total_month(self.building.heating_consumption[cte.HOUR])
-      self.building.heating_consumption[cte.YEAR] = [sum(self.building.heating_consumption[cte.MONTH])]
-      disaggregated_consumption = {}
-      for energy_system in self.building.energy_systems:
-        if cte.HEATING in energy_system.demand_types:
-          for generation_system in energy_system.generation_systems:
-            if generation_system.system_type == cte.HEAT_PUMP:
-              generation_system.heat_supply_temperature = supply_temperature
-            disaggregated_consumption[generation_system.fuel_type] = {}
-            if generation_system.fuel_type == cte.ELECTRICITY:
-              disaggregated_consumption[generation_system.fuel_type][
-                cte.HOUR] = building_heating_electricity_consumption
-              disaggregated_consumption[generation_system.fuel_type][cte.MONTH] = MonthlyValues.get_total_month(
-                disaggregated_consumption[generation_system.fuel_type][cte.HOUR])
-              disaggregated_consumption[generation_system.fuel_type][cte.YEAR] = [
-                sum(disaggregated_consumption[generation_system.fuel_type][cte.MONTH])]
-            else:
-              disaggregated_consumption[generation_system.fuel_type][cte.HOUR] = building_heating_gas_consumption
-              disaggregated_consumption[generation_system.fuel_type][cte.MONTH] = MonthlyValues.get_total_month(
-                disaggregated_consumption[generation_system.fuel_type][cte.HOUR])
-              disaggregated_consumption[generation_system.fuel_type][cte.YEAR] = [
-                sum(disaggregated_consumption[generation_system.fuel_type][cte.MONTH])]
-      self.building.heating_fuel_consumption_disaggregated = disaggregated_consumption
-    return self.building
-
-