diff --git a/src/__init__.py b/src/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/__init__.py b/src/materials/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py b/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py
new file mode 100644
index 0000000000000000000000000000000000000000..ff5897e704531a1e1244bd0c9053ccfa27b06100
--- /dev/null
+++ b/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py
@@ -0,0 +1,103 @@
+import numpy as np
+import sympy as sp
+import pystencils as ps
+from dataclasses import dataclass
+from pystencils.typing import PointerType
+from src.materials.constants import Constants
+from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
+
+
+# Material Data Class
+@dataclass
+class Ti6Al4V:                                          # 90% Ti, 6% Al, 4% V
+    latent_heat =                           290000      # m²/s²
+    temperature_solidus_1 =                 1876.0      # k         Ref [1]
+    temperature_solidus_2 =                 1880.0      # k         Ref [2]
+    temperature_liquidus =                  1928.0      # k         Ref [2]
+    density_room =                          4430.0      # kg/m³     Ref [5]
+    thermal_expansion_coefficient =         10e-6       # /k        Ref[6]
+    heat_capacity_solidus =                 700.0       # J/kg-k
+    heat_capacity_liquidus =                1126.0      # J/kg-k
+    heat_conductivity_solidus =             21.58       # W/m-k
+    heat_conductivity_liquidus =            33.60       # W/m-k
+    density_temp_array = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+    density_data_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+    density_temp_base =                     1878.0
+    density_temp_incr =                     8.0
+    thermal_diffusivity_solidus =           727.7
+    thermal_diffusivity_liquidus =          708.6
+    # electrical_resistivity = [1, 2, 5, 10, 20, 50, 100]
+    # electrical_resistivity_temperatures = [800, 900, 1200, 1500, 1700, 1900, 2000]
+    # Define temperature and density arrays
+    density_temp_array_lookup = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
+    density_data_array_lookup = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
+
+
+'''
+@dataclass
+class Constants:
+    temperature_room =                      298.15          # k
+    N_a =                                   6.022141e23     # /mol
+    u =                                     1.660538e-27    # kg
+'''
+
+# Instantiate the Ti6Al4V class
+ti6al4v = Ti6Al4V()
+
+Ti6Al4V_latent_heat = ti6al4v.latent_heat
+Ti6Al4V_thermal_expansion_coefficient = ti6al4v.thermal_expansion_coefficient
+
+'''
+@dataclass
+class Solidus_Liquidus_Parameter:
+    sol: float
+    liq: float
+
+
+@dataclass
+class Functional_Parameter:
+    func: str
+    args: list
+'''
+
+# Parameters for Ti6Al4V
+Ti6Al4V_temperatures = Solidus_Liquidus_Parameter(
+    sol=0.5 * (ti6al4v.temperature_solidus_1 + ti6al4v.temperature_solidus_2),
+    liq=ti6al4v.temperature_liquidus)
+
+Ti6Al4V_heat_conductivity = Solidus_Liquidus_Parameter(
+    sol=ti6al4v.heat_conductivity_solidus,
+    liq=ti6al4v.heat_conductivity_liquidus)
+
+Ti6Al4V_specific_heat_capacity = Solidus_Liquidus_Parameter(
+    sol=ti6al4v.heat_capacity_solidus,
+    liq=ti6al4v.heat_capacity_liquidus)
+
+Ti6Al4V_density_tec = Functional_Parameter(
+    func='thermal_expansion',
+    args=[Constants.temperature_room, ti6al4v.density_room, ti6al4v.thermal_expansion_coefficient])
+
+Ti6Al4V_density_int = Functional_Parameter(
+    func='interpolation',
+    args=[ti6al4v.density_temp_array, ti6al4v.density_data_array])
+
+Ti6Al4V_density_in_array = Functional_Parameter(
+    func='interpolation_in_array',  # TODO
+    args=[ti6al4v.density_temp_base, ti6al4v.density_temp_incr, ti6al4v.density_data_array,
+          sp.IndexedBase(ps.TypedSymbol('density_data_symbol', PointerType(sp.Basic('float32'))),
+                         shape=len(ti6al4v.density_data_array))])
+# PointerType(create_type(np.float32)) PointerType(base_type=float32) PointerType(sp.Basic('float32'))
+
+Ti6Al4V_density_lookup = Functional_Parameter(
+    func='lookup',
+    args=[ti6al4v.density_temp_array_lookup, ti6al4v.density_data_array_lookup])
+
+Ti6Al4V_thermal_diffusivity1 = Solidus_Liquidus_Parameter(
+    sol=ti6al4v.thermal_diffusivity_solidus,
+    liq=ti6al4v.thermal_diffusivity_liquidus
+)
+
+Ti6Al4V_thermal_diffusivity = Functional_Parameter(
+    func='thermal_diffusivity',
+    args=[Ti6Al4V_heat_conductivity, Ti6Al4V_density_int, Ti6Al4V_specific_heat_capacity]
+)
diff --git a/src/materials/alloys/Ti6Al4V/__init__.py b/src/materials/alloys/Ti6Al4V/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/alloys/__init__.py b/src/materials/alloys/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/constants.py b/src/materials/constants.py
new file mode 100644
index 0000000000000000000000000000000000000000..9c804c57558a058feee8e95d2923c5bad65c8131
--- /dev/null
+++ b/src/materials/constants.py
@@ -0,0 +1,8 @@
+from dataclasses import dataclass
+
+
+@dataclass
+class Constants:
+    temperature_room =                      298.15          # k
+    N_a =                                   6.022141e23     # /mol
+    u =                                     1.660538e-27    # kg
\ No newline at end of file
diff --git a/src/materials/dataclasses.py b/src/materials/dataclasses.py
new file mode 100644
index 0000000000000000000000000000000000000000..a67834b2d70628044a2da0986dbed7731c23b4a2
--- /dev/null
+++ b/src/materials/dataclasses.py
@@ -0,0 +1,13 @@
+from dataclasses import dataclass
+
+
+@dataclass
+class Solidus_Liquidus_Parameter:
+    sol: float
+    liq: float
+
+
+@dataclass
+class Functional_Parameter:
+    func: str
+    args: list
diff --git a/src/materials/elements/Ti.py b/src/materials/elements/Ti.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/elements/__init__.py b/src/materials/elements/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/materials/get_arguments.py b/src/materials/get_arguments.py
new file mode 100644
index 0000000000000000000000000000000000000000..d6a3bc2e0aca0c3be9ba58ac4b9d87ef2dfcc602
--- /dev/null
+++ b/src/materials/get_arguments.py
@@ -0,0 +1,20 @@
+# from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
+from get_parameters import get_parameters
+from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
+
+
+def get_arguments(T, temperatures, args):
+    """
+    Recursively get arguments for parameter functions.
+
+    :param T: Temperature value.
+    :param temperatures: Solidus and liquidus temperatures.
+    :param args: Arguments to be processed.
+    """
+    converted_args = []
+    for arg in args:
+        if isinstance(arg, (Solidus_Liquidus_Parameter, Functional_Parameter)):
+            converted_args.append(get_parameters(T, temperatures, arg))
+        else:
+            converted_args.append(arg)
+    return converted_args
diff --git a/src/materials/get_parameters.py b/src/materials/get_parameters.py
new file mode 100644
index 0000000000000000000000000000000000000000..311ac0eda1cf9470340166e322c26aa3be0a3868
--- /dev/null
+++ b/src/materials/get_parameters.py
@@ -0,0 +1,44 @@
+import sympy as sp
+from src.materials.interpolators import interpolate, interpolate_property, density_lookup
+from src.materials.models import density_by_thermal_expansion, thermal_diffusivity_func
+# from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
+from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
+
+
+def get_parameters(T, temperatures, parameter):
+    """
+    Resolve parameters based on temperature and parameter type.
+
+    :param T: Temperature value.
+    :param temperatures: Solidus and liquidus temperatures.
+    :param parameter: Parameter to be resolved.
+    """
+    if isinstance(parameter, Solidus_Liquidus_Parameter):
+        return sp.Piecewise(
+            (parameter.sol, T <= temperatures.sol),
+            (parameter.liq, T >= temperatures.liq),
+            (interpolate(T, temperatures.sol, temperatures.liq, parameter.sol, parameter.liq), True)
+        )
+    elif isinstance(parameter, Functional_Parameter):
+        resolved_args = [get_parameters(T, temperatures, arg) for arg in parameter.args]
+        if parameter.func == 'thermal_expansion':
+            return density_by_thermal_expansion(T, *resolved_args)  # *get_parameters(T, temperatures, parameter.args)
+        elif parameter.func == 'interpolation':
+            return interpolate_property(T, *resolved_args)
+        elif parameter.func == 'interpolation_in_array':  # TODO
+            T_base, T_incr, data_array, data_symbol = resolved_args
+            idx = sp.floor((T - T_base) / T_incr)
+            mod = (T - T_base) - idx * T_incr  # mod replaces (T - T_base) % T_incr
+            return sp.Piecewise(
+                (data_array[0], T < T_base),
+                (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
+                (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True)
+            )
+        elif parameter.func == 'lookup':
+            return density_lookup(T, *resolved_args)
+        elif parameter.func == 'thermal_diffusivity':
+            return thermal_diffusivity_func(*resolved_args)
+        else:
+            raise ValueError("No known format in the dict: parameter = ", parameter)
+    else:
+        return parameter  # assume everything which is not a dict is a float, int, etc.
diff --git a/src/materials/interpolators.py b/src/materials/interpolators.py
new file mode 100644
index 0000000000000000000000000000000000000000..242254280a2e6bf712862ac077744f67c22fda12
--- /dev/null
+++ b/src/materials/interpolators.py
@@ -0,0 +1,114 @@
+import numpy as np
+import sympy as sp
+
+
+def interpolate(x, xmin, xmax, ymin, ymax):
+    """
+    Linear interpolation between two points: y = y_1 + ((y2 - y1) / (x2 - x1)) * (x - x1)
+    """
+    return (ymin * (xmax - x) + ymax * (x - xmin)) / (xmax - xmin)
+
+
+def interpolate_property(T, base_temperature, base_values):  # TODO
+    """
+    Interpolates properties based on temperature for both symbolic and numeric cases.
+
+    :param T: Temperature value, either symbolic or numeric.
+    :param base_temperature: List or array of base temperatures.
+    :param base_values: List or array of property values corresponding to base temperatures.
+    """
+    # Ensure base_temperature and base_values are lists or arrays
+    if not isinstance(base_temperature, (list, np.ndarray)) or not isinstance(base_values, (list, np.ndarray)):
+        raise TypeError("base_temperature and base_values must be of type list or np.ndarray")
+    # Convert to lists if they are numpy arrays
+    if isinstance(base_temperature, np.ndarray):
+        base_temperature = base_temperature.tolist()
+    if isinstance(base_values, np.ndarray):
+        base_values = base_values.tolist()
+
+    # Debugging: Print lengths and types
+    # print(f"Length of base_temperature: {len(base_temperature)}")
+    # print(f"Length of base_values: {len(base_values)}")
+    # print(f"Type of base_temperature: {type(base_temperature)}")
+    # print(f"Type of base_values: {type(base_values)}")
+
+    # Check if base_temperature and base_values have the same length
+    if len(base_temperature) != len(base_values):
+        raise ValueError("base_temperature and base_values must have the same length")
+    # assert len(base_temperature) == len(base_values), "base_temperature and base_values must have the same length"
+
+    if isinstance(T, sp.Expr):  # Symbolic expression
+        base_temperature = [sp.Float(t) for t in base_temperature]
+        base_values = [sp.Float(val) for val in base_values]
+        conditions = [
+            (base_values[0], T < base_temperature[0]),
+            (base_values[-1], T >= base_temperature[-1])]
+        for i in range(len(base_temperature) - 1):
+            interp_expr = (base_values[i] +
+                           (base_values[i+1] - base_values[i]) / (base_temperature[i+1] - base_temperature[i]) *
+                           (T - base_temperature[i]))
+            conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
+        return sp.Piecewise(*conditions)
+    else:  # Numeric expression
+        return np.interp(T, base_temperature, base_values)
+
+
+# Example usage:
+Temp = sp.Symbol('T')
+density_temp_array1 = [200.0, 300.0, 400.0, 500.0, 600.0]
+density_data_array1 = [50.0, 40.0, 30.0, 20.0, 10.0]
+density_temp_array2 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+density_data_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+
+# Symbolic interpolation
+symbolic_result = interpolate_property(Temp, density_temp_array2, density_data_array2)
+print("Symbolic interpolation result:", symbolic_result)
+
+# Numeric interpolation
+numeric_result = interpolate_property(450.0, density_temp_array2, density_data_array2)
+print("Numeric interpolation result:", numeric_result)
+
+
+def density_lookup(T, temperature_array, density_array):
+    """
+    Look up density based on temperature using interpolation.
+
+    :param T: Temperature value.
+    :param temperature_array: Array of temperatures.
+    :param density_array: Array of densities.
+    """
+    if len(temperature_array) != len(density_array):
+        raise ValueError("temperature_array and density_array must have the same length")
+    if isinstance(T, sp.Expr):  # Symbolic expression
+        conditions = [
+            (density_array[0], T < temperature_array[0]),
+            (density_array[-1], T >= temperature_array[-1])
+        ]
+        for i in range(len(temperature_array) - 1):
+            interp_expr = (density_array[i] +
+                           (density_array[i+1] - density_array[i]) / (temperature_array[i+1] - temperature_array[i]) *
+                           (T - temperature_array[i]))
+            conditions.append((interp_expr, sp.And(T >= temperature_array[i], T < temperature_array[i + 1])))
+        return sp.Piecewise(*conditions)
+    else:  # Numeric expression
+        return np.interp(T, temperature_array, density_array)
+
+
+# Example usage:
+Temp = sp.Symbol('T')
+density_temp_array3 = [200.0, 300.0, 400.0, 500.0, 600.0]
+density_data_array3 = [50.0, 40.0, 30.0, 20.0, 10.0]
+density_temp_array4 = np.array([200.0, 300.0, 400.0, 500.0, 600.0])
+density_data_array4 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
+density_temp_array5 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+density_data_array5 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+density_temp_array6 = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
+density_data_array6 = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
+
+# Symbolic interpolation
+symbolic_result = density_lookup(Temp, density_temp_array4, density_data_array4)
+print("Symbolic interpolation result:", symbolic_result)
+
+# Numeric interpolation
+numeric_result = density_lookup(450.0, density_temp_array4, density_data_array4)
+print("Numeric interpolation result:", numeric_result)
diff --git a/src/materials/models.py b/src/materials/models.py
new file mode 100644
index 0000000000000000000000000000000000000000..cd098fec05fdc4899f95e48450874ce3386d1a8f
--- /dev/null
+++ b/src/materials/models.py
@@ -0,0 +1,21 @@
+def density_by_thermal_expansion(T, base_temperature, base_density, thermal_expansion_coefficient):
+    """
+    Calculate density based on thermal expansion: rho_0 / (1 + tec * (T - T_0))
+
+    :param T: Temperature value.
+    :param base_temperature: Reference temperature.
+    :param base_density: Density at reference temperature.
+    :param thermal_expansion_coefficient: Thermal expansion coefficient.
+    """
+    return base_density / (1 + thermal_expansion_coefficient * (T - base_temperature)) ** 3
+
+
+def thermal_diffusivity_func(heat_conductivity, density, specific_heat_capacity):
+    """
+    Calculate thermal diffusivity: k(T) / (rho(T) * c_p(T))
+
+    :param heat_conductivity: Thermal conductivity.
+    :param density: Density.
+    :param specific_heat_capacity: Specific heat capacity.
+    """
+    return heat_conductivity / (density * specific_heat_capacity)