diff --git a/src/topupopt/data/buildings/dk/heat.py b/src/topupopt/data/buildings/dk/heat.py
index 72625f9f7495a517fdb3cf3ca9019cc650c3da5c..8e91f96dd18ac67ec4a16e32ffc84f096a62f516 100644
--- a/src/topupopt/data/buildings/dk/heat.py
+++ b/src/topupopt/data/buildings/dk/heat.py
@@ -4,9 +4,9 @@
 import numpy as np
 from geopandas import GeoDataFrame
 from ...misc.utils import discrete_sinusoid_matching_integral
-# from ...misc.utils import create_profile_using_time_weighted_state
+from ...misc.utils import create_profile_using_time_weighted_state
 from ...misc.utils import generate_manual_state_correlated_profile
-from ...misc.utils import generate_state_correlated_profile
+from ...misc.utils import generate_state_correlated_profile, generate_profile
 from .bbr import label_bbr_entrance_id, label_bbr_housing_area
 
 # *****************************************************************************
@@ -62,14 +62,15 @@ def heat_demand_dict_by_building_entrance(
     number_intervals: int,
     time_interval_durations: list,
     bdg_specific_demand: dict,
-    bdg_ratio_min_max: dict,
+    bdg_min_max_ratio: dict,
     bdg_demand_phase_shift: dict = None,
     key_osm_entr_id: str = label_osm_entrance_id,
     key_bbr_entr_id: str = label_bbr_entrance_id,
     avg_state: list = None,
     state_correlates_with_output: bool = False,
     deviation_gain: float = 1.0,
-    solver: str = 'glpk'
+    solver: str = 'glpk',
+    **kwargs
 ) -> dict:
     
     # initialise dict for each building entrance
@@ -89,7 +90,7 @@ def heat_demand_dict_by_building_entrance(
         for building_index in building_indexes:
             # get relevant data
             area = gdf_buildings.loc[building_index][label_bbr_housing_area]
-            # generate the profile
+            # generate the profile            
             if type(avg_state) == type(None):
                 # ignore states
                 heat_demand_profiles.append(
@@ -97,7 +98,7 @@ def heat_demand_dict_by_building_entrance(
                         discrete_sinusoid_matching_integral(
                             bdg_specific_demand[building_index] * area,
                             time_interval_durations=time_interval_durations,
-                            min_to_max_ratio=bdg_ratio_min_max[building_index],
+                            min_max_ratio=bdg_min_max_ratio[building_index],
                             phase_shift_radians=(
                                 bdg_demand_phase_shift[building_index]
                             ),
@@ -115,7 +116,7 @@ def heat_demand_dict_by_building_entrance(
                 #             ),
                 #             avg_state=avg_state,
                 #             time_interval_durations=time_interval_durations,
-                #             min_to_max_ratio=bdg_ratio_min_max[building_index],
+                #             min_max_ratio=bdg_min_max_ratio[building_index],
                 #             state_correlates_with_output=state_correlates_with_output,
                 #         )
                 #     )
@@ -129,7 +130,7 @@ def heat_demand_dict_by_building_entrance(
                             states=avg_state, 
                             time_interval_durations=time_interval_durations,
                             states_correlate_profile=state_correlates_with_output,
-                            min_max_ratio=bdg_ratio_min_max[building_index],
+                            min_max_ratio=bdg_min_max_ratio[building_index],
                             solver=solver
                             )
                         )
@@ -161,10 +162,74 @@ def heat_demand_dict_by_building_entrance(
     # return
     return demand_dict
 
-
 # *****************************************************************************
 # *****************************************************************************
 
+# TODO: allow reusing the gain
+
+def heat_demand_profiles(
+    gdf_osm: GeoDataFrame,
+    gdf_buildings: GeoDataFrame,
+    time_interval_durations: list,
+    assessments: list,
+    annual_heat_demand: dict,
+    air_temperature: dict = None, 
+    reuse_deviation_gain: bool = True,
+    **kwargs
+) -> dict:    
+    
+    # calculate the total area
+    total_area = total_heating_area(gdf_osm, gdf_buildings)
+    # initialise data dict
+    heat_demand_dict = {}
+    # for each building entrance
+    for osm_index in gdf_osm.index:
+        # initialise dict for each building entrance
+        bdg_entr_dict = {}
+        # find the indexes for each building leading to the curr. cons. point
+        building_indexes = gdf_buildings[
+            gdf_buildings[label_bbr_entrance_id] == gdf_osm.loc[osm_index][label_osm_entrance_id]
+        ].index
+        for q in assessments:
+            # define the specific heat demand
+            specific_demand = annual_heat_demand[q]/total_area
+            # 
+            # initialise dict for each building consumption point
+            bdg_profile_list = []
+            # for each building
+            for building_index in building_indexes:
+                # get relevant data
+                area = gdf_buildings.loc[building_index][label_bbr_housing_area]
+                # handle states
+                if type(air_temperature) != type(None):
+                    kwargs['states'] = air_temperature[q]
+                # append the profile for each building to the list
+                _profile = generate_profile(
+                    integration_result=specific_demand*area, 
+                    time_interval_durations=time_interval_durations, 
+                    **kwargs
+                    # min_max_ratio, 
+                    # states, 
+                    # states_correlate_profile, 
+                    # solver, 
+                    # deviation_gain
+                    )
+                bdg_profile_list.append(np.array(_profile))
+            # aggregate profiles for the same building entrance
+            bdg_entr_profile = (
+                [] 
+                if len(bdg_profile_list) == 0 else 
+                sum(profile for profile in bdg_profile_list)
+                )
+            # store the demand profile
+            bdg_entr_dict[q] = bdg_entr_profile
+        # add to the main dict
+        heat_demand_dict[osm_index] = bdg_entr_dict
+        
+    return heat_demand_dict, total_area
+
+# *****************************************************************************
+# *****************************************************************************
 
 def total_heating_area(
     gdf_osm: GeoDataFrame,
@@ -185,6 +250,5 @@ def total_heating_area(
             area += gdf_buildings.loc[building_index][label_bbr_housing_area]
     return area
 
-
 # *****************************************************************************
 # *****************************************************************************
\ No newline at end of file
diff --git a/src/topupopt/data/misc/utils.py b/src/topupopt/data/misc/utils.py
index 0077e519758cc5caab0fd9bb0119c1576f67e1cc..4536d1f27157b72d494589eddf06dff633e29c41 100644
--- a/src/topupopt/data/misc/utils.py
+++ b/src/topupopt/data/misc/utils.py
@@ -5,7 +5,6 @@
 import uuid
 import math
 from statistics import mean
-
 # local, external
 import pyomo.environ as pyo
 
@@ -38,7 +37,7 @@ def generate_pseudo_unique_key(key_list: tuple, max_iterations: int = 10) -> str
 def discrete_sinusoid_matching_integral(
     integration_result: float,
     time_interval_durations: list,
-    min_to_max_ratio: float,
+    min_max_ratio: float,
     phase_shift_radians: float = None,
 ) -> list:
     """
@@ -55,7 +54,7 @@ def discrete_sinusoid_matching_integral(
 
     where:
 
-    a = b*(1-min_to_max_ratio)/(1+min_to_max_ratio)
+    a = b*(1-min_max_ratio)/(1+min_max_ratio)
 
     b = integration_result/integration_period
 
@@ -69,7 +68,7 @@ def discrete_sinusoid_matching_integral(
         The result of integrating the sinusoidal function for one period.
     time_interval_durations : list
         The time interval durations for each sample.
-    min_to_max_ratio : float
+    min_max_ratio : float
         The ratio between the minimum and maximum values of the function.
     phase_shift_radians : float, optional
         The phase shift for the sinusoidal function. The default is None, which
@@ -88,7 +87,7 @@ def discrete_sinusoid_matching_integral(
 
     b = integration_result / integration_period
 
-    a = b * (1 - min_to_max_ratio) / (1 + min_to_max_ratio)
+    a = b * (1 - min_max_ratio) / (1 + min_max_ratio)
 
     alpha = 2 * math.pi / integration_period
 
@@ -186,10 +185,10 @@ def synch_profile(profile: list, reference_profile: list, synch: bool = True) ->
 
 def create_profile_using_time_weighted_state(
     integration_result: float,
-    avg_state: list,
+    states: list,
     time_interval_durations: list,
-    min_to_max_ratio: float,
-    state_correlates_with_output: bool = True,
+    min_max_ratio: float,
+    states_correlate_profile: bool = True,
 ) -> list:
     """
     Returns a profile that approximates a sinusoidal function in discrete time.
@@ -208,7 +207,7 @@ def create_profile_using_time_weighted_state(
 
     where:
 
-    a = b*(1-min_to_max_ratio)/(1+min_to_max_ratio)
+    a = b*(1-min_max_ratio)/(1+min_max_ratio)
 
     b = integration_result/integration_period
 
@@ -220,13 +219,13 @@ def create_profile_using_time_weighted_state(
     ----------
     integration_result : float
         The result of integrating the sinusoidal function for one period.
-    avg_state : list
+    states : list
         The average state during each time interval.
     time_interval_durations : list
         The time interval durations for each sample.
-    min_to_max_ratio : float
+    min_max_ratio : float
         The ratio between the minimum and maximum values of the function.
-    state_correlates_with_output : bool, optional
+    states_correlate_profile : bool, optional
         If True, the peak should happen when the state is at its highest point.
         If False, the peak should happen when the state is at its lowest point.
         The default is True.
@@ -244,26 +243,26 @@ def create_profile_using_time_weighted_state(
 
     """
 
-    if len(avg_state) != len(time_interval_durations):
+    if len(states) != len(time_interval_durations):
         raise ValueError("The inputs are inconsistent.")
 
     period = sum(time_interval_durations)
 
     avg_time_interval_duration = mean(time_interval_durations)
 
-    avg_state_weighted = [
+    states_weighted = [
         (
             x_k * delta_k / avg_time_interval_duration
-            if state_correlates_with_output
+            if states_correlate_profile
             else -x_k * delta_k / avg_time_interval_duration
         )
-        for delta_k, x_k in zip(time_interval_durations, avg_state)
+        for delta_k, x_k in zip(time_interval_durations, states)
     ]
 
     # find the peak
 
     _sorted = sorted(
-        ((state, index) for index, state in enumerate(avg_state_weighted)), reverse=True
+        ((state, index) for index, state in enumerate(states_weighted)), reverse=True
     )
 
     # create new list for time durations starting with that of the peak
@@ -278,7 +277,7 @@ def create_profile_using_time_weighted_state(
     new_profile = discrete_sinusoid_matching_integral(
         integration_result=integration_result,
         time_interval_durations=swapped_time_durations,
-        min_to_max_ratio=min_to_max_ratio,
+        min_max_ratio=min_max_ratio,
         phase_shift_radians=(
             math.pi / 2
             - 0.5 * (time_interval_durations[_sorted[0][1]] / period) * 2 * math.pi
@@ -506,6 +505,14 @@ def generate_state_correlated_profile(
     f = [dt_k/dt_mean for dt_k in time_interval_durations]
     set_I = tuple(i for i in range(number_time_steps))
     
+    if len(set(states)) == 1:
+        # alpha = 0, return trivial profile
+        return (
+            [fi*beta for fi in f], 
+            0,
+            beta
+            )
+    
     # create a dictionary with the data (using pyomo conventions)
     data_dict = {
         None: {
@@ -541,7 +548,7 @@ def max_min_sinusoidal_profile(
     integration_result: float or int,
     period: float or int,
     time_interval_duration: float or int,
-    min_to_max_ratio: float,
+    min_max_ratio: float,
 ) -> tuple:
     """
     Returns the maximum and minimum amount for a given time interval, according
@@ -558,7 +565,7 @@ def max_min_sinusoidal_profile(
 
     where:
 
-    a = b*(1-min_to_max_ratio)/(1+min_to_max_ratio)
+    a = b*(1-min_max_ratio)/(1+min_max_ratio)
 
     b = integration_result/integration_period
 
@@ -572,7 +579,7 @@ def max_min_sinusoidal_profile(
         The result of integrating the sinusoidal function for one period.
     time_interval_durations : list
         The time interval durations for each sample.
-    min_to_max_ratio : float
+    min_max_ratio : float
         The ratio between the minimum and maximum values of the function.
     phase_shift_radians : float, optional
         The phase shift for the sinusoidal function. The default is None, which
@@ -586,7 +593,7 @@ def max_min_sinusoidal_profile(
     """
 
     b = integration_result / period
-    a = b * (1 - min_to_max_ratio) / (1 + min_to_max_ratio)
+    a = b * (1 - min_max_ratio) / (1 + min_max_ratio)
     alpha = 2 * math.pi / period
     amplitude = a * (2 / alpha) * math.sin(alpha * time_interval_duration / 2)
 
@@ -595,6 +602,87 @@ def max_min_sinusoidal_profile(
         b * time_interval_duration - amplitude,
     )
 
+# *****************************************************************************
+# *****************************************************************************
+
+def generate_profile(
+        integration_result: float,
+        time_interval_durations: list,
+        **kwargs
+    ) -> tuple:
+    """
+    Returns a profile according to a variety of methods.
+
+    Parameters
+    ----------
+    integration_result : float
+        The value which must be obtained by adding up all samples.
+    time_interval_durations : list
+        A list with the durations of each time interval.
+    **kwargs : 
+        A sequence of key and value pairs for use in subsequent methods.
+
+    Returns
+    -------
+    numpy array
+        Returns the desired profile.
 
+    """
+    
+    # generate the profile
+    if 'states' not in kwargs: 
+        # min_max_ratio is necessary
+        # phase_shift_radians is optional
+        # states play no role
+        
+        # # remove unnecessary arguments
+        # if 'deviation_gain' in kwargs:
+        #     kwargs.pop('deviation_gain')
+        # if 'solver' in kwargs:
+        #     kwargs.pop('solver')
+        return discrete_sinusoid_matching_integral(
+            integration_result=integration_result,
+            time_interval_durations=time_interval_durations,
+            **kwargs
+            )
+    elif 'deviation_gain' not in kwargs:
+        # states matter but the gain must be determined
+        
+        if 'solver' in kwargs:
+            # - states_correlate_profile is necessary
+            # - min_max_ratio is necessary
+            # - solver is necessary
+            return generate_state_correlated_profile(
+                integration_result=integration_result, 
+                time_interval_durations=time_interval_durations,
+                **kwargs
+                )[0]
+        else:
+            # - states_correlate_profile is necessary
+            # - min_max_ratio is necessary
+            # - solver is not necessary
+            return create_profile_using_time_weighted_state(
+                integration_result=integration_result,
+                time_interval_durations=time_interval_durations,
+                **kwargs)
+    else:
+        # states matter and the gain is predefined
+        # states are necessary
+        # deviation gain is necessary
+        # # remove unnecessary arguments
+        # if 'phase_shift_radians' in kwargs:
+        #     kwargs.pop('phase_shift_radians')
+        return generate_manual_state_correlated_profile(
+            integration_result=integration_result, 
+            time_interval_durations=time_interval_durations,
+            **kwargs
+            )
+    
+    # integration_result: float,
+    # states: list,
+    # time_interval_durations: list,
+    # min_max_ratio: float,
+    # states_correlate_profile: bool = True,
+    
 # *****************************************************************************
 # *****************************************************************************
diff --git a/tests/test_data_buildings_dk.py b/tests/test_data_buildings_dk.py
index 5f0e69375c36f2f61629978294dc73345c212ab1..c19ba389d18ba4aad08666c500d44ea84bebfb54 100644
--- a/tests/test_data_buildings_dk.py
+++ b/tests/test_data_buildings_dk.py
@@ -3,9 +3,7 @@
 # standard
 
 import math
-import random
 from numbers import Real
-from statistics import mean
 
 import geopandas as gpd
 
@@ -15,7 +13,6 @@ import geopandas as gpd
 # local, internal
 from src.topupopt.data.gis.utils import read_gdf_file
 from src.topupopt.data.buildings.dk import heat
-from src.topupopt.data.buildings.dk import bbr
 
 # *****************************************************************************
 # *****************************************************************************
@@ -38,9 +35,11 @@ class TestDataBuildingsDK:
             30*24*3600
             for i in range(number_time_intervals)
             ]
-        annual_heat_demand_scenario = 1000
-        total_area = 1000
-        air_temperature_scenario = [10 for i in range(number_time_intervals)]
+        total_demand_true = 1000
+        total_area_true = 4563  # 5%: 4563 # 100%: 100882
+        assessments = ['q']
+        annual_heat_demand = {'q': 1000}
+        air_temperature =  {'q': [5+i for i in range(number_time_intervals)]}
         
         gdf_osm = gpd.read_file(osm_data_filename)
         gdf_osm.set_index(['element_type', 'osmid'], drop=True, inplace=True)
@@ -51,55 +50,259 @@ class TestDataBuildingsDK:
             index='index'
             )
         
-        # order by state
+        def verify_result(
+                out_dict, 
+                out_area, 
+                total_demand_true, 
+                total_area_true,
+                # assessments,
+                # number_time_intervals
+                ):
+            assert type(out_dict) == dict
+            assert isinstance(out_area, Real)
+            assert len(out_dict) == len(gdf_osm)
+            assert math.isclose(out_area, total_area_true, abs_tol=1e-3) # 5%: 4563 # 100%: 100882
+            for q in assessments:
+                assert math.isclose(
+                    sum(sum(v[q]) for k, v in out_dict.items() if len(v[q]) != 0),
+                    total_demand_true, 
+                    abs_tol=1e-3
+                    )
+            # output dict must be keyed by entrance id and then by scenario
+            for k, v in out_dict.items():
+                assert k in gdf_osm.index
+                if len(v) == 0:
+                    continue
+                for q in assessments:
+                    assert q in v
+                    assert len(v[q]) == number_time_intervals or len(v[q]) == 0
         
-        heat_demand_dict = heat.heat_demand_dict_by_building_entrance(
+        # drop entries to keep things fast
+        share_keeper_osm_entries = 0.05 
+        number_osm_entries = len(gdf_osm)
+        for index in gdf_osm.index:
+            if len(gdf_osm) < round(share_keeper_osm_entries*number_osm_entries):
+                break
+            gdf_osm.drop(index=index, inplace=True)
+        
+        # create profiles in accordance with a set of states and a positive gain
+        
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
+            gdf_osm=gdf_osm,
+            gdf_buildings=gdf_buildings,
+            time_interval_durations=intraperiod_time_interval_duration,
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            air_temperature=air_temperature,
+            deviation_gain=1
+            )
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+            
+        # create profiles in accordance with a set of states and a negative gain
+        
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
             gdf_osm=gdf_osm,
             gdf_buildings=gdf_buildings,
-            number_intervals=number_time_intervals,
             time_interval_durations=intraperiod_time_interval_duration,
-            bdg_ratio_min_max={
-                index: min_to_max_ratio for index in gdf_buildings.index
-                },
-            bdg_specific_demand={
-                index: annual_heat_demand_scenario/total_area 
-                for index in gdf_buildings.index
-                },
-            bdg_demand_phase_shift=None,
-            avg_state=air_temperature_scenario,
-            state_correlates_with_output=False
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            air_temperature=air_temperature,
+            deviation_gain=-1
             )
-        assert type(heat_demand_dict) == dict
-        assert len(heat_demand_dict) == len(gdf_osm)
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+         
+        # create profiles in accordance with a sinusoidal function (no phase shift)
         
-        # no state preference, use phase shift
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
+            gdf_osm=gdf_osm,
+            gdf_buildings=gdf_buildings,
+            time_interval_durations=intraperiod_time_interval_duration,
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            min_max_ratio=min_to_max_ratio,
+            # air_temperature=air_temperature,
+            # state_correlates_with_output=False
+            # deviation_gain=1
+            )
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+        
+        # create profiles in accordance with a sinusoidal function (with phase shift)
+        
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
+            gdf_osm=gdf_osm,
+            gdf_buildings=gdf_buildings,
+            time_interval_durations=intraperiod_time_interval_duration,
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            min_max_ratio=min_to_max_ratio,
+            phase_shift_radians=math.pi/2
+            # air_temperature=air_temperature,
+            # state_correlates_with_output=False
+            # deviation_gain=1
+            )
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+            
+        # create profiles in accordance with states but without a predefined gain
     
-        heat_demand_dict2 = heat.heat_demand_dict_by_building_entrance(
+        # create profile (no optimisation)
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
             gdf_osm=gdf_osm,
             gdf_buildings=gdf_buildings,
-            number_intervals=number_time_intervals,
             time_interval_durations=intraperiod_time_interval_duration,
-            bdg_ratio_min_max={
-                index: min_to_max_ratio for index in gdf_buildings.index
-                },
-            bdg_specific_demand={
-                index: annual_heat_demand_scenario/total_area 
-                for index in gdf_buildings.index
-                },
-            bdg_demand_phase_shift={
-                index: 2*math.pi*random.random() for index in gdf_buildings.index
-                },
-            avg_state=None,
-            state_correlates_with_output=False
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            air_temperature=air_temperature,
+            min_max_ratio=min_to_max_ratio,
+            states_correlate_profile=True,
             )
-        assert type(heat_demand_dict2) == dict
-        assert len(heat_demand_dict2) == len(gdf_osm)
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+            
+        # create profiles in accordance with states but without a predefined gain (optimisation)
+        
+        # remove all but one osm entry (to keep things light)
+        for index in gdf_osm.index:
+            if len(gdf_osm) <= 1:
+                break
+            gdf_osm.drop(index=index, inplace=True)
+        
+        # create profile
+        heat_demand_dict, total_area = heat.heat_demand_profiles(
+            gdf_osm=gdf_osm,
+            gdf_buildings=gdf_buildings,
+            time_interval_durations=intraperiod_time_interval_duration,
+            assessments=assessments,
+            annual_heat_demand=annual_heat_demand,
+            air_temperature=air_temperature,
+            min_max_ratio=min_to_max_ratio,
+            states_correlate_profile=True,
+            solver='glpk'
+            )
+        total_area_true = 200
+        verify_result(heat_demand_dict, total_area, total_demand_true, total_area_true)
+
+    # *************************************************************************
+    # *************************************************************************
+    
+    # def test_demand_dict3(self):
         
-        # total heating area
+    #     # heat_demand_dict_by_building_entrance
         
-        heating_area = heat.total_heating_area(gdf_osm, gdf_buildings)
-        assert isinstance(heating_area, Real)
-        assert math.isclose(heating_area, 100882, abs_tol=1e-3)
+    #     osm_data_filename = 'tests/data/gdf_osm.gpkg'
+    #     building_data_filename = 'tests/data/gdf_buildings.gpkg'
+    #     bdg_gdf_container_columns = ('ejerskaber','koordinater','bygningspunkt')
+    #     number_time_intervals = 12
+    #     min_to_max_ratio = 0.1
+    #     intraperiod_time_interval_duration = [
+    #         30*24*3600
+    #         for i in range(number_time_intervals)
+    #         ]
+    #     annual_heat_demand_scenario = 1000
+    #     total_area = 1000
+    #     states = [10 for i in range(number_time_intervals)]
+        
+    #     gdf_osm = gpd.read_file(osm_data_filename)
+    #     gdf_osm.set_index(['element_type', 'osmid'], drop=True, inplace=True)
+    
+    #     gdf_buildings = read_gdf_file(
+    #         filename=building_data_filename,
+    #         packed_columns=bdg_gdf_container_columns,
+    #         index='index'
+    #         )
+        
+    #     # sinusoidal
+        
+    #     heat_demand_dict = heat.heat_demand_dict_by_building_entrance2(
+    #         gdf_osm=gdf_osm,
+    #         gdf_buildings=gdf_buildings,
+    #         number_intervals=number_time_intervals,
+    #         time_interval_durations=intraperiod_time_interval_duration,
+    #         min_max_ratio=min_to_max_ratio,
+    #         specific_demand=annual_heat_demand_scenario/total_area,
+    #         )
+    #     assert type(heat_demand_dict) == dict
+    #     assert len(heat_demand_dict) == len(gdf_osm)
+    #     assert math.isclose(
+    #         annual_heat_demand_scenario, 
+    #         sum(sum(value) for value in heat_demand_dict.values()),
+    #         abs_tol=1e-3,
+    #         )
+        
+    #     # sinusoidal with phase shift
+    
+    #     heat_demand_dict = heat.heat_demand_dict_by_building_entrance2(
+    #         gdf_osm=gdf_osm,
+    #         gdf_buildings=gdf_buildings,
+    #         number_intervals=number_time_intervals,
+    #         time_interval_durations=intraperiod_time_interval_duration,
+    #         min_max_ratio=min_to_max_ratio,
+    #         specific_demand=annual_heat_demand_scenario/total_area ,
+    #         phase_shift_radians=math.pi,
+    #         )
+    #     assert type(heat_demand_dict) == dict
+    #     assert len(heat_demand_dict) == len(gdf_osm)
+    #     assert math.isclose(
+    #         annual_heat_demand_scenario, 
+    #         sum(sum(value) for value in heat_demand_dict.values()),
+    #         abs_tol=1e-3,
+    #         )
+        
+    #     # predefined deviation gain, positive
+    
+    #     heat_demand_dict = heat.heat_demand_dict_by_building_entrance2(
+    #         gdf_osm=gdf_osm,
+    #         gdf_buildings=gdf_buildings,
+    #         number_intervals=number_time_intervals,
+    #         time_interval_durations=intraperiod_time_interval_duration,
+    #         states=states,
+    #         specific_demand=annual_heat_demand_scenario/total_area ,
+    #         deviation_gain=3,
+    #         )
+    #     assert type(heat_demand_dict) == dict
+    #     assert len(heat_demand_dict) == len(gdf_osm)
+    #     assert math.isclose(
+    #         annual_heat_demand_scenario, 
+    #         sum(sum(value) for value in heat_demand_dict.values()),
+    #         abs_tol=1e-3,
+    #         )
+        
+    #     # predefined deviation gain, negative
+    
+    #     heat_demand_dict = heat.heat_demand_dict_by_building_entrance2(
+    #         gdf_osm=gdf_osm,
+    #         gdf_buildings=gdf_buildings,
+    #         number_intervals=number_time_intervals,
+    #         time_interval_durations=intraperiod_time_interval_duration,
+    #         states=states,
+    #         specific_demand=annual_heat_demand_scenario/total_area ,
+    #         deviation_gain=-3,
+    #         )
+    #     assert type(heat_demand_dict) == dict
+    #     assert len(heat_demand_dict) == len(gdf_osm)
+    #     assert math.isclose(
+    #         annual_heat_demand_scenario, 
+    #         sum(sum(value) for value in heat_demand_dict.values()),
+    #         abs_tol=1e-3,
+    #         )
+        
+    #     # optimisation
+    
+    #     heat_demand_dict = heat.heat_demand_dict_by_building_entrance2(
+    #         gdf_osm=gdf_osm,
+    #         gdf_buildings=gdf_buildings,
+    #         number_intervals=number_time_intervals,
+    #         time_interval_durations=intraperiod_time_interval_duration,
+    #         states=states,
+    #         specific_demand=annual_heat_demand_scenario/total_area,
+    #         states_correlate_profile=True,
+    #         solver='glpk'
+    #         )
+    #     assert type(heat_demand_dict) == dict
+    #     assert len(heat_demand_dict) == len(gdf_osm)
+    #     assert math.isclose(
+    #         annual_heat_demand_scenario, 
+    #         sum(sum(value) for value in heat_demand_dict.values()),
+    #         abs_tol=1e-3,
+    #         )
     
     # *************************************************************************
     # *************************************************************************
diff --git a/tests/test_data_utils.py b/tests/test_data_utils.py
index 32aebc674834939cc590726aa1d7f5ffbd038f65..d653fe46c327ec833a8bae0b3f50ceab8b7cd2a2 100644
--- a/tests/test_data_utils.py
+++ b/tests/test_data_utils.py
@@ -13,9 +13,9 @@ class TestDataUtils:
     def test_profile_synching2(self):
         integration_result = 10446
         ratio_min_avg = 0.2
-        min_to_max_ratio = ratio_min_avg / (2 - ratio_min_avg)
+        min_max_ratio = ratio_min_avg / (2 - ratio_min_avg)
 
-        avg_state = [
+        states = [
             2.66,
             2.34,
             3.54,
@@ -52,10 +52,10 @@ class TestDataUtils:
 
         new_profile = utils.create_profile_using_time_weighted_state(
             integration_result=integration_result,
-            avg_state=avg_state,
+            states=states,
             time_interval_durations=time_interval_durations,
-            min_to_max_ratio=min_to_max_ratio,
-            state_correlates_with_output=False,
+            min_max_ratio=min_max_ratio,
+            states_correlate_profile=False,
         )
 
         expected_result = [
@@ -89,10 +89,10 @@ class TestDataUtils:
 
         new_profile = utils.create_profile_using_time_weighted_state(
             integration_result=integration_result,
-            avg_state=avg_state,
+            states=states,
             time_interval_durations=time_interval_durations,
-            min_to_max_ratio=min_to_max_ratio,
-            state_correlates_with_output=True,
+            min_max_ratio=min_max_ratio,
+            states_correlate_profile=True,
         )
 
         expected_result = [
@@ -126,7 +126,7 @@ class TestDataUtils:
             integration_result=integration_result,
             period=sum(time_interval_durations),
             time_interval_duration=mean(time_interval_durations),
-            min_to_max_ratio=min_to_max_ratio,
+            min_max_ratio=min_max_ratio,
         )
 
         expected_pmax, expected_pmin = 1558.972133279683, 182.02786672031687
@@ -145,10 +145,10 @@ class TestDataUtils:
         try:
             new_profile = utils.create_profile_using_time_weighted_state(
                 integration_result=integration_result,
-                avg_state=avg_state,
+                states=states,
                 time_interval_durations=time_interval_durations,
-                min_to_max_ratio=min_to_max_ratio,
-                state_correlates_with_output=True,
+                min_max_ratio=min_max_ratio,
+                states_correlate_profile=True,
             )
         except ValueError:
             error_triggered = True
@@ -209,12 +209,12 @@ class TestDataUtils:
 
             integration_result = 100
 
-            min_to_max_ratio = 0.2
+            min_max_ratio = 0.2
 
             profile = utils.discrete_sinusoid_matching_integral(
                 integration_result,
                 time_interval_durations,
-                min_to_max_ratio,
+                min_max_ratio,
                 phase_shift_radians=phase_shift_radians,
             )
 
@@ -257,12 +257,12 @@ class TestDataUtils:
 
             integration_result = 100
 
-            min_to_max_ratio = 0.2
+            min_max_ratio = 0.2
 
             profile = utils.discrete_sinusoid_matching_integral(
                 integration_result,
                 time_interval_durations,
-                min_to_max_ratio,
+                min_max_ratio,
                 phase_shift_radians=phase_shift_radians,
             )
 
@@ -305,10 +305,10 @@ class TestDataUtils:
 
         integration_result = 100
 
-        min_to_max_ratio = 0.2
+        min_max_ratio = 0.2
 
         profile = utils.discrete_sinusoid_matching_integral(
-            integration_result, time_interval_durations, min_to_max_ratio
+            integration_result, time_interval_durations, min_max_ratio
         )
 
         assert math.isclose(sum(profile), integration_result, abs_tol=0.01)
@@ -521,6 +521,190 @@ class TestDataUtils:
 
     # *************************************************************************
     # *************************************************************************
+    
+    def test_create_profile_sinusoidal(self):
+        
+        number_intervals = 10
+        integration_result = 100
+        min_max_ratio = 0.25
+        
+        # sinusoidal profile
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            min_max_ratio=min_max_ratio, 
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)
+        
+        # sinusoidal profile with phase shift
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            min_max_ratio=min_max_ratio, 
+            phase_shift_radians=math.pi/2
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)
+        
+        # use incorrect parameter
+        
+        error_raised = False
+        try:
+            profile = utils.generate_profile(
+                integration_result=integration_result, 
+                time_interval_durations=[1 for i in range(number_intervals)], 
+                min_max_ratio=min_max_ratio, 
+                deviation_gain=-1,
+                )
+        except TypeError:
+            error_raised = True
+        assert error_raised
+
+    # *************************************************************************
+    # *************************************************************************
+        
+    def test_create_profile_predefined_gain(self):
+        
+        number_intervals = 10
+        integration_result = 100
+        deviation_gain = 5
+        states = [number_intervals-i*0.5 for i in range(number_intervals)]
+        
+        # predefined gain
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            states=states, 
+            deviation_gain=deviation_gain
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)      
+        
+        # predefined gain, opposite sign
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            states=states, 
+            deviation_gain=-deviation_gain
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)    
+        
+        # use incorrect parameter
+        
+        error_raised = False
+        try:
+            profile = utils.generate_profile(
+                integration_result=integration_result, 
+                time_interval_durations=[1 for i in range(number_intervals)], 
+                states=states, 
+                deviation_gain=-deviation_gain,
+                phase_shift_radians=math.pi
+                )
+        except TypeError:
+            error_raised = True
+        assert error_raised
+
+    # *************************************************************************
+    # *************************************************************************
+        
+    def test_create_profile_via_sorting_sinusoid(self):
+        
+        number_intervals = 10
+        integration_result = 100
+        states_correlate_profile = True
+        min_max_ratio = 0.25
+        states = [number_intervals-i*0.5 for i in range(number_intervals)]
+        
+        # sorting and sinusoidal function        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            min_max_ratio=min_max_ratio, 
+            states=states, 
+            states_correlate_profile=states_correlate_profile, 
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)    
+        
+
+    # *************************************************************************
+    # *************************************************************************
+        
+    def test_create_profile_via_optimisation(self):
+        
+        number_intervals = 10
+        integration_result = 100
+        states_correlate_profile = True
+        min_max_ratio = 0.25
+        solver = 'glpk'
+        states = [number_intervals-i*0.5 for i in range(number_intervals)]
+        
+        # optimisation
+        
+        # states_correlate_profile is necessary
+        # min_max_ratio is necessary
+        # solver is necessary
+        # states matter but the gain must be determined
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            min_max_ratio=min_max_ratio, 
+            states=states, 
+            states_correlate_profile=states_correlate_profile, 
+            solver=solver
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)    
+        assert math.isclose(min(profile),max(profile)*min_max_ratio, abs_tol=1e-3)
+        
+        # optimisation but with states that do no warrant it
+        states = [5 for i in range(number_intervals)]
+        
+        profile = utils.generate_profile(
+            integration_result=integration_result, 
+            time_interval_durations=[1 for i in range(number_intervals)], 
+            min_max_ratio=min_max_ratio, 
+            states=states, 
+            states_correlate_profile=states_correlate_profile, 
+            solver=solver
+            )
+        
+        assert len(profile) == number_intervals
+        assert math.isclose(sum(profile), integration_result, abs_tol=1e-3)   
+        # the min to max ratio cannot be observed if the states do not change
+        assert math.isclose(min(profile), max(profile), abs_tol=1e-3)
+        
+        # use incorrect parameter
+        error_raised = False
+        try:
+            profile = utils.generate_profile(
+                integration_result=integration_result, 
+                time_interval_durations=[1 for i in range(number_intervals)], 
+                min_max_ratio=min_max_ratio, 
+                states=states, 
+                states_correlate_profile=states_correlate_profile, 
+                solver=solver,
+                phase_shift_radians=math.pi
+                )
+        except TypeError:
+            error_raised = True
+        assert error_raised
+        
+    # *************************************************************************
+    # *************************************************************************
 
 # *****************************************************************************
 # *****************************************************************************
\ No newline at end of file