Repository: MassimoCimmino/pygfunction (language: python)

Issue:
Inclined boreholes 
Does the most current version account for inclined boreholes? 

Given that the segment-to-segment response factor is calculated by equation 8 (Cimmino and Bernier 2014), how would the distance formula `d_{ij}` be implemented?
![equation_8](https://user-images.githubusercontent.com/39248734/88466475-36e8e080-ce92-11ea-990c-167d8a5ffe7f.png)

For example, currently calculation of `d_{ij}` is given by the equation 11.
![equation_11](https://user-images.githubusercontent.com/39248734/88466551-2edd7080-ce93-11ea-9d68-22d31b5b760b.png)

Would equation 11 need to become `d_{uv}?` Given the implementation of inclined boreholes, could it now be that `x_i == x_u is false`, and `x_j == x_v is false`?

Additionally, the location of segment u,  `(x_u, y_u)` could vary slightly along the segment. Would the method to calculate the distance `d_{uv}` result in the distance from top-to-top, center-to-center, bottom-to-bottom or some other approach?

I'm assuming this implementation would occur in  <a href="https://github.com/MassimoCimmino/pygfunction/blob/982ab6c47dae4931415f1a0f9a43bbfcfb34580c/pygfunction/gfunction.py#L709">_borehole_segments</a>. There has been extensive validation of the analytical implementation of the uniform borehole wall temperature (UBWHT) boundary condition when compared to Eskilson's numerical method for vertical boreholes. Given implementation, inclined boreholes would need compared to some of Eskilson's inclined borehole g-functions. 

---------------------------------------------------------------------------------------------------------------------------------------------------------------- 
Cimmino, Massimo, and Michel Bernier. "A Semi-analytical Method to Generate G-functions for Geothermal Bore Fields." International Journal of Heat and Mass Transfer 70 (2014): 641-50. Web.



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c6b0f887..239fb7b6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### New features
 
+* [Issue 50](https://github.com/MassimoCimmino/pygfunction/issues/50) - Implemented inclined boreholes for the evaluation of *g*-functions. The implementation includes an approximation of the FLS solution for inclined boreholes based on the method of Cimmino (2021) (see [Issue 138](https://github.com/MassimoCimmino/pygfunction/issues/138)). The `'equivalent'` solver is not yet supported.
 * [Issue 138](https://github.com/MassimoCimmino/pygfunction/issues/138) - Implemented the approximation of the finite line source solution of Cimmino (2021). The approximation avoids the numerical evaluation of integrals. This speeds up the calculation of g-functions when enabled.
 * [Issue 148](https://github.com/MassimoCimmino/pygfunction/issues/148) - Implemented `effective_borehole_thermal_resistance()` and `local_borehole_thermal_resistance()` methods for all pipe classes. Deprecated `pipes.borehole_thermal_resistance()`, which computed the effective borehole thermal resistance. It will be removed in `v3.0.0`. Implemented a new `update_thermal_resistances()` method to all pipe classes. This method allows to update the delta-circuit of thermal resistance of the boreholes based on provided values for the fluid thermal resistances. This allows simulations with time-variable fluid thermal resistances.
 
diff --git a/doc/source/example_index.rst b/doc/source/example_index.rst
index b77ffec0..17e869c3 100644
--- a/doc/source/example_index.rst
+++ b/doc/source/example_index.rst
@@ -25,3 +25,4 @@ Examples
    examples/multipole_temperature
    examples/bore_field_thermal_resistance
    examples/mixed_inlet_conditions
+   examples/inclined_boreholes
diff --git a/doc/source/examples/inclined_boreholes.rst b/doc/source/examples/inclined_boreholes.rst
new file mode 100644
index 00000000..77b16b0e
--- /dev/null
+++ b/doc/source/examples/inclined_boreholes.rst
@@ -0,0 +1,27 @@
+.. examples:
+
+**************************************************
+Calculation of g-functions with inclined boreholes
+**************************************************
+
+This example demonstrates the use of the :doc:`g-function <../modules/gfunction>` module
+to calculate *g*-functions of fields of inclined boreholes using a boundary condition of
+uniform and equal borehole wall temperature for all boreholes. The total rate of heat
+extraction in the bore field is constant.
+
+The following script generates the *g*-functions of two bore fields. The first field
+corresponds to the optimal configuration presented by Claesson and Eskilson [1]_. The
+second field corresponds to the configuration comprised of 8 boreholes in a circle
+presented by the same authors.
+
+The script is located in:
+`pygfunction/examples/inclined_boreholes.py`
+
+.. literalinclude:: ../../../examples/inclined_boreholes.py
+   :language: python
+   :linenos:
+
+.. rubric:: References
+.. [1] Claesson J, and Eskilson, P. (1987). Conductive heat extraction by
+   thermally interacting deep boreholes, in "Thermal analysis of heat
+   extraction boreholes". Ph.D. Thesis, University of Lund, Lund, Sweden. 
diff --git a/examples/inclined_boreholes.py b/examples/inclined_boreholes.py
new file mode 100644
index 00000000..40020724
--- /dev/null
+++ b/examples/inclined_boreholes.py
@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+""" Example of calculation of g-functions with inclined boreholesusing uniform
+    and equal borehole wall temperatures.
+
+    The g-functions of two fields of 8 boreholes are calculated for boundary
+    condition of uniform borehole wall temperature along the boreholes, equal
+    for all boreholes. The first field corresponds to the "optimum"
+    configuration presented by Claesson and Eskilson (1987) and the second
+    field corresponds to the configuration comprised of 8 boreholes in a
+    circle.
+    
+    Claesson J, and Eskilson, P. (1987). Conductive heat extraction by
+    thermally interacting deep boreholes, in "Thermal analysis of heat
+    extraction boreholes". Ph.D. Thesis, University of Lund, Lund, Sweden. 
+"""
+
+import pygfunction as gt
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def main():
+    # -------------------------------------------------------------------------
+    # Simulation parameters
+    # -------------------------------------------------------------------------
+
+    # Borehole dimensions
+    D = 4.0             # Borehole buried depth (m)
+    H = 150.0           # Borehole length (m)
+    r_b = 0.075         # Borehole radius (m)
+    tilt = np.radians(20.)  # Borehole inclination (rad)
+
+    # Thermal properties
+    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
+
+    # g-Function calculation options
+    options = {'disp': True}
+
+    # Geometrically expanding time vector.
+    dt = 100*3600.                  # Time step
+    tmax = 3000. * 8760. * 3600.    # Maximum time
+    Nt = 25                         # Number of time steps
+    ts = H**2/(9.*alpha)            # Bore field characteristic time
+    time = gt.utilities.time_geometric(dt, tmax, Nt)
+    lntts = np.log(time/ts)
+
+    # -------------------------------------------------------------------------
+    # Borehole fields
+    # -------------------------------------------------------------------------
+    """
+    Bore field #1
+
+    This field corresponds to the optimal configuration presented by
+    Claesson and Eskilson (1987). The field is built using the `cardinal_point`
+    function to define the orientation of each borehole, individually.
+    """
+    B = 7.5 # Borehole spacing (m)
+    # Orientation of the boreholes
+    borehole_orientations = [
+        gt.utilities.cardinal_point('W'),
+        gt.utilities.cardinal_point('NW'),
+        gt.utilities.cardinal_point('SW'),
+        gt.utilities.cardinal_point('N'),
+        gt.utilities.cardinal_point('S'),
+        gt.utilities.cardinal_point('NE'),
+        gt.utilities.cardinal_point('SE'),
+        gt.utilities.cardinal_point('E')]
+
+    # "Optimal" field of 8 boreholes
+    boreField1 = []
+    for i, orientation in enumerate(borehole_orientations):
+        borehole = gt.boreholes.Borehole(
+            H, D, r_b, i * B, 0., tilt=tilt, orientation=orientation)
+        boreField1.append(borehole)
+
+    # Visualize the borehole field
+    fig1 = gt.boreholes.visualize_field(boreField1)
+
+    """
+    Bore field #2
+
+    This field corresponds to the configuration comprised of 8 boreholes in
+    a circle presented by Claesson and Eskilson (1987). The field is built
+    using the `circle_field` function.
+    """
+    N = 8   # Number of boreholes
+    R = 3.  # Borehole spacing from the center of the field (m)
+
+    # Field of 6 boreholes in a circle
+    boreField2 = gt.boreholes.circle_field(N, R, H, D, r_b, tilt=tilt)
+
+    # Visualize the borehole field
+    fig2 = gt.boreholes.visualize_field(boreField2)
+
+    # -------------------------------------------------------------------------
+    # Evaluate g-functions for all fields
+    # -------------------------------------------------------------------------
+    # Bore field #1
+    gfunc1 = gt.gfunction.gFunction(
+        boreField1, alpha, time=time, options=options, method='similarities')
+    fig3 = gfunc1.visualize_g_function()
+    fig3.suptitle('"Optimal" field of 8 boreholes')
+    fig3.tight_layout()
+    # Bore field #2
+    gfunc2 = gt.gfunction.gFunction(
+        boreField2, alpha, time=time, options=options, method='similarities')
+    fig4 = gfunc2.visualize_g_function()
+    fig4.suptitle(f'Field of {N} boreholes in a circle')
+    fig4.tight_layout()
+
+
+# Main function
+if __name__ == '__main__':
+    main()
diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py
index d372d739..7befb901 100644
--- a/pygfunction/boreholes.py
+++ b/pygfunction/boreholes.py
@@ -35,8 +35,12 @@ def __init__(self, H, D, r_b, x, y, tilt=0., orientation=0.):
         self.r_b = float(r_b)  # Borehole radius
         self.x = float(x)      # Borehole x coordinate position
         self.y = float(y)      # Borehole y coordinate position
+        # Borehole inclination
         self.tilt = float(tilt)
+        # Borehole orientation
         self.orientation = float(orientation)
+        # Check if borehole is inclined
+        self._is_tilted = np.abs(self.tilt) > 1.0e-6
 
     def __repr__(self):
         s = (f'Borehole(H={self.H}, D={self.D}, r_b={self.r_b}, x={self.x},'
@@ -76,6 +80,30 @@ def distance(self, target):
                   np.sqrt((self.x - target.x)**2 + (self.y - target.y)**2))
         return dis
 
+    def is_tilted(self):
+        """
+        Returns true if the borehole is inclined.
+
+        Returns
+        -------
+        bool
+            True if borehole is inclined.
+
+        """
+        return self._is_tilted
+
+    def is_vertical(self):
+        """
+        Returns true if the borehole is vertical.
+
+        Returns
+        -------
+        bool
+            True if borehole is vertical.
+
+        """
+        return not self._is_tilted
+
     def position(self):
         """
         Returns the position of the borehole.
@@ -224,7 +252,7 @@ class _EquivalentBorehole(object):
     ----------
     boreholes : list of Borehole objects, or tuple
         Boreholes to be represented by the equivalent borehole. Alternatively,
-        tuple of attributes (H, D, r_b, x, y).
+        tuple of attributes (H, D, r_b, x, y, tilt=0., orientation=0.).
 
     Attributes
     ----------
@@ -238,6 +266,10 @@ class _EquivalentBorehole(object):
         Position (in meters) of the head of each borehole along the x-axis.
     y : (nBoreholes,) array
         Position (in meters) of the head of each borehole along the y-axis.
+    tilt : float
+        Angle (in radians) from vertical of the axis of the borehole.
+    orientation : (nBoreholes,) array
+        Direction (in radians) of the tilt of the borehole.
     nBoreholes : int
         Number of boreholes represented by the equivalent borehole.
     
@@ -256,18 +288,27 @@ def __init__(self, boreholes):
             self.r_b = boreholes[0].r_b
             self.x = np.array([b.x for b in boreholes])
             self.y = np.array([b.y for b in boreholes])
+            self.tilt = boreholes[0].tilt
+            self.orientation = np.array([b.orientation for b in boreholes])
         elif isinstance(boreholes[0], _EquivalentBorehole):
             self.H = boreholes[0].H
             self.D = boreholes[0].D
             self.r_b = boreholes[0].r_b
             self.x = np.concatenate([b.x for b in boreholes])
             self.y = np.concatenate([b.y for b in boreholes])
+            self.tilt = boreholes[0].tilt
+            self.orientation = np.concatenate(
+                [b.orientation for b in boreholes])
         elif type(boreholes) is tuple:
-            self.H, self.D, self.r_b, self.x, self.y = boreholes
+            self.H, self.D, self.r_b, self.x, self.y = boreholes[:5]
             self.x = np.atleast_1d(self.x)
             self.y = np.atleast_1d(self.y)
+            if len(boreholes)==7:
+                self.tilt, self.orientation = boreholes[5:]
 
         self.nBoreholes = len(self.x)
+        # Check if borehole is inclined
+        self._is_tilted = np.abs(self.tilt) > 1.0e-6
 
     def distance(self, target):
         """
@@ -303,6 +344,30 @@ def distance(self, target):
             self.r_b)
         return dis
 
+    def is_tilted(self):
+        """
+        Returns true if the borehole is inclined.
+
+        Returns
+        -------
+        bool
+            True if borehole is inclined.
+
+        """
+        return self._is_tilted
+
+    def is_vertical(self):
+        """
+        Returns true if the borehole is vertical.
+
+        Returns
+        -------
+        bool
+            True if borehole is vertical.
+
+        """
+        return not self._is_tilted
+
     def position(self):
         """
         Returns the position of the boreholes represented by the equivalent
@@ -350,7 +415,16 @@ def segments(self, nSegments, segment_ratios=None):
         if segment_ratios is None:
             segment_ratios = np.full(nSegments, 1. / nSegments)
         z = self._segment_edges(nSegments, segment_ratios=segment_ratios)[:-1]
-        return [_EquivalentBorehole((ratios * self.H, z_i + self.D, self.r_b, self.x, self.y)) for z_i, ratios in zip(z, segment_ratios)]
+        segments = [_EquivalentBorehole(
+            (ratios * self.H,
+             self.D + z_i * np.cos(self.tilt),
+             self.r_b,
+             self.x + z_i * np.sin(self.tilt) * np.cos(self.orientation),
+             self.y + z_i * np.sin(self.tilt) * np.sin(self.orientation),
+             self.tilt,
+             self.orientation)
+            ) for z_i, ratios in zip(z, segment_ratios)]
+        return segments
 
     def unique_distance(self, target, disTol=0.01):
         """
@@ -562,7 +636,7 @@ def remove_duplicates(boreField, disp=False):
     return new_boreField
 
 
-def rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b):
+def rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b, tilt=0., origin=None):
     """
     Build a list of boreholes in a rectangular bore field configuration.
 
@@ -582,12 +656,24 @@ def rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b):
         Borehole buried depth (in meters).
     r_b : float
         Borehole radius (in meters).
+    tilt : float, optional
+        Angle (in radians) from vertical of the axis of the borehole. The
+        orientation of the tilt is orthogonal to the origin coordinate.
+        Default is 0.
+    origin : tuple, optional
+        A coordinate indicating the origin of reference for orientation of
+        boreholes.
+        Default is the center of the rectangle.
 
     Returns
     -------
     boreField : list of Borehole objects
         List of boreholes in the rectangular bore field.
 
+    Notes
+    -----
+    Boreholes located at the origin will remain vertical.
+
     Examples
     --------
     >>> boreField = gt.boreholes.rectangle_field(N_1=3, N_2=2, B_1=5., B_2=5.,
@@ -603,14 +689,31 @@ def rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b):
     """
     borefield = []
 
+    if origin is None:
+        # When no origin is supplied, compute the origin to be at the center of
+        # the rectangle
+        x0 = (N_1 - 1) / 2 * B_1
+        y0 = (N_2 - 1) / 2 * B_2
+    else:
+        x0, y0 = origin
+
     for j in range(N_2):
         for i in range(N_1):
-            borefield.append(Borehole(H, D, r_b, x=i*B_1, y=j*B_2))
+            x = i * B_1
+            y = j * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
 
     return borefield
 
 
-def L_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
+def L_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b, tilt=0., origin=None):
     """
     Build a list of boreholes in a L-shaped bore field configuration.
 
@@ -630,12 +733,24 @@ def L_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
         Borehole buried depth (in meters).
     r_b : float
         Borehole radius (in meters).
+    tilt : float, optional
+        Angle (in radians) from vertical of the axis of the borehole. The
+        orientation of the tilt is orthogonal to the origin coordinate.
+        Default is 0.
+    origin : tuple, optional
+        A coordinate indicating the origin of reference for orientation of
+        boreholes.
+        Default is origin is placed at the center of an assumed rectangle.
 
     Returns
     -------
     boreField : list of Borehole objects
         List of boreholes in the L-shaped bore field.
 
+    Notes
+    -----
+    Boreholes located at the origin will remain vertical.
+
     Examples
     --------
     >>> boreField = gt.boreholes.L_shaped_field(N_1=3, N_2=2, B_1=5., B_2=5.,
@@ -651,15 +766,41 @@ def L_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
     """
     borefield = []
 
+    if origin is None:
+        # When no origin is supplied, compute the origin to be at the center of
+        # the rectangle
+        x0 = (N_1 - 1) / 2 * B_1
+        y0 = (N_2 - 1) / 2 * B_2
+    else:
+        x0, y0 = origin
+
     for i in range(N_1):
-        borefield.append(Borehole(H, D, r_b, x=i*B_1, y=0.))
+        x = i * B_1
+        y = 0.
+        # The borehole is inclined only if it does not lie on the origin
+        if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+            orientation = np.arctan2(y - y0, x - x0)
+            borefield.append(
+                Borehole(
+                    H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+        else:
+            borefield.append(Borehole(H, D, r_b, x, y))
     for j in range(1, N_2):
-        borefield.append(Borehole(H, D, r_b, x=0., y=j*B_2))
+        x = 0.
+        y = j * B_2
+        # The borehole is inclined only if it does not lie on the origin
+        if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+            orientation = np.arctan2(y - y0, x - x0)
+            borefield.append(
+                Borehole(
+                    H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+        else:
+            borefield.append(Borehole(H, D, r_b, x, y))
 
     return borefield
 
 
-def U_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
+def U_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b, tilt=0., origin=None):
     """
     Build a list of boreholes in a U-shaped bore field configuration.
 
@@ -679,12 +820,24 @@ def U_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
         Borehole buried depth (in meters).
     r_b : float
         Borehole radius (in meters).
+    tilt : float, optional
+        Angle (in radians) from vertical of the axis of the borehole. The
+        orientation of the tilt is orthogonal to the origin coordinate.
+        Default is 0.
+    origin : tuple, optional
+        A coordinate indicating the origin of reference for orientation of
+        boreholes.
+        Default is the center considering an outer rectangle.
 
     Returns
     -------
     boreField : list of Borehole objects
         List of boreholes in the U-shaped bore field.
 
+    Notes
+    -----
+    Boreholes located at the origin will remain vertical.
+
     Examples
     --------
     >>> boreField = gt.boreholes.U_shaped_field(N_1=3, N_2=2, B_1=5., B_2=5.,
@@ -700,19 +853,55 @@ def U_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
     """
     borefield = []
 
+    if origin is None:
+        # When no origin is supplied, compute the origin to be at the center of
+        # the rectangle
+        x0 = (N_1 - 1) / 2 * B_1
+        y0 = (N_2 - 1) / 2 * B_2
+    else:
+        x0, y0 = origin
+
     if N_1 > 2 and N_2 > 1:
         for i in range(N_1):
-            borefield.append(Borehole(H, D, r_b, x=i*B_1, y=0.))
+            x = i * B_1
+            y = 0.
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
         for j in range(1, N_2):
-            borefield.append(Borehole(H, D, r_b, x=0, y=j*B_2))
-            borefield.append(Borehole(H, D, r_b, x=(N_1-1)*B_1, y=j*B_2))
+            x = 0.
+            y = j * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
+            x = (N_1 - 1) * B_1
+            y = j * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
     else:
-        borefield = rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b)
+        borefield = rectangle_field(
+            N_1, N_2, B_1, B_2, H, D, r_b, tilt=tilt, origin=origin)
 
     return borefield
 
 
-def box_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
+def box_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b, tilt=0, origin=None):
     """
     Build a list of boreholes in a box-shaped bore field configuration.
 
@@ -732,12 +921,24 @@ def box_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
         Borehole buried depth (in meters).
     r_b : float
         Borehole radius (in meters).
+    tilt : float, optional
+        Angle (in radians) from vertical of the axis of the borehole. The
+        orientation of the tilt is orthogonal to the origin coordinate.
+        Default is 0.
+    origin : tuple, optional
+        A coordinate indicating the origin of reference for orientation of
+        boreholes.
+        Default is the center of the box.
 
     Returns
     -------
     boreField : list of Borehole objects
         List of boreholes in the box-shaped bore field.
 
+    Notes
+    -----
+    Boreholes located at the origin will remain vertical.
+
     Examples
     --------
     >>> boreField = gt.boreholes.box_shaped_field(N_1=4, N_2=3, B_1=5., B_2=5.,
@@ -755,22 +956,65 @@ def box_shaped_field(N_1, N_2, B_1, B_2, H, D, r_b):
     """
     borefield = []
 
+    if origin is None:
+        # When no origin is supplied, compute the origin to be at the center of
+        # the rectangle
+        x0 = (N_1 - 1) / 2 * B_1
+        y0 = (N_2 - 1) / 2 * B_2
+    else:
+        x0, y0 = origin
 
     if N_1 > 2 and N_2 > 2:
         for i in range(N_1):
-            borefield.append(Borehole(H, D, r_b, x=i*B_1, y=0.))
+            x = i * B_1
+            y = 0.
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
+            x = i * B_1
+            y = (N_2 - 1) * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
         for j in range(1, N_2-1):
-            borefield.append(Borehole(H, D, r_b, x=0., y=j*B_2))
-            borefield.append(Borehole(H, D, r_b, x=(N_1-1)*B_1, y=j*B_2))
-        for i in range(N_1):
-            borefield.append(Borehole(H, D, r_b, x=i*B_1, y=(N_2-1)*B_2))
+            x = 0.
+            y = j * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
+            x = (N_1 - 1) * B_1
+            y = j * B_2
+            # The borehole is inclined only if it does not lie on the origin
+            if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+                orientation = np.arctan2(y - y0, x - x0)
+                borefield.append(
+                    Borehole(
+                        H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+            else:
+                borefield.append(Borehole(H, D, r_b, x, y))
     else:
-        borefield = rectangle_field(N_1, N_2, B_1, B_2, H, D, r_b)
+        borefield = rectangle_field(
+            N_1, N_2, B_1, B_2, H, D, r_b, tilt=tilt, origin=origin)
 
     return borefield
 
 
-def circle_field(N, R, H, D, r_b):
+def circle_field(N, R, H, D, r_b, tilt=0., origin=None):
     """
     Build a list of boreholes in a circular field configuration.
 
@@ -786,12 +1030,24 @@ def circle_field(N, R, H, D, r_b):
         Borehole buried depth (in meters).
     r_b : float
         Borehole radius (in meters).
+    tilt : float, optional
+        Angle (in radians) from vertical of the axis of the borehole. The
+        orientation of the tilt is towards the exterior of the bore field.
+        Default is 0.
+    origin : tuple
+        A coordinate indicating the origin of reference for orientation of
+        boreholes.
+        Default is the origin (0, 0).
 
     Returns
     -------
     boreField : list of Borehole objects
         List of boreholes in the circular shaped bore field.
 
+    Notes
+    -----
+    Boreholes located at the origin will remain vertical.
+
     Examples
     --------
     >>> boreField = gt.boreholes.circle_field(N=8, R = 5., H=100., D=2.5,
@@ -811,9 +1067,26 @@ def circle_field(N, R, H, D, r_b):
     """
     borefield = []
 
+    if origin is None:
+        # When no origin is supplied, compute the origin to be at the center of
+        # the rectangle
+        x0 = 0.
+        y0 = 0.
+    else:
+        x0, y0 = origin
+
     for i in range(N):
-        borefield.append(Borehole(H, D, r_b, x=R*np.cos(2*pi*i/N),
-                                  y=R*np.sin(2*pi*i/N)))
+        x = R * np.cos(2 * pi * i / N)
+        y = R * np.sin(2 * pi * i / N)
+        orientation = np.arctan2(y - y0, x - x0)
+        # The borehole is inclined only if it does not lie on the origin
+        if np.sqrt((x - x0)**2 + (y - y0)**2) > r_b:
+            orientation = np.arctan2(y - y0, x - x0)
+            borefield.append(
+                Borehole(
+                    H, D, r_b, x, y, tilt=tilt, orientation=orientation))
+        else:
+            borefield.append(Borehole(H, D, r_b, x, y))
 
     return borefield
 
@@ -837,12 +1110,12 @@ def field_from_file(filename):
     -----
     The text file should be formatted as follows::
 
-        # x   y     H     D     r_b
-        0.    0.    100.  2.5   0.075
-        5.    0.    100.  2.5   0.075
-        0.    5.    100.  2.5   0.075
-        0.    10.   100.  2.5   0.075
-        0.    20.   100.  2.5   0.075
+        # x   y     H     D     r_b     tilt   orientation
+        0.    0.    100.  2.5   0.075   0.     0.
+        5.    0.    100.  2.5   0.075   0.     0.
+        0.    5.    100.  2.5   0.075   0.     0.
+        0.    10.   100.  2.5   0.075   0.     0.
+        0.    20.   100.  2.5   0.075   0.     0.
 
     """
     # Load data from file
@@ -855,12 +1128,22 @@ def field_from_file(filename):
         H = line[2]
         D = line[3]
         r_b = line[4]
-        borefield.append(Borehole(H, D, r_b, x=x, y=y))
+        # Previous versions of pygfunction only required up to line[4].
+        # Now check to see if tilt and orientation exist.
+        if len(line) == 6:
+            tilt = line[5]
+            orientation = line[6]
+        else:
+            tilt = 0.
+            orientation = 0.
+        borefield.append(
+            Borehole(H, D, r_b, x=x, y=y, tilt=tilt, orientation=orientation))
 
     return borefield
 
 
-def visualize_field(borefield, viewTop=True, view3D=True, labels=True):
+def visualize_field(
+        borefield, viewTop=True, view3D=True, labels=True, showTilt=True):
     """
     Plot the top view and 3D view of borehole positions.
 
@@ -868,15 +1151,18 @@ def visualize_field(borefield, viewTop=True, view3D=True, labels=True):
     ----------
     borefield : list
         List of boreholes in the bore field.
-    viewTop : bool
+    viewTop : bool, optional
         Set to True to plot top view.
         Default is True
-    view3D : bool
+    view3D : bool, optional
         Set to True to plot 3D view.
         Default is True
-    labels : bool
+    labels : bool, optional
         Set to True to annotate borehole indices to top view plot.
         Default is True
+    showTilt : bool, optional
+        Set to True to show borehole inclination on top view plot.
+        Default is True
 
     Returns
     -------
@@ -913,8 +1199,17 @@ def visualize_field(borefield, viewTop=True, view3D=True, labels=True):
     if viewTop:
         i = 0   # Initialize borehole index
         for borehole in borefield:
-            (x, y) = borehole.position()    # Extract borehole position
+            # Extract borehole parameters
+            (x, y) = borehole.position()
+            H = borehole.H
+            tilt = borehole.tilt
+            orientation = borehole.orientation
             # Add current borehole to the figure
+            if showTilt:
+                ax1.plot(
+                    [x, x + H * np.sin(tilt) * np.cos(orientation)],
+                    [y, y + H * np.sin(tilt) * np.sin(orientation)],
+                    'k--')
             ax1.plot(x, y, 'ko')
             if labels: ax1.text(x, y,
                                 f' {i}',
diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py
index 98384d80..99e1d8ad 100644
--- a/pygfunction/gfunction.py
+++ b/pygfunction/gfunction.py
@@ -10,7 +10,8 @@
 
 from .boreholes import Borehole, _EquivalentBorehole, find_duplicates
 from .heat_transfer import finite_line_source, finite_line_source_vectorized, \
-    finite_line_source_equivalent_boreholes_vectorized
+    finite_line_source_equivalent_boreholes_vectorized, \
+    finite_line_source_inclined_vectorized
 from .networks import Network, _EquivalentNetwork, network_thermal_resistance
 from .utilities import _initialize_figure, _format_axes
 from . import utilities
@@ -110,6 +111,10 @@ class gFunction(object):
                 Number of terms in the approximation of the FLS solution. This
                 parameter is unused if `approximate_FLS` is set to False.
                 Default is 10. Maximum is 25.
+            mQuad : int, optional
+                Number of Gauss-Legendre sample points for the integral over
+                :math:`u` in the inclined FLS solution.
+                Default is 11.
             disp : bool, optional
                 Set to true to print progression messages.
                 Default is False.
@@ -164,6 +169,12 @@ class gFunction(object):
                 method.
                 Default is 1.
 
+    Notes
+    -----
+    - The 'equivalent' solver does not support the 'MIFT' boundary condition
+      when boreholes are connected in series.
+    - The 'equivalent' solver does not support inclined boreholes.
+
     References
     ----------
     .. [#gFunction-CimBer2014] Cimmino, M., & Bernier, M. (2014). A
@@ -290,7 +301,7 @@ def visualize_g_function(self):
         plt.tight_layout()
         return fig
 
-    def visualize_heat_extraction_rates(self, iBoreholes=None):
+    def visualize_heat_extraction_rates(self, iBoreholes=None, showTilt=True):
         """
         Plot the time-variation of the average heat extraction rates.
 
@@ -301,6 +312,9 @@ def visualize_heat_extraction_rates(self, iBoreholes=None):
             If iBoreholes is None, heat extraction rates are plotted for all
             boreholes.
             Default is None.
+        showTilt : bool
+            Set to True to show borehole inclination.
+            Default is True
 
         Returns
         -------
@@ -337,6 +351,13 @@ def visualize_heat_extraction_rates(self, iBoreholes=None):
                 line = ax2.plot(lntts, Q_t[iBoreholes.index(i)])
                 color = line[-1]._color
                 # Draw colored marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color=color)
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -344,6 +365,13 @@ def visualize_heat_extraction_rates(self, iBoreholes=None):
                          color=color)
             else:
                 # Draw black marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color='k')
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -355,7 +383,7 @@ def visualize_heat_extraction_rates(self, iBoreholes=None):
         return fig
 
     def visualize_heat_extraction_rate_profiles(
-            self, time=None, iBoreholes=None):
+            self, time=None, iBoreholes=None, showTilt=True):
         """
         Plot the heat extraction rate profiles at chosen time.
 
@@ -371,6 +399,9 @@ def visualize_heat_extraction_rate_profiles(
             If iBoreholes is None, heat extraction rates are plotted for all
             boreholes.
             Default is None.
+        showTilt : bool
+            Set to True to show borehole inclination.
+            Default is True
 
         Returns
         -------
@@ -405,6 +436,13 @@ def visualize_heat_extraction_rate_profiles(
                     Q_b[iBoreholes.index(i)], z[iBoreholes.index(i)])
                 color = line[-1]._color
                 # Draw colored marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color=color)
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -412,6 +450,13 @@ def visualize_heat_extraction_rate_profiles(
                          color=color)
             else:
                 # Draw black marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color='k')
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -421,7 +466,7 @@ def visualize_heat_extraction_rate_profiles(
         plt.tight_layout()
         return fig
 
-    def visualize_temperatures(self, iBoreholes=None):
+    def visualize_temperatures(self, iBoreholes=None, showTilt=True):
         """
         Plot the time-variation of the average borehole wall temperatures.
 
@@ -431,6 +476,9 @@ def visualize_temperatures(self, iBoreholes=None):
             Borehole indices to plot temperatures.
             If iBoreholes is None, temperatures are plotted for all boreholes.
             Default is None.
+        showTilt : bool
+            Set to True to show borehole inclination.
+            Default is True
 
         Returns
         -------
@@ -467,6 +515,13 @@ def visualize_temperatures(self, iBoreholes=None):
                 line = ax2.plot(lntts, T_b[iBoreholes.index(i)])
                 color = line[-1]._color
                 # Draw colored marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color=color)
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -474,6 +529,13 @@ def visualize_temperatures(self, iBoreholes=None):
                          color=color)
             else:
                 # Draw black marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color='k')
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -484,7 +546,8 @@ def visualize_temperatures(self, iBoreholes=None):
         plt.tight_layout()
         return fig
 
-    def visualize_temperature_profiles(self, time=None, iBoreholes=None):
+    def visualize_temperature_profiles(
+            self, time=None, iBoreholes=None, showTilt=True):
         """
         Plot the borehole wall temperature profiles at chosen time.
 
@@ -498,6 +561,9 @@ def visualize_temperature_profiles(self, time=None, iBoreholes=None):
             Borehole indices to plot temperature profiles.
             If iBoreholes is None, temperatures are plotted for all boreholes.
             Default is None.
+        showTilt : bool
+            Set to True to show borehole inclination.
+            Default is True
 
         Returns
         -------
@@ -532,6 +598,13 @@ def visualize_temperature_profiles(self, time=None, iBoreholes=None):
                     T_b[iBoreholes.index(i)], z[iBoreholes.index(i)])
                 color = line[-1]._color
                 # Draw colored marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color=color)
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -539,6 +612,13 @@ def visualize_temperature_profiles(self, time=None, iBoreholes=None):
                          color=color)
             else:
                 # Draw black marker for borehole position
+                if showTilt:
+                    ax1.plot(
+                        [borehole.x, borehole.x + borehole.H*np.sin(borehole.tilt)*np.cos(borehole.orientation)],
+                        [borehole.y, borehole.y + borehole.H*np.sin(borehole.tilt)*np.sin(borehole.orientation)],
+                         linestyle='--',
+                         marker='None',
+                         color='k')
                 ax1.plot(borehole.x,
                          borehole.y,
                          linestyle='None',
@@ -1327,6 +1407,10 @@ class _BaseSolver(object):
         Number of terms in the approximation of the FLS solution. This
         parameter is unused if `approximate_FLS` is set to False.
         Default is 10. Maximum is 25.
+    mQuad : int, optional
+        Number of Gauss-Legendre sample points for the integral over :math:`u`
+        in the inclined FLS solution.
+        Default is 11.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -1346,8 +1430,9 @@ class _BaseSolver(object):
     """
     def __init__(self, boreholes, network, time, boundary_condition,
                  nSegments=8, segment_ratios=utilities.segment_ratios,
-                 approximate_FLS=False, nFLS=10, disp=False, profiles=False,
-                 kind='linear', dtype=np.double, **other_options):
+                 approximate_FLS=False, mQuad=11, nFLS=10, disp=False,
+                 profiles=False, kind='linear', dtype=np.double,
+                 **other_options):
         self.boreholes = boreholes
         self.network = network
         # Convert time to a 1d array
@@ -1382,6 +1467,7 @@ def __init__(self, boreholes, network, time, boundary_condition,
         self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)])
                             for i in range(nBoreholes)]
         self.approximate_FLS = approximate_FLS
+        self.mQuad = mQuad
         self.nFLS = nFLS
         self.disp = disp
         self.profiles = profiles
@@ -1703,9 +1789,6 @@ def _check_inputs(self):
         assert np.all([isinstance(b, Borehole) for b in self.boreholes]), \
             "The list of boreholes contains elements that are not Borehole " \
             "objects."
-        assert np.all(np.abs([b.tilt for b in self.boreholes]) <= 1e-6), \
-            "The current version of pygfunction only supports vertical " \
-            "boreholes."
         assert self.network is None or isinstance(self.network, Network), \
             "The network is not a valid 'Network' object."
         assert self.network is None or (self.network.m_flow_network is not None and self.network.cp_f is not None), \
@@ -1813,6 +1896,10 @@ class _Detailed(_BaseSolver):
         Number of terms in the approximation of the FLS solution. This
         parameter is unused if `approximate_FLS` is set to False.
         Default is 10. Maximum is 25.
+    mQuad : int, optional
+        Number of Gauss-Legendre sample points for the integral over :math:`u`
+        in the inclined FLS solution.
+        Default is 11.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -1871,7 +1958,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
         of thermal response factors, containing a copy of the matrix accessible
         by h_ij.y[:nSources,:nSources,:nt+1]. The first index along the
         third axis corresponds to time t=0. The interp1d object can be used to
-        obtain thermal response factors at any intermediat time by
+        obtain thermal response factors at any intermediate time by
         h_ij(t)[:nSources,:nSources].
 
         Attributes
@@ -1908,39 +1995,29 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
         # Segment-to-segment thermal response factors for same-borehole
         # thermal interactions
         # ---------------------------------------------------------------------
-        h_ii, i_segment, j_segment = \
+        h, i_segment, j_segment = \
             self._thermal_response_factors_borehole_to_self(time, alpha)
         # Broadcast values to h_ij matrix
-        h_ij[j_segment, i_segment, 1:] = h_ii
+        h_ij[j_segment, i_segment, 1:] = h
         # ---------------------------------------------------------------------
         # Segment-to-segment thermal response factors for
         # borehole-to-borehole thermal interactions
         # ---------------------------------------------------------------------
-        for i, (borehole, nSegments, ratios, i0, i1) in enumerate(
-                zip(self.boreholes,
-                    self.nBoreSegments,
-                    self.segment_ratios,
-                    self._i0Segments,
-                    self._i1Segments)):
+        for i, (i0, i1) in enumerate(zip(self._i0Segments, self._i1Segments)):
             # Segments of the receiving borehole
-            b2 = borehole.segments(nSegments, segment_ratios=ratios)
+            b2 = self.boreSegments[i0:i1]
             if i+1 < nBoreholes:
                 # Segments of the emitting borehole
                 b1 = self.boreSegments[i1:]
                 h = finite_line_source(
                     time, alpha, b1, b2, approximation=self.approximate_FLS,
-                    N=self.nFLS)
+                    N=self.nFLS, M=self.mQuad)
                 # Broadcast values to h_ij matrix
-                for j, (j0, j1) in enumerate(
-                        zip(self._i0Segments[i+1:], self._i1Segments[i+1:]),
-                        start=i+1):
-                    h_ij[i0:i1, j0:j1, 1:] = h[:, j0-i1:j1-i1, :]
-                    if j > i:
-                        h_ij[j0:j1, i0:i1, 1:] = np.einsum(
-                            'ijk,i,j->jik',
-                            h[:, j0-i1:j1-i1, :],
-                            segment_lengths[i0:i1],
-                            1/segment_lengths[j0:j1])
+                h_ij[i0:i1, i1:, 1:] = h
+                h_ij[i1:, i0:i1, 1:] = \
+                    np.swapaxes(h, 0, 1) * np.divide.outer(
+                        segment_lengths[i0:i1],
+                        segment_lengths[i1:]).T[:,:,np.newaxis]
 
         # Return 2d array if time is a scalar
         if np.isscalar(time):
@@ -1968,7 +2045,7 @@ def _thermal_response_factors_borehole_to_self(self, time, alpha):
 
         Returns
         -------
-        h_ii : array
+        h : array
             Finite line source solution.
         i_segment : list
             Indices of the emitting segments in the bore field.
@@ -1997,11 +2074,22 @@ def _thermal_response_factors_borehole_to_self(self, time, alpha):
             np.sqrt((x[i_segment] - x[j_segment])**2 + (y[i_segment] - y[j_segment])**2),
             r_b[i_segment])
         # FlS solution
-        h_ii = finite_line_source_vectorized(
-            time, alpha,
-            dis, H[i_segment], D[i_segment], H[j_segment], D[j_segment],
-            approximation=self.approximate_FLS, N=self.nFLS)
-        return h_ii, i_segment, j_segment
+        if np.all([b.is_vertical() for b in self.boreholes]):
+            h = finite_line_source_vectorized(
+                time, alpha,
+                dis, H[i_segment], D[i_segment], H[j_segment], D[j_segment],
+                approximation=self.approximate_FLS, N=self.nFLS)
+        else:
+            tilt = np.array([b.tilt for b in self.boreSegments])
+            orientation = np.array([b.orientation for b in self.boreSegments])
+            h = finite_line_source_inclined_vectorized(
+                time, alpha,
+                r_b[i_segment], x[i_segment], y[i_segment], H[i_segment],
+                D[i_segment], tilt[i_segment], orientation[i_segment],
+                x[j_segment], y[j_segment], H[j_segment], D[j_segment],
+                tilt[j_segment], orientation[j_segment], M=self.mQuad,
+                approximation=self.approximate_FLS, N=self.nFLS)
+        return h, i_segment, j_segment
         
 
 
@@ -2064,6 +2152,10 @@ class _Similarities(_BaseSolver):
         Number of terms in the approximation of the FLS solution. This
         parameter is unused if `approximate_FLS` is set to False.
         Default is 10. Maximum is 25.
+    mQuad : int, optional
+        Number of Gauss-Legendre sample points for the integral over :math:`u`
+        in the inclined FLS solution.
+        Default is 11.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -2175,49 +2267,43 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
         tic = perf_counter()
         # Initialize segment-to-segment response factors
         h_ij = np.zeros((self.nSources, self.nSources, nt+1), dtype=self.dtype)
-        segment_lengths = self.segment_lengths()
 
         # ---------------------------------------------------------------------
         # Segment-to-segment thermal response factors for same-borehole thermal
-        # interactions
+        # interactions (vertical boreholes)
         # ---------------------------------------------------------------------
-        for group in self.borehole_to_self:
-            # Index of first borehole in group
-            i = group[0]
-            # Find segment-to-segment similarities
-            H1, D1, H2, D2, i_pair, j_pair, k_pair = \
-                self._map_axial_segment_pairs(i, i)
-            # Locate thermal response factors in the h_ij matrix
-            i_segment, j_segment, k_segment, l_segment = \
-                self._map_segment_pairs(
-                    i_pair, j_pair, k_pair, [(n, n) for n in group], [0])
-            # Evaluate FLS at all time steps
-            H1 = H1.reshape(1, -1)
-            H2 = H2.reshape(1, -1)
-            D1 = D1.reshape(1, -1)
-            D2 = D2.reshape(1, -1)
-            dis = self.boreholes[i].r_b
-            h = finite_line_source_vectorized(
-                time, alpha, dis, H1, D1, H2, D2,
-                approximation=self.approximate_FLS, N=self.nFLS)
-            # Broadcast values to h_ij matrix
-            h_ij[j_segment, i_segment, 1:] = h[0, k_segment, :]
+        # Evaluate FLS at all time steps
+        h, i_segment, j_segment, k_segment = \
+            self._thermal_response_factors_borehole_to_self_vertical(
+                time, alpha)
+        # Broadcast values to h_ij matrix
+        h_ij[j_segment, i_segment, 1:] = h[k_segment, :]
+        # ---------------------------------------------------------------------
+        # Segment-to-segment thermal response factors for same-borehole thermal
+        # interactions (inclined boreholes)
+        # ---------------------------------------------------------------------
+        # Evaluate FLS at all time steps
+        h, i_segment, j_segment, k_segment = \
+            self._thermal_response_factors_borehole_to_self_inclined(
+                time, alpha)
+        # Broadcast values to h_ij matrix
+        h_ij[j_segment, i_segment, 1:] = h[k_segment, :]
         # ---------------------------------------------------------------------
         # Segment-to-segment thermal response factors for borehole-to-borehole
-        # thermal interactions
+        # thermal interactions (vertical boreholes)
         # ---------------------------------------------------------------------
         for pairs, distances, distance_indices in zip(
-                self.borehole_to_borehole,
-                self.borehole_to_borehole_distances,
-                self.borehole_to_borehole_indices):
+                self.borehole_to_borehole_vertical,
+                self.borehole_to_borehole_distances_vertical,
+                self.borehole_to_borehole_indices_vertical):
             # Index of first borehole pair in group
             i, j = pairs[0]
             # Find segment-to-segment similarities
             H1, D1, H2, D2, i_pair, j_pair, k_pair = \
-                self._map_axial_segment_pairs(i, j)
+                self._map_axial_segment_pairs_vertical(i, j)
             # Locate thermal response factors in the h_ij matrix
             i_segment, j_segment, k_segment, l_segment = \
-                self._map_segment_pairs(
+                self._map_segment_pairs_vertical(
                     i_pair, j_pair, k_pair, pairs, distance_indices)
             # Evaluate FLS at all time steps
             dis = np.reshape(distances, (-1, 1))
@@ -2237,6 +2323,17 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
                 h_ij[i_segment, j_segment, 1:] = h[l_segment, k_segment, :]
             else:
                 h_ij[i_segment, j_segment, 1:] = (h * H2.T / H1.T)[l_segment, k_segment, :]
+        # ---------------------------------------------------------------------
+        # Segment-to-segment thermal response factors for borehole-to-borehole
+        # thermal interactions (inclined boreholes)
+        # ---------------------------------------------------------------------
+        # Evaluate FLS at all time steps
+        h, hT, i_segment, j_segment, k_segment = \
+            self._thermal_response_factors_borehole_to_borehole_inclined(
+                time, alpha)
+        # Broadcast values to h_ij matrix
+        h_ij[j_segment, i_segment, 1:] = h[k_segment, :]
+        h_ij[i_segment, j_segment, 1:] = hT[k_segment, :]
 
         # Return 2d array if time is a scalar
         if np.isscalar(time):
@@ -2251,6 +2348,227 @@ def thermal_response_factors(self, time, alpha, kind='linear'):
 
         return h_ij
 
+    def _thermal_response_factors_borehole_to_borehole_inclined(
+            self, time, alpha):
+        """
+        Evaluate the segment-to-segment thermal response factors for all pairs
+        of inclined segments.
+
+        Attributes
+        ----------
+        time : float or array
+            Values of time (in seconds) for which the g-function is evaluated.
+        alpha : float
+            Soil thermal diffusivity (in m2/s).
+
+        Returns
+        -------
+        h : array
+            Finite line source solution.
+        hT : array
+            Reciprocal finite line source solution.
+        i_segment : list
+            Indices of the emitting segments in the bore field.
+        j_segment : list
+            Indices of the receiving segments in the bore field.
+        k_segment : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair) in the bore field.
+        """
+        rb1 = np.array([])
+        x1 = np.array([])
+        y1 = np.array([])
+        H1 = np.array([])
+        D1 = np.array([])
+        tilt1 = np.array([])
+        orientation1 = np.array([])
+        x2 = np.array([])
+        y2 = np.array([])
+        H2 = np.array([])
+        D2 = np.array([])
+        tilt2 = np.array([])
+        orientation2 = np.array([])
+        i_segment = np.array([], dtype=np.uint)
+        j_segment = np.array([], dtype=np.uint)
+        k_segment = np.array([], dtype=np.uint)
+        k0 = 0
+        for pairs in self.borehole_to_borehole_inclined:
+            # Index of first borehole pair in group
+            i, j = pairs[0]
+            # Find segment-to-segment similarities
+            rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \
+                x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \
+                i_pair, j_pair, k_pair = \
+                    self._map_axial_segment_pairs_inclined(i, j)
+            # Locate thermal response factors in the h_ij matrix
+            i_segment_i, j_segment_i, k_segment_i = \
+                self._map_segment_pairs_inclined(i_pair, j_pair, k_pair, pairs)
+            # Append lists
+            rb1 = np.append(rb1, rb1_i)
+            x1 = np.append(x1, x1_i)
+            y1 = np.append(y1, y1_i)
+            H1 = np.append(H1, H1_i)
+            D1 = np.append(D1, D1_i)
+            tilt1 = np.append(tilt1, tilt1_i)
+            orientation1 = np.append(orientation1, orientation1_i)
+            x2 = np.append(x2, x2_i)
+            y2 = np.append(y2, y2_i)
+            H2 = np.append(H2, H2_i)
+            D2 = np.append(D2, D2_i)
+            tilt2 = np.append(tilt2, tilt2_i)
+            orientation2 = np.append(orientation2, orientation2_i)
+            i_segment = np.append(i_segment, i_segment_i)
+            j_segment = np.append(j_segment, j_segment_i)
+            k_segment = np.append(k_segment, k_segment_i + k0)
+            k0 += len(k_pair)
+        # Evaluate FLS at all time steps
+        h = finite_line_source_inclined_vectorized(
+            time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1,
+            x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad,
+            approximation=self.approximate_FLS, N=self.nFLS)
+        hT = (h.T * H2 / H1).T
+        return h, hT, i_segment, j_segment, k_segment
+
+    def _thermal_response_factors_borehole_to_self_inclined(self, time, alpha):
+        """
+        Evaluate the segment-to-segment thermal response factors for all pairs
+        of segments between each inclined borehole and itself.
+
+        Attributes
+        ----------
+        time : float or array
+            Values of time (in seconds) for which the g-function is evaluated.
+        alpha : float
+            Soil thermal diffusivity (in m2/s).
+
+        Returns
+        -------
+        h : array
+            Finite line source solution.
+        i_segment : list
+            Indices of the emitting segments in the bore field.
+        j_segment : list
+            Indices of the receiving segments in the bore field.
+        k_segment : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair) in the bore field.
+        """
+        rb1 = np.array([])
+        x1 = np.array([])
+        y1 = np.array([])
+        H1 = np.array([])
+        D1 = np.array([])
+        tilt1 = np.array([])
+        orientation1 = np.array([])
+        x2 = np.array([])
+        y2 = np.array([])
+        H2 = np.array([])
+        D2 = np.array([])
+        tilt2 = np.array([])
+        orientation2 = np.array([])
+        i_segment = np.array([], dtype=np.uint)
+        j_segment = np.array([], dtype=np.uint)
+        k_segment = np.array([], dtype=np.uint)
+        k0 = 0
+        for group in self.borehole_to_self_inclined:
+            # Index of first borehole in group
+            i = group[0]
+            # Find segment-to-segment similarities
+            rb1_i, x1_i, y1_i, H1_i, D1_i, tilt1_i, orientation1_i, \
+                x2_i, y2_i, H2_i, D2_i, tilt2_i, orientation2_i, \
+                i_pair, j_pair, k_pair = \
+                    self._map_axial_segment_pairs_inclined(i, i)
+            # Locate thermal response factors in the h_ij matrix
+            i_segment_i, j_segment_i, k_segment_i = \
+                self._map_segment_pairs_inclined(
+                    i_pair, j_pair, k_pair, [(n, n) for n in group])
+            # Append lists
+            rb1 = np.append(rb1, rb1_i)
+            x1 = np.append(x1, x1_i)
+            y1 = np.append(y1, y1_i)
+            H1 = np.append(H1, H1_i)
+            D1 = np.append(D1, D1_i)
+            tilt1 = np.append(tilt1, tilt1_i)
+            orientation1 = np.append(orientation1, orientation1_i)
+            x2 = np.append(x2, x2_i)
+            y2 = np.append(y2, y2_i)
+            H2 = np.append(H2, H2_i)
+            D2 = np.append(D2, D2_i)
+            tilt2 = np.append(tilt2, tilt2_i)
+            orientation2 = np.append(orientation2, orientation2_i)
+            i_segment = np.append(i_segment, i_segment_i)
+            j_segment = np.append(j_segment, j_segment_i)
+            k_segment = np.append(k_segment, k_segment_i + k0)
+            k0 += len(k_pair)
+        # Evaluate FLS at all time steps
+        h = finite_line_source_inclined_vectorized(
+            time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1,
+            x2, y2, H2, D2, tilt2, orientation2, M=self.mQuad,
+            approximation=self.approximate_FLS, N=self.nFLS)
+        return h, i_segment, j_segment, k_segment
+
+    def _thermal_response_factors_borehole_to_self_vertical(self, time, alpha):
+        """
+        Evaluate the segment-to-segment thermal response factors for all pairs
+        of segments between each vertical borehole and itself.
+
+        Attributes
+        ----------
+        time : float or array
+            Values of time (in seconds) for which the g-function is evaluated.
+        alpha : float
+            Soil thermal diffusivity (in m2/s).
+
+        Returns
+        -------
+        h : array
+            Finite line source solution.
+        i_segment : list
+            Indices of the emitting segments in the bore field.
+        j_segment : list
+            Indices of the receiving segments in the bore field.
+        k_segment : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair) in the bore field.
+        """
+        H1 = np.array([])
+        D1 = np.array([])
+        H2 = np.array([])
+        D2 = np.array([])
+        dis = np.array([])
+        i_segment = np.array([], dtype=np.uint)
+        j_segment = np.array([], dtype=np.uint)
+        k_segment = np.array([], dtype=np.uint)
+        k0 = 0
+        for group in self.borehole_to_self_vertical:
+            # Index of first borehole in group
+            i = group[0]
+            # Find segment-to-segment similarities
+            H1_i, D1_i, H2_i, D2_i, i_pair, j_pair, k_pair = \
+                self._map_axial_segment_pairs_vertical(i, i)
+            # Locate thermal response factors in the h_ij matrix
+            i_segment_i, j_segment_i, k_segment_i, l_segment_i = \
+                self._map_segment_pairs_vertical(
+                    i_pair, j_pair, k_pair, [(n, n) for n in group], [0])
+            # Append lists
+            H1 = np.append(H1, H1_i)
+            D1 = np.append(D1, D1_i)
+            H2 = np.append(H2, H2_i)
+            D2 = np.append(D2, D2_i)
+            if len(self.borehole_to_self_vertical) > 1:
+                dis = np.append(dis, np.full(len(H1_i), self.boreholes[i].r_b))
+            else:
+                dis = self.boreholes[i].r_b
+            i_segment = np.append(i_segment, i_segment_i)
+            j_segment = np.append(j_segment, j_segment_i)
+            k_segment = np.append(k_segment, k_segment_i + k0)
+            k0 += np.max(k_pair) + 1
+        # Evaluate FLS at all time steps
+        h = finite_line_source_vectorized(
+            time, alpha, dis, H1, D1, H2, D2,
+            approximation=self.approximate_FLS, N=self.nFLS)
+        return h, i_segment, j_segment, k_segment
+
     def find_similarities(self):
         """
         Find similarities in the FLS solution for groups of boreholes.
@@ -2265,12 +2583,13 @@ def find_similarities(self):
 
         # Find similar pairs of boreholes
         # Boreholes can only be similar if their segments are similar
-        self.borehole_to_self, self.borehole_to_borehole = \
-            self._find_axial_borehole_pairs(self.boreholes)
-        # Find distances for each similar pairs
-        self.borehole_to_borehole_distances, self.borehole_to_borehole_indices = \
+        self.borehole_to_self_vertical, self.borehole_to_self_inclined, \
+            self.borehole_to_borehole_vertical, self.borehole_to_borehole_inclined = \
+                self._find_axial_borehole_pairs(self.boreholes)
+        # Find distances for each similar pairs of vertical boreholes
+        self.borehole_to_borehole_distances_vertical, self.borehole_to_borehole_indices_vertical = \
             self._find_distances(
-                self.boreholes, self.borehole_to_borehole)
+                self.boreholes, self.borehole_to_borehole_vertical)
 
         # Stop chrono
         toc = perf_counter()
@@ -2281,7 +2600,7 @@ def find_similarities(self):
     def _compare_boreholes(self, borehole1, borehole2):
         """
         Compare two boreholes and checks if they have the same dimensions :
-        H, D, and r_b.
+        H, D, r_b, and tilt.
 
         Parameters
         ----------
@@ -2299,16 +2618,17 @@ def _compare_boreholes(self, borehole1, borehole2):
         # Compare lengths (H), buried depth (D) and radius (r_b)
         if (abs((borehole1.H - borehole2.H)/borehole1.H) < self.tol and
             abs((borehole1.r_b - borehole2.r_b)/borehole1.r_b) < self.tol and
-            abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol):
+            abs((borehole1.D - borehole2.D)/(borehole1.D + 1e-30)) < self.tol and
+            abs(abs(borehole1.tilt) - abs(borehole2.tilt))/(abs(borehole1.tilt) + 1e-30) < self.tol):
             similarity = True
         else:
             similarity = False
         return similarity
 
-    def _compare_real_pairs(self, pair1, pair2):
+    def _compare_real_pairs_vertical(self, pair1, pair2):
         """
-        Compare two pairs of boreholes or segments and return True if the two
-        pairs have the same FLS solution for real sources.
+        Compare two pairs of vertical boreholes or segments and return True if
+        the two pairs have the same FLS solution for real sources.
 
         Parameters
         ----------
@@ -2342,10 +2662,10 @@ def _compare_real_pairs(self, pair1, pair2):
             similarity = False
         return similarity
 
-    def _compare_image_pairs(self, pair1, pair2):
+    def _compare_image_pairs_vertical(self, pair1, pair2):
         """
-        Compare two pairs of boreholes or segments and return True if the two
-        pairs have the same FLS solution for mirror sources.
+        Compare two pairs of vertical boreholes or segments and return True if
+        the two pairs have the same FLS solution for mirror sources.
 
         Parameters
         ----------
@@ -2374,10 +2694,11 @@ def _compare_image_pairs(self, pair1, pair2):
             similarity = False
         return similarity
 
-    def _compare_realandimage_pairs(self, pair1, pair2):
+    def _compare_realandimage_pairs_vertical(self, pair1, pair2):
         """
-        Compare two pairs of boreholes or segments and return True if the two
-        pairs have the same FLS solution for both real and mirror sources.
+        Compare two pairs of vertical boreholes or segments and return True if
+        the two pairs have the same FLS solution for both real and mirror
+        sources.
 
         Parameters
         ----------
@@ -2392,8 +2713,159 @@ def _compare_realandimage_pairs(self, pair1, pair2):
             True if the two pairs have the same FLS solution.
 
         """
-        if (self._compare_real_pairs(pair1, pair2)
-            and self._compare_image_pairs(pair1, pair2)):
+        if (self._compare_real_pairs_vertical(pair1, pair2)
+            and self._compare_image_pairs_vertical(pair1, pair2)):
+            similarity = True
+        else:
+            similarity = False
+        return similarity
+
+    def _compare_real_pairs_inclined(self, pair1, pair2):
+        """
+        Compare two pairs of inclined boreholes or segments and return True if
+        the two pairs have the same FLS solution for real sources.
+
+        Parameters
+        ----------
+        pair1 : Tuple of Borehole objects
+            First pair of boreholes or segments.
+        pair2 : Tuple of Borehole objects
+            Second pair of boreholes or segments.
+
+        Returns
+        -------
+        similarity : bool
+            True if the two pairs have the same FLS solution.
+
+        """
+        dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x
+        dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y
+        dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2)
+        theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2)
+        deltaD1 = pair1[0].D - pair1[1].D; deltaD2 = pair2[0].D - pair2[1].D
+        # Equality of lengths between pairs
+        cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol
+            and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol)
+        # Equality of buried depths differences
+        cond_deltaD = abs(deltaD1 - deltaD2)/(abs(deltaD1) + 1e-30) < self.tol
+        # Equality of distances
+        cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol
+        # Equality of tilts
+        cond_beta = (
+            abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol
+            and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol)
+        # Equality of relative orientations
+        sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation)
+        sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation)
+        sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation)
+        sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation)
+        cond_theta = (
+            abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol
+            and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol)
+        if cond_H and cond_deltaD and cond_dis and cond_beta and cond_theta:
+            similarity = True
+        else:
+            similarity = False
+        return similarity
+
+    def _compare_image_pairs_inclined(self, pair1, pair2):
+        """
+        Compare two pairs of inclined boreholes or segments and return True if
+        the two pairs have the same FLS solution for mirror sources.
+
+        Parameters
+        ----------
+        pair1 : Tuple of Borehole objects
+            First pair of boreholes or segments.
+        pair2 : Tuple of Borehole objects
+            Second pair of boreholes or segments.
+
+        Returns
+        -------
+        similarity : bool
+            True if the two pairs have the same FLS solution.
+
+        """
+        dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x
+        dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y
+        dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2)
+        theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2)
+        sumD1 = pair1[0].D + pair1[1].D; sumD2 = pair2[0].D + pair2[1].D
+        # Equality of lengths between pairs
+        cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol
+            and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol)
+        # Equality of buried depths sums
+        cond_sumD = abs(sumD1 - sumD2)/(abs(sumD1) + 1e-30) < self.tol
+        # Equality of distances
+        cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol
+        # Equality of tilts
+        cond_beta = (
+            abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol
+            and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol)
+        # Equality of relative orientations
+        sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation)
+        sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation)
+        sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation)
+        sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation)
+        cond_theta = (
+            abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) > self.tol
+            and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) > self.tol)
+        if cond_H and cond_sumD and cond_dis and cond_beta and cond_theta:
+            similarity = True
+        else:
+            similarity = False
+        return similarity
+
+    def _compare_realandimage_pairs_inclined(self, pair1, pair2):
+        """
+        Compare two pairs of inclined boreholes or segments and return True if
+        the two pairs have the same FLS solution for both real and mirror
+        sources.
+
+        Parameters
+        ----------
+        pair1 : Tuple of Borehole objects
+            First pair of boreholes or segments.
+        pair2 : Tuple of Borehole objects
+            Second pair of boreholes or segments.
+
+        Returns
+        -------
+        similarity : bool
+            True if the two pairs have the same FLS solution.
+
+        Notes
+        -----
+        For inclined boreholes the similarity condition is the same for real
+        and image parts of the solution.
+
+        """
+        dx1 = pair1[0].x - pair1[1].x; dx2 = pair2[0].x - pair2[1].x
+        dy1 = pair1[0].y - pair1[1].y; dy2 = pair2[0].y - pair2[1].y
+        dis1 = np.sqrt(dx1**2 + dy1**2); dis2 = np.sqrt(dx2**2 + dy2**2)
+        theta_12_1 = np.arctan2(dy1, dx1); theta_12_2 = np.arctan2(dy2, dx2)
+        # Equality of lengths between pairs
+        cond_H = (abs((pair1[0].H - pair2[0].H)/pair1[0].H) < self.tol
+            and abs((pair1[1].H - pair2[1].H)/pair1[1].H) < self.tol)
+        # Equality of buried depths
+        cond_D = (
+            abs(pair1[0].D - pair2[0].D)/(abs(pair1[0].D) + 1e-30) < self.tol
+            and abs(pair1[1].D - pair2[1].D)/(abs(pair1[1].D) + 1e-30) < self.tol)
+        # Equality of distances
+        cond_dis = abs(dis1 - dis2)/(abs(dis1) + 1e-30) < self.disTol
+        # Equality of tilts
+        cond_beta = (
+            abs(abs(pair1[0].tilt) - abs(pair2[0].tilt))/(abs(pair1[0].tilt) + 1e-30) < self.tol
+            and abs(abs(pair1[1].tilt) - abs(pair2[1].tilt))/(abs(pair1[1].tilt) + 1e-30) < self.tol)
+        # Equality of relative orientations
+        sin_b1_cos_dt1_1 = np.sin(pair1[0].tilt) * np.cos(theta_12_1 - pair1[0].orientation)
+        sin_b2_cos_dt2_1 = np.sin(pair1[1].tilt) * np.cos(theta_12_1 - pair1[1].orientation)
+        sin_b1_cos_dt1_2 = np.sin(pair2[0].tilt) * np.cos(theta_12_2 - pair2[0].orientation)
+        sin_b2_cos_dt2_2 = np.sin(pair2[1].tilt) * np.cos(theta_12_2 - pair2[1].orientation)
+        cond_theta = (
+            abs(sin_b1_cos_dt1_1 - sin_b1_cos_dt1_2) / (abs(sin_b1_cos_dt1_1) + 1e-30) < self.tol
+            and abs(sin_b2_cos_dt2_1 - sin_b2_cos_dt2_2) / (abs(sin_b2_cos_dt2_1) + 1e-30) < self.tol)
+        if cond_H and cond_D and cond_dis and cond_beta and cond_theta:
             similarity = True
         else:
             similarity = False
@@ -2419,17 +2891,22 @@ def _find_axial_borehole_pairs(self, boreholes):
             boreholes that share the same (pairwise) dimensions (H, D).
 
         """
-        # Compare for the full (real + image) FLS solution
-        compare_pairs = self._compare_realandimage_pairs
-
         nBoreholes = len(boreholes)
-        borehole_to_self = []
+        borehole_to_self_vertical = []
+        borehole_to_self_inclined = []
         # Only check for similarities if there is more than one borehole
         if nBoreholes > 1:
-            borehole_to_borehole = []
+            borehole_to_borehole_vertical = []
+            borehole_to_borehole_inclined = []
             for i, (borehole_i, nSegments_i, ratios_i) in enumerate(
                     zip(boreholes, self.nBoreSegments, self.segment_ratios)):
                 # Compare the borehole to all known unique sets of dimensions
+                if borehole_i.is_vertical():
+                    borehole_to_self = borehole_to_self_vertical
+                    compare_pairs = self._compare_realandimage_pairs_vertical
+                else:
+                    borehole_to_self = borehole_to_self_inclined
+                    compare_pairs = self._compare_realandimage_pairs_inclined
                 for k, borehole_set in enumerate(borehole_to_self):
                     m = borehole_set[0]
                     # Add the borehole to the group if a similar borehole is
@@ -2454,6 +2931,12 @@ def _find_axial_borehole_pairs(self, boreholes):
                     pair0 = (borehole_i, borehole_j) # pair
                     pair1 = (borehole_j, borehole_i) # reciprocal pair
                     # Compare pairs of boreholes to known unique pairs
+                    if borehole_i.is_vertical() and borehole_j.is_vertical():
+                        borehole_to_borehole = borehole_to_borehole_vertical
+                        compare_pairs = self._compare_realandimage_pairs_vertical
+                    else:
+                        borehole_to_borehole = borehole_to_borehole_inclined
+                        compare_pairs = self._compare_realandimage_pairs_inclined
                     for pairs in borehole_to_borehole:
                         m, n = pairs[0]
                         pair_ref = (boreholes[m], boreholes[n])
@@ -2486,11 +2969,19 @@ def _find_axial_borehole_pairs(self, boreholes):
                     # If no similar pairs are known, append the groups
                     else:
                         borehole_to_borehole.append([(i, j)])
+                        
         else:
             # Outputs for a single borehole
-            borehole_to_self = [[0]]
-            borehole_to_borehole = []
-        return borehole_to_self, borehole_to_borehole
+            if boreholes[0].is_vertical:
+                borehole_to_self_vertical = [[0]]
+                borehole_to_self_inclined = []
+            else:
+                borehole_to_self_vertical = []
+                borehole_to_self_inclined = [[0]]
+            borehole_to_borehole_vertical = []
+            borehole_to_borehole_inclined = []
+        return borehole_to_self_vertical, borehole_to_self_inclined, \
+            borehole_to_borehole_vertical, borehole_to_borehole_inclined
 
     def _find_distances(self, boreholes, borehole_to_borehole):
         """
@@ -2559,8 +3050,8 @@ def _find_distances(self, boreholes, borehole_to_borehole):
                 nDis += 1
         return borehole_to_borehole_distances, borehole_to_borehole_indices
 
-    def _map_axial_segment_pairs(self, i, j,
-                                 reaSource=True, imgSource=True):
+    def _map_axial_segment_pairs_vertical(
+            self, i, j, reaSource=True, imgSource=True):
         """
         Find axial (i.e. disregarding the radial distance) similarities between
         segment pairs along two boreholes to simplify the evaluation of the
@@ -2602,13 +3093,13 @@ def _map_axial_segment_pairs(self, i, j,
             "At least one of reaSource and imgSource must be True."
         if reaSource and imgSource:
             # Find segment pairs for the full (real + image) FLS solution
-            compare_pairs = self._compare_realandimage_pairs
+            compare_pairs = self._compare_realandimage_pairs_vertical
         elif reaSource:
             # Find segment pairs for the real FLS solution
-            compare_pairs = self._compare_real_pairs
+            compare_pairs = self._compare_real_pairs_vertical
         elif imgSource:
             # Find segment pairs for the image FLS solution
-            compare_pairs = self._compare_image_pairs
+            compare_pairs = self._compare_image_pairs_vertical
         # Dive both boreholes into segments
         segments1 = borehole1.segments(
             self.nBoreSegments[i], segment_ratios=self.segment_ratios[i])
@@ -2657,8 +3148,143 @@ def _map_axial_segment_pairs(self, i, j,
                 p += 1
         return np.array(H1), np.array(D1), np.array(H2), np.array(D2), i_pair, j_pair, k_pair
 
-    def _map_segment_pairs(self, i_pair, j_pair, k_pair, borehole_to_borehole,
-                           borehole_to_borehole_indices):
+    def _map_axial_segment_pairs_inclined(
+            self, i, j, reaSource=True, imgSource=True):
+        """
+        Find axial similarities between segment pairs along two boreholes to
+        simplify the evaluation of the FLS solution.
+
+        The returned H1, D1, H2, and D2 can be used to evaluate the segment-to-
+        segment response factors using scipy.integrate.quad_vec.
+
+        Parameters
+        ----------
+        i : int
+            Index of the first borehole.
+        j : int
+            Index of the second borehole.
+
+        Returns
+        -------
+        rb1 : array
+            Radii of the emitting heat sources.
+        x1 : array
+            x-Positions of the emitting heat sources.
+        y1 : array
+            y-Positions of the emitting heat sources.
+        H1 : array
+            Lengths of the emitting heat sources.
+        D1 : array
+            Buried depths of the emitting heat sources.
+        tilt1 : array
+            Angles (in radians) from vertical of the emitting heat sources.
+        orientation1 : array
+            Directions (in radians) of the tilt the emitting heat sources.
+        x2 : array
+            x-Positions of the receiving heat sources.
+        y2 : array
+            y-Positions of the receiving heat sources.
+        H2 : array
+            Lengths of the receiving heat sources.
+        D2 : array
+            Buried depths of the receiving heat sources.
+        tilt2 : array
+            Angles (in radians) from vertical of the receiving heat sources.
+        orientation2 : array
+            Directions (in radians) of the tilt the receiving heat sources.
+        i_pair : list
+            Indices of the emitting segments along a borehole.
+        j_pair : list
+            Indices of the receiving segments along a borehole.
+        k_pair : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair).
+        """
+        # Initialize local variables
+        borehole1 = self.boreholes[i]
+        borehole2 = self.boreholes[j]
+        assert reaSource or imgSource, \
+            "At least one of reaSource and imgSource must be True."
+        if reaSource and imgSource:
+            # Find segment pairs for the full (real + image) FLS solution
+            compare_pairs = self._compare_realandimage_pairs_inclined
+        elif reaSource:
+            # Find segment pairs for the real FLS solution
+            compare_pairs = self._compare_real_pairs_inclined
+        elif imgSource:
+            # Find segment pairs for the image FLS solution
+            compare_pairs = self._compare_image_pairs_inclined
+        # Dive both boreholes into segments
+        segments1 = borehole1.segments(
+            self.nBoreSegments[i], segment_ratios=self.segment_ratios[i])
+        segments2 = borehole2.segments(
+            self.nBoreSegments[j], segment_ratios=self.segment_ratios[j])
+        # Prepare lists of FLS-inclined arguments
+        rb1 = []
+        x1 = []
+        y1 = []
+        H1 = []
+        D1 = []
+        tilt1 = []
+        orientation1 = []
+        x2 = []
+        y2 = []
+        H2 = []
+        D2 = []
+        tilt2 = []
+        orientation2 = []
+        # All possible pairs (i, j) of indices between segments
+        i_pair = np.repeat(np.arange(self.nBoreSegments[i], dtype=np.uint),
+                           self.nBoreSegments[j])
+        j_pair = np.tile(np.arange(self.nBoreSegments[j], dtype=np.uint),
+                         self.nBoreSegments[i])
+        # Empty list of indices for unique pairs
+        k_pair = np.empty(self.nBoreSegments[i] * self.nBoreSegments[j],
+                          dtype=np.uint)
+        unique_pairs = []
+        nPairs = 0
+
+        p = 0
+        for ii, segment_i in enumerate(segments1):
+            for jj, segment_j in enumerate(segments2):
+                pair = (segment_i, segment_j)
+                # Compare the segment pairs to all known unique pairs
+                for k, pair_k in enumerate(unique_pairs):
+                    m, n = pair_k[0], pair_k[1]
+                    pair_ref = (segments1[m], segments2[n])
+                    # Stop if a similar pair is found and assign the index
+                    if compare_pairs(pair, pair_ref):
+                        k_pair[p] = k
+                        break
+                # If no similar pair is found : add a new pair, increment the
+                # number of unique pairs, and extract the associated buried
+                # depths
+                else:
+                    k_pair[p] = nPairs
+                    rb1.append(segment_i.r_b)
+                    x1.append(segment_i.x)
+                    y1.append(segment_i.y)
+                    H1.append(segment_i.H)
+                    D1.append(segment_i.D)
+                    tilt1.append(segment_i.tilt)
+                    orientation1.append(segment_i.orientation)
+                    x2.append(segment_j.x)
+                    y2.append(segment_j.y)
+                    H2.append(segment_j.H)
+                    D2.append(segment_j.D)
+                    tilt2.append(segment_j.tilt)
+                    orientation2.append(segment_j.orientation)
+                    unique_pairs.append((ii, jj))
+                    nPairs += 1
+                p += 1
+        return np.array(rb1), np.array(x1), np.array(y1), np.array(H1), \
+            np.array(D1), np.array(tilt1), np.array(orientation1), \
+            np.array(x2), np.array(y2), np.array(H2), np.array(D2), \
+            np.array(tilt2), np.array(orientation2), i_pair, j_pair, k_pair
+
+    def _map_segment_pairs_vertical(
+            self, i_pair, j_pair, k_pair, borehole_to_borehole,
+            borehole_to_borehole_indices):
         """
         Return the maping of the unique segment-to-segment thermal response
         factors (h) to the complete h_ij array of the borefield, such that:
@@ -2705,6 +3331,47 @@ def _map_segment_pairs(self, i_pair, j_pair, k_pair, borehole_to_borehole,
             [np.repeat(i, len(k_pair)) for i in borehole_to_borehole_indices])
         return i_segment, j_segment, k_segment, l_segment
 
+    def _map_segment_pairs_inclined(
+            self, i_pair, j_pair, k_pair, borehole_to_borehole):
+        """
+        Return the maping of the unique segment-to-segment thermal response
+        factors (h) to the complete h_ij array of the borefield, such that:
+
+            h_ij[j_segment, i_segment, :nt] = h[:nt, k_segment].T,
+
+        where h is the array of unique segment-to-segment thermal response
+        factors for a given unique pair of boreholes at all unique distances.
+
+        Parameters
+        ----------
+        i_pair : list
+            Indices of the emitting segments.
+        j_pair : list
+            Indices of the receiving segments.
+        k_pair : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair).
+        borehole_to_borehole : list
+            Tuples of borehole indexes.
+
+        Returns
+        -------
+        i_segment : list
+            Indices of the emitting segments in the bore field.
+        j_segment : list
+            Indices of the receiving segments in the bore field.
+        k_segment : list
+            Indices of unique segment pairs in the (H1, D1, H2, D2) dimensions
+            corresponding to all pairs in (i_pair, j_pair) in the bore field.
+
+        """
+        i_segment = np.concatenate(
+            [i_pair + self._i0Segments[i] for (i, j) in borehole_to_borehole])
+        j_segment = np.concatenate(
+            [j_pair + self._i0Segments[j] for (i, j) in borehole_to_borehole])
+        k_segment = np.tile(k_pair, len(borehole_to_borehole))
+        return i_segment, j_segment, k_segment
+
     def _check_solver_specific_inputs(self):
         """
         This method ensures that solver specific inputs to the Solver object
@@ -2782,6 +3449,10 @@ class _Equivalent(_BaseSolver):
         Number of terms in the approximation of the FLS solution. This
         parameter is unused if `approximate_FLS` is set to False.
         Default is 10. Maximum is 25.
+    mQuad : int, optional
+        Number of Gauss-Legendre sample points for the integral over :math:`u`
+        in the inclined FLS solution.
+        Default is 11.
     disp : bool, optional
         Set to true to print progression messages.
         Default is False.
@@ -3487,6 +4158,8 @@ def _check_solver_specific_inputs(self):
             "Solver 'equivalent' can only handle equal numbers of segments."
         assert np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios]), \
             "Solver 'equivalent' can only handle identical segment_ratios for all boreholes."
+        assert not np.any([b.is_tilted() for b in self.boreholes]), \
+            "Solver 'equivalent' can only handle vertical boreholes."
         if self.boundary_condition == 'MIFT':
             assert np.all(np.array(self.network.c, dtype=int) == -1), \
                 "Solver 'equivalent' is only valid for parallel-connected " \
diff --git a/pygfunction/heat_transfer.py b/pygfunction/heat_transfer.py
index 14aada3a..323b3e08 100644
--- a/pygfunction/heat_transfer.py
+++ b/pygfunction/heat_transfer.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 from scipy.integrate import quad, quad_vec
-from scipy.special import erfc
+from scipy.special import erfc, erf, roots_legendre
 
 from .boreholes import Borehole
 from .utilities import erfint, exp1, _erf_coeffs
@@ -9,14 +9,15 @@
 
 def finite_line_source(
         time, alpha, borehole1, borehole2, reaSource=True, imgSource=True,
-        approximation=False, N=10):
+        approximation=False, M=11, N=10):
     """
     Evaluate the Finite Line Source (FLS) solution.
 
     This function uses a numerical quadrature to evaluate the one-integral form
-    of the FLS solution, as proposed by Claesson and Javed [#FLS-ClaJav2011]_
-    and extended to boreholes with different vertical positions by Cimmino and
-    Bernier [#FLS-CimBer2014]_. The FlS solution is given by:
+    of the FLS solution. For vertical boreholes, the FLS solution was proposed
+    by Claesson and Javed [#FLS-ClaJav2011]_ and extended to boreholes with
+    different vertical positions by Cimmino and Bernier [#FLS-CimBer2014]_.
+    The FlS solution is given by:
 
         .. math::
             h_{1\\rightarrow2}(t) &= \\frac{1}{2H_2}
@@ -24,6 +25,9 @@ def finite_line_source(
             e^{-d_{12}^2s^2}(I_{real}(s)+I_{imag}(s))ds
 
 
+            d_{12} &= \\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}
+
+
             I_{real}(s) &= erfint((D_2-D_1+H_2)s) - erfint((D_2-D_1)s)
 
             &+ erfint((D_2-D_1-H_1)s) - erfint((D_2-D_1+H_2-H_1)s)
@@ -37,6 +41,73 @@ def finite_line_source(
 
                       &= Xerf(X) - \\frac{1}{\\sqrt{\\pi}}(1-e^{-X^2})
 
+    For inclined boreholes, the FLS solution was proposed by Lazzarotto
+    [#FLS-Lazzar2016]_ and Lazzarotto and Björk [#FLS-LazBjo2016]_.
+    The FLS solution is given by:
+
+        .. math::
+            h_{1\\rightarrow2}(t) &= \\frac{H_1}{2H_2}
+            \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty}
+            \\frac{1}{s}
+            \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds
+
+
+            I_{real}(u, s) &=
+            e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 - D_2)^2) s^2}
+
+            &\\cdot (erf((u H_1 k_{0,real} + k_{2,real}) s)
+             - erf((u H_1 k_{0,real} + k_{2,real} - H_2) s))
+
+            &\\cdot  e^{(u^2 H_1^2 (k_{0,real}^2 - 1)
+                + 2 u H_1 (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2}
+            du ds
+
+
+            I_{imag}(u, s) &=
+            -e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 + D_2)^2) s^2}
+
+            &\\cdot (erf((u H_1 k_{0,imag} + k_{2,imag}) s)
+             - erf((u H_1 k_{0,imag} + k_{2,imag} - H_2) s))
+
+            &\\cdot e^{(u^2 H_1^2 (k_{0,imag}^2 - 1)
+                + 2 u H_1 (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2}
+            du ds
+
+
+            k_{0,real} &=
+            sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2)
+            + cos(\\beta_1) cos(\\beta_2)
+
+
+            k_{0,imag} &=
+            sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2)
+            - cos(\\beta_1) cos(\\beta_2)
+
+
+            k_{1,real} &= sin(\\beta_1)
+            (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2))
+            + cos(\\beta_1) (D_1 - D_2)
+
+
+            k_{1,imag} &= sin(\\beta_1)
+            (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2))
+            + cos(\\beta_1) (D_1 + D_2)
+
+
+            k_{2,real} &= sin(\\beta_2)
+            (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2))
+            + cos(\\beta_2) (D_1 - D_2)
+
+
+            k_{2,imag} &= sin(\\beta_2)
+            (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2))
+            - cos(\\beta_2) (D_1 + D_2)
+
+    where :math:`\\beta_1` and :math:`\\beta_2` are the tilt angle of the
+    boreholes (relative to vertical), and :math:`\\theta_1` and
+    :math:`\\theta_2` are the orientation of the boreholes (relative to the
+    x-axis).
+
         .. Note::
             The reciprocal thermal response factor
             :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by:
@@ -66,6 +137,10 @@ def finite_line_source(
         (2021) [#FLS-Cimmin2021]_. This approximation does not require
         the numerical evaluation of any integral.
         Default is False.
+    M : int, optional
+        Number of Gauss-Legendre sample points for the quadrature over
+        :math:`u`. This is only used for inclined boreholes.
+        Default is 11.
     N : int, optional
         Number of terms in the approximation of the FLS solution. This
         parameter is unused if `approximation` is set to False.
@@ -97,6 +172,11 @@ def finite_line_source(
     >>> h = gt.heat_transfer.finite_line_source(
         4*168*3600., 1.0e-6, b1, b2, approximation=True, N=10)
     h = 0.0110474667731
+    >>> b3 = gt.boreholes.Borehole(
+        H=150., D=4., r_b=0.075, x=5., y=0., tilt=3.1415/15, orientation=0.)
+    >>> h = gt.heat_transfer.finite_line_source(
+        4*168*3600., 1.0e-6, b1, b3, M=21)
+    h = 0.0002017450051
 
     References
     ----------
@@ -110,44 +190,94 @@ def finite_line_source(
        finite line source solution to model thermal interactions between
        geothermal boreholes. International Communications in Heat and Mass
        Transfer, 127, 105496.
+    .. [#FLS-Lazzar2016] Lazzarotto, A. (2016). A methodology for the
+       calculation of response functions for geothermal fields with
+       arbitrarily oriented boreholes – Part 1, Renewable Energy, 86,
+       1380-1393.
+    .. [#FLS-LazBjo2016] Lazzarotto, A., & Björk, F. (2016). A methodology for
+       the calculation of response functions for geothermal fields with
+       arbitrarily oriented boreholes – Part 2, Renewable Energy, 86,
+       1353-1361.
 
     """
     if isinstance(borehole1, Borehole) and isinstance(borehole2, Borehole):
         # Unpack parameters
-        dis = borehole1.distance(borehole2)
         H1, D1 = borehole1.H, borehole1.D
         H2, D2 = borehole2.H, borehole2.D
-        # Integrand of the finite line source solution
-        f = _finite_line_source_integrand(
-            dis, H1, D1, H2, D2, reaSource, imgSource)
-
-        # Evaluate integral
-        if time is np.inf:
-            h = _finite_line_source_steady_state(
-                dis, H1, D1, H2, D2, reaSource, imgSource)
-        elif isinstance(time, (np.floating, float)):
-            if not approximation:
-                # Lower bound of integration
-                a = 1.0 / np.sqrt(4.0*alpha*time)
-                h = 0.5 / H2 * quad(f, a, np.inf)[0]
-            else:
+        if borehole1.is_vertical() and borehole2.is_vertical():
+            # Boreholes are vertical
+            dis = borehole1.distance(borehole2)
+            if time is np.inf:
+                # Steady-state solution
+                h = _finite_line_source_steady_state(
+                    dis, H1, D1, H2, D2, reaSource, imgSource)
+            elif approximation:
+                # Approximation
                 h = finite_line_source_approximation(
                     time, alpha, dis, H1, D1, H2, D2,
                     reaSource=reaSource, imgSource=imgSource, N=N)
+            else:
+                # Integrand of the finite line source solution
+                f = _finite_line_source_integrand(
+                    dis, H1, D1, H2, D2, reaSource, imgSource)
+                # Evaluate integral
+                if isinstance(time, (np.floating, float)):
+                    # Lower bound of integration
+                    a = 1.0 / np.sqrt(4.0*alpha*time)
+                    h = 0.5 / H2 * quad(f, a, np.inf)[0]
+                else:
+                    # Lower bound of integration
+                    a = 1.0 / np.sqrt(4.0*alpha*time)
+                    # Upper bound of integration
+                    b = np.concatenate(([np.inf], a[:-1]))
+                    h = np.cumsum(np.stack(
+                        [0.5 / H2 * quad(f, a_i, b_i)[0]
+                         for t, a_i, b_i in zip(time, a, b)],
+                        axis=-1), axis=-1)
         else:
-            if not approximation:
-                # Lower bound of integration
-                a = 1.0 / np.sqrt(4.0*alpha*time)
-                # Upper bound of integration
-                b = np.concatenate(([np.inf], a[:-1]))
-                h = np.cumsum(np.stack(
-                    [0.5 / H2 * quad(f, a_i, b_i)[0]
-                     for t, a_i, b_i in zip(time, a, b)],
-                    axis=-1), axis=-1)
+            # At least one borehole is tilted
+            # Unpack parameters
+            x1, y1 = borehole1.position()
+            rb1 = borehole1.r_b
+            tilt1 = borehole1.tilt
+            orientation1 = borehole1.orientation
+            x2, y2 = borehole2.position()
+            tilt2 = borehole2.tilt
+            orientation2 = borehole2.orientation
+            if time is np.inf:
+                # Steady-state solution
+                h = _finite_line_source_inclined_steady_state(
+                    rb1, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource, imgSource, M=M)
+            elif approximation:
+                # Approximation
+                h = finite_line_source_inclined_approximation(
+                    time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource=reaSource, imgSource=imgSource, M=M, N=N)
             else:
-                h = finite_line_source_approximation(
-                    time, alpha, dis, H1, D1, H2, D2,
-                    reaSource=reaSource, imgSource=imgSource, N=N)
+                # Integrand of the inclined finite line source solution
+                f = _finite_line_source_inclined_integrand(
+                    rb1, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource, imgSource, M)
+
+                # Evaluate integral
+                if isinstance(time, (np.floating, float)):
+                    # Lower bound of integration
+                    a = 1.0 / np.sqrt(4.0*alpha*time)
+                    h = 0.5 / H2 * quad(f, a, np.inf)[0]
+                else:
+                    # Lower bound of integration
+                    a = 1.0 / np.sqrt(4.0*alpha*time)
+                    # Upper bound of integration
+                    b = np.concatenate(([np.inf], a[:-1]))
+                    h = np.cumsum(np.stack(
+                        [0.5 / H2 * quad(f, a_i, b_i)[0]
+                         for t, a_i, b_i in zip(time, a, b)],
+                        axis=-1), axis=-1)
+
     else:
         # Unpack parameters
         if isinstance(borehole1, Borehole): borehole1 = [borehole1]
@@ -157,25 +287,62 @@ def finite_line_source(
         x2 = np.array([b.x for b in borehole2])
         y2 = np.array([b.y for b in borehole2])
         r_b = np.array([b.r_b for b in borehole1])
-        dis = np.maximum(np.sqrt(np.add.outer(x2, -x1)**2 + np.add.outer(y2, -y1)**2), r_b)
+        dis = np.maximum(
+            np.sqrt(np.add.outer(x2, -x1)**2 + np.add.outer(y2, -y1)**2),
+            r_b)
         D1 = np.array([b.D for b in borehole1]).reshape(1, -1)
         H1 = np.array([b.H for b in borehole1]).reshape(1, -1)
         D2 = np.array([b.D for b in borehole2]).reshape(-1, 1)
         H2 = np.array([b.H for b in borehole2]).reshape(-1, 1)
 
-        if time is np.inf:
-            h = _finite_line_source_steady_state(
-                dis, H1, D1, H2, D2, reaSource, imgSource)
-        else:
-            if not approximation:
+        if (np.all([b.is_vertical() for b in borehole1])
+            and np.all([b.is_vertical() for b in borehole2])):
+            # All boreholes are vertical
+            if time is np.inf:
+                # Steady-state solution
+                h = _finite_line_source_steady_state(
+                    dis, H1, D1, H2, D2, reaSource, imgSource)
+            elif approximation:
+                # Approximation
+                h = finite_line_source_approximation(
+                    time, alpha, dis, H1, D1, H2, D2,
+                    reaSource=reaSource, imgSource=imgSource, N=N)
+            else:
                 # Evaluate integral
                 h = finite_line_source_vectorized(
                     time, alpha, dis, H1, D1, H2, D2,
                     reaSource=reaSource, imgSource=imgSource)
+        else:
+            # At least one borehole is tilted
+            # Unpack parameters
+            x1 = x1.reshape(1, -1)
+            y1 = y1.reshape(1, -1)
+            tilt1 = np.array([b.tilt for b in borehole1]).reshape(1, -1)
+            orientation1 = np.array([b.orientation for b in borehole1]).reshape(1, -1)
+            x2 = x2.reshape(-1, 1)
+            y2 = y2.reshape(-1, 1)
+            tilt2 = np.array([b.tilt for b in borehole2]).reshape(-1, 1)
+            orientation2 = np.array([b.orientation for b in borehole2]).reshape(-1, 1)
+            r_b = r_b.reshape(1, -1)
+            if time is np.inf:
+                # Steady-state solution
+                h = _finite_line_source_inclined_steady_state(
+                    r_b, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource, imgSource, M=M)
+            elif approximation:
+                # Approximation
+                h = finite_line_source_inclined_approximation(
+                    time, alpha, r_b, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource=reaSource, imgSource=imgSource, M=M, N=N)
             else:
-                h = finite_line_source_approximation(
-                    time, alpha, dis, H1, D1, H2, D2,
-                    reaSource=reaSource, imgSource=imgSource, N=N)
+                # Evaluate integral
+                h = finite_line_source_inclined_vectorized(
+                    time, alpha,
+                    r_b, x1, y1, H1, D1, tilt1, orientation1,
+                    x2, y2, H2, D2, tilt2, orientation2,
+                    reaSource=reaSource, imgSource=imgSource, M=M)
     return h
 
 
@@ -202,10 +369,10 @@ def finite_line_source_approximation(
         Lengths of the receiving heat sources.
     D2 : float or array
         Buried depths of the receiving heat sources.
-    reaSource : bool
+    reaSource : bool, optional
         True if the real part of the FLS solution is to be included.
         Default is True.
-    imgSource : bool
+    imgSource : bool, optional
         True if the image part of the FLS solution is to be included.
         Default is True.
     N : int, optional
@@ -288,6 +455,201 @@ def finite_line_source_approximation(
     return h
 
 
+def finite_line_source_inclined_approximation(
+        time, alpha, 
+        rb1, x1, y1, H1, D1, tilt1, orientation1,
+        x2, y2, H2, D2, tilt2, orientation2,
+        reaSource=True, imgSource=True, M=11, N=10):
+    """
+    Evaluate the inclined Finite Line Source (FLS) solution using the
+    approximation method of Cimmino (2021) [#IncFLSApprox-Cimmin2021]_.
+
+    Parameters
+    ----------
+    time : float or array, shape (K)
+        Value of time (in seconds) for which the FLS solution is evaluated.
+    alpha : float
+        Soil thermal diffusivity (in m2/s).
+    rb1 : array
+        Radii of the emitting heat sources.
+    x1 : float or array
+        x-Positions of the emitting heat sources.
+    y1 : float or array
+        y-Positions of the emitting heat sources.
+    H1 : float or array
+        Lengths of the emitting heat sources.
+    D1 : float or array
+        Buried depths of the emitting heat sources.
+    tilt1 : float or array
+        Angles (in radians) from vertical of the emitting heat sources.
+    orientation1 : float or array
+        Directions (in radians) of the tilt the emitting heat sources.
+    x2 : array
+        x-Positions of the receiving heat sources.
+    y2 : array
+        y-Positions of the receiving heat sources.
+    H2 : float or array
+        Lengths of the receiving heat sources.
+    D2 : float or array
+        Buried depths of the receiving heat sources.
+    tilt2 : float or array
+        Angles (in radians) from vertical of the receiving heat sources.
+    orientation2 : float or array
+        Directions (in radians) of the tilt the receiving heat sources.
+    reaSource : bool, optional
+        True if the real part of the FLS solution is to be included.
+        Default is True.
+    imgSource : bool, optional
+        True if the image part of the FLS solution is to be included.
+        Default is true.
+    M : int, optional
+        Number of points for the Gauss-Legendre quadrature rule along the
+        receiving heat sources.
+        Default is 21.
+    N : int, optional
+        Number of terms in the approximation of the FLS solution.
+        Default is 10. Maximum is 25.
+
+    Returns
+    -------
+    h : float
+        Value of the FLS solution. The average (over the length) temperature
+        drop on the wall of borehole2 due to heat extracted from borehole1 is:
+
+        .. math:: \\Delta T_{b,2} = T_g - \\frac{Q_1}{2\\pi k_s H_2} h
+
+    References
+    ----------
+    .. [#IncFLSApprox-Cimmin2021] Cimmino, M. (2021). An approximation of the
+       finite line source solution to model thermal interactions between
+       geothermal boreholes. International Communications in Heat and Mass
+       Transfer, 127, 105496.
+
+    """
+    # Expected output shape of h, excluding time
+    output_shape = np.broadcast_shapes(
+            *[np.shape(arg) for arg in (
+                rb1, x1, y1, H1, D1, tilt1, orientation1,
+                x2, y2, H2, D2, tilt2, orientation2)])
+    # Number of dimensions of the output, excluding time
+    ouput_ndim = len(output_shape)
+    # Shape of the time variable
+    time_shape = np.shape(time)
+    # Number of dimensions of the time variable
+    time_ndim = len(time_shape)
+    # Roots for Gauss-Legendre quadrature
+    x, w = roots_legendre(M)
+    u = (0.5 * x + 0.5).reshape((-1, 1) + (1,) * ouput_ndim)
+    w = w / 2
+    # Coefficients of the approximation of the error function
+    a, b = _erf_coeffs(N)
+    b = b.reshape((1, -1) + (1,) * ouput_ndim)
+    # Sines and cosines of tilt (b: beta) and orientation (t: theta)
+    sb1 = np.sin(tilt1)
+    sb2 = np.sin(tilt2)
+    cb1 = np.cos(tilt1)
+    cb2 = np.cos(tilt2)
+    st1 = np.sin(orientation1)
+    st2 = np.sin(orientation2)
+    ct1 = np.cos(orientation1)
+    ct2 = np.cos(orientation2)
+    ct12 = np.cos(orientation1 - orientation2)
+    # Horizontal distances
+    dx = x1 - x2
+    dy = y1 - y2
+    rr = dx**2 + dy**2  # Squared radial distance
+    # Length ratios
+    H_ratio = H1 / H2
+    H_ratio = np.reshape(H_ratio, np.shape(H_ratio) + (1,) * time_ndim)
+    # Approximation
+    ss = 1. / (4 * alpha * time)
+    if reaSource and imgSource:
+        # Full (real + image) FLS solution
+        # Axial distances
+        dzRea = D1 - D2
+        dzImg = D1 + D2
+        # FLS-inclined coefficients
+        kRea_0 = sb1 * sb2 * ct12 + cb1 * cb2
+        kImg_0 = sb1 * sb2 * ct12 - cb1 * cb2
+        kRea_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzRea
+        kImg_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzImg
+        kRea_2 = sb2 * (ct2 * dx + st2 * dy) + cb2 * dzRea
+        kImg_2 = sb2 * (ct2 * dx + st2 * dy) - cb2 * dzImg
+        dRea_1 = u * H1 * kRea_0 + kRea_2
+        dImg_1 = u * H1 * kImg_0 + kImg_2
+        dRea_2 = u * H1 * kRea_0 + kRea_2 - H2
+        dImg_2 = u * H1 * kImg_0 + kImg_2 - H2
+        cRea = np.maximum(
+            rr + dzRea**2 - dRea_1**2 + (kRea_1 + u * H1)**2 - kRea_1**2,
+            rb1**2)
+        cImg = np.maximum(
+            rr + dzImg**2 - dImg_1**2 + (kImg_1 + u * H1)**2 - kImg_1**2,
+            rb1**2)
+        # Signs for summation
+        pRea_1 = np.sign(dRea_1)
+        pRea_1 = np.reshape(pRea_1, np.shape(pRea_1) + (1,) * time_ndim)
+        pRea_2 = np.sign(dRea_2)
+        pRea_2 = np.reshape(pRea_2, np.shape(pRea_2) + (1,) * time_ndim)
+        pImg_1 = np.sign(dImg_1)
+        pImg_1 = np.reshape(pImg_1, np.shape(pImg_1) + (1,) * time_ndim)
+        pImg_2 = np.sign(dImg_2)
+        pImg_2 = np.reshape(pImg_2, np.shape(pImg_2) + (1,) * time_ndim)
+        # FLS-inclined approximation
+        h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a,
+            (pRea_1 * exp1(np.multiply.outer(cRea + b * dRea_1**2, ss)) \
+            - pRea_2 * exp1(np.multiply.outer(cRea + b * dRea_2**2, ss)) \
+            - pImg_1 * exp1(np.multiply.outer(cImg + b * dImg_1**2, ss)) \
+            + pImg_2 * exp1(np.multiply.outer(cImg + b * dImg_2**2, ss))) )
+    elif reaSource:
+        # Real FLS solution
+        # Axial distance
+        dzRea = D1 - D2
+        # FLS-inclined coefficients
+        kRea_0 = sb1 * sb2 * ct12 + cb1 * cb2
+        kRea_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzRea
+        kRea_2 = sb2 * (ct2 * dx + st2 * dy) + cb2 * dzRea
+        dRea_1 = u * H1 * kRea_0 + kRea_2
+        dRea_2 = u * H1 * kRea_0 + kRea_2 - H2
+        cRea = np.maximum(
+            rr + dzRea**2 - dRea_1**2 + (kRea_1 + u * H1)**2 - kRea_1**2,
+            rb1**2)
+        # Signs for summation
+        pRea_1 = np.sign(dRea_1)
+        pRea_1 = np.reshape(pRea_1, np.shape(pRea_1) + (1,) * time_ndim)
+        pRea_2 = np.sign(dRea_2)
+        pRea_2 = np.reshape(pRea_2, np.shape(pRea_2) + (1,) * time_ndim)
+        # FLS-inclined approximation
+        h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a,
+            (pRea_1 * exp1(np.multiply.outer(cRea + b * dRea_1**2, ss)) \
+            - pRea_2 * exp1(np.multiply.outer(cRea + b * dRea_2**2, ss))) )
+    elif imgSource:
+        # Image FLS solution
+        # Axial distance
+        dzImg = D1 + D2
+        # FLS-inclined coefficients
+        kImg_0 = sb1 * sb2 * ct12 - cb1 * cb2
+        kImg_1 = sb1 * (ct1 * dx + st1 * dy) + cb1 * dzImg
+        kImg_2 = sb2 * (ct2 * dx + st2 * dy) - cb2 * dzImg
+        dImg_1 = u * H1 * kImg_0 + kImg_2
+        dImg_2 = u * H1 * kImg_0 + kImg_2 - H2
+        cImg = np.maximum(
+            rr + dzImg**2 - dImg_1**2 + (kImg_1 + u * H1)**2 - kImg_1**2,
+            rb1**2)
+        # Signs for summation
+        pImg_1 = np.sign(dImg_1)
+        pImg_1 = np.reshape(pImg_1, np.shape(pImg_1) + (1,) * time_ndim)
+        pImg_2 = np.sign(dImg_2)
+        pImg_2 = np.reshape(pImg_2, np.shape(pImg_2) + (1,) * time_ndim)
+        # FLS-inclined approximation
+        h = 0.25 * H_ratio * np.einsum('i,j,ij...', w, a,
+            (-pImg_1 * exp1(np.multiply.outer(cImg + b * dImg_1**2, ss)) \
+            + pImg_2 * exp1(np.multiply.outer(cImg + b * dImg_2**2, ss))) )
+    else:
+        # No heat source
+        h = np.zeros(output_shape + np.shape(time))
+    return h
+
+
 def finite_line_source_vectorized(
         time, alpha, dis, H1, D1, H2, D2, reaSource=True, imgSource=True,
         approximation=False, N=10):
@@ -526,6 +888,195 @@ def finite_line_source_equivalent_boreholes_vectorized(
     return h
 
 
+def finite_line_source_inclined_vectorized(
+        time, alpha,
+        rb1, x1, y1, H1, D1, tilt1, orientation1,
+        x2, y2, H2, D2, tilt2, orientation2,
+        reaSource=True, imgSource=True, M=11, approximation=False, N=10):
+    """
+    Evaluate the inclined Finite Line Source (FLS) solution.
+
+    This function uses a numerical quadrature to evaluate the inclined FLS
+    solution, as proposed by Lazzarotto [#incFLSVec-Lazzar2016]_.
+    The inclined FLS solution is given by:
+
+        .. math::
+            h_{1\\rightarrow2}(t) &= \\frac{H_1}{2H_2}
+            \\int_{\\frac{1}{\\sqrt{4\\alpha t}}}^{\\infty}
+            \\frac{1}{s}
+            \\int_{0}^{1} (I_{real}(u, s)+I_{imag}(u, s)) du ds
+
+
+            I_{real}(u, s) &=
+            e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 - D_2)^2) s^2}
+
+            &\\cdot (erf((u H_1 k_{0,real} + k_{2,real}) s)
+             - erf((u H_1 k_{0,real} + k_{2,real} - H_2) s))
+
+            &\\cdot  e^{(u^2 H_1^2 (k_{0,real}^2 - 1)
+                + 2 u H_1 (k_{0,real} k_{2,real} - k_{1,real}) + k_{2,real}^2) s^2}
+            du ds
+
+
+            I_{imag}(u, s) &=
+            -e^{-((x_1 - x_2)^2 + (y_1 - y_2)^2 + (D_1 + D_2)^2) s^2}
+
+            &\\cdot (erf((u H_1 k_{0,imag} + k_{2,imag}) s)
+             - erf((u H_1 k_{0,imag} + k_{2,imag} - H_2) s))
+
+            &\\cdot e^{(u^2 H_1^2 (k_{0,imag}^2 - 1)
+                + 2 u H_1 (k_{0,imag} k_{2,imag} - k_1) + k_{2,imag}^2) s^2}
+            du ds
+
+
+            k_{0,real} &=
+            sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2)
+            + cos(\\beta_1) cos(\\beta_2)
+
+
+            k_{0,imag} &=
+            sin(\\beta_1) sin(\\beta_2) cos(\\theta_1 - \\theta_2)
+            - cos(\\beta_1) cos(\\beta_2)
+
+
+            k_{1,real} &= sin(\\beta_1)
+            (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2))
+            + cos(\\beta_1) (D_1 - D_2)
+
+
+            k_{1,imag} &= sin(\\beta_1)
+            (cos(\\theta_1) (x_1 - x_2) + sin(\\theta_1) (y_1 - y_2))
+            + cos(\\beta_1) (D_1 + D_2)
+
+
+            k_{2,real} &= sin(\\beta_2)
+            (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2))
+            + cos(\\beta_2) (D_1 - D_2)
+
+
+            k_{2,imag} &= sin(\\beta_2)
+            (cos(\\theta_2) (x_1 - x_2) + sin(\\theta_2) (y_1 - y_2))
+            - cos(\\beta_2) (D_1 + D_2)
+
+    where :math:`\\beta_1` and :math:`\\beta_2` are the tilt angle of the
+    boreholes (relative to vertical), and :math:`\\theta_1` and
+    :math:`\\theta_2` are the orientation of the boreholes (relative to the
+    x-axis).
+
+        .. Note::
+            The reciprocal thermal response factor
+            :math:`h_{2\\rightarrow1}(t)` can be conveniently calculated by:
+
+                .. math::
+                    h_{2\\rightarrow1}(t) = \\frac{H_2}{H_1}
+                    h_{1\\rightarrow2}(t)
+
+    Parameters
+    ----------
+    time : float or array, shape (K)
+        Value of time (in seconds) for which the FLS solution is evaluated.
+    alpha : float
+        Soil thermal diffusivity (in m2/s).
+    rb1 : array
+        Radii of the emitting heat sources.
+    x1 : float or array
+        x-Positions of the emitting heat sources.
+    y1 : float or array
+        y-Positions of the emitting heat sources.
+    H1 : float or array
+        Lengths of the emitting heat sources.
+    D1 : float or array
+        Buried depths of the emitting heat sources.
+    tilt1 : float or array
+        Angles (in radians) from vertical of the emitting heat sources.
+    orientation1 : float or array
+        Directions (in radians) of the tilt the emitting heat sources.
+    x2 : array
+        x-Positions of the receiving heat sources.
+    y2 : array
+        y-Positions of the receiving heat sources.
+    H2 : float or array
+        Lengths of the receiving heat sources.
+    D2 : float or array
+        Buried depths of the receiving heat sources.
+    tilt2 : float or array
+        Angles (in radians) from vertical of the receiving heat sources.
+    orientation2 : float or array
+        Directions (in radians) of the tilt the receiving heat sources.
+    reaSource : bool, optional
+        True if the real part of the FLS solution is to be included.
+        Default is True.
+    imgSource : bool, optional
+        True if the image part of the FLS solution is to be included.
+        Default is true.
+    M : int, optional
+        Number of points for the Gauss-Legendre quadrature rule along the
+        receiving heat sources.
+        Default is 21.
+    approximation : bool, optional
+        Set to true to use the approximation of the FLS solution of Cimmino
+        (2021) [#FLSVec-Cimmin2021]_. This approximation does not require
+        the numerical evaluation of any integral.
+        Default is False.
+    N : int, optional
+        Number of terms in the approximation of the FLS solution. This
+        parameter is unused if `approximation` is set to False.
+        Default is 10. Maximum is 25.
+
+    Returns
+    -------
+    f : callable
+        Integrand of the finite line source solution. Can be vector-valued.
+
+    Notes
+    -----
+    This is a vectorized version of the :func:`finite_line_source` function
+    using scipy.integrate.quad_vec to speed up calculations. All arrays
+    (x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2,
+    orientation2) must follow numpy array broadcasting rules.
+
+    References
+    ----------
+    .. [#incFLSVec-Lazzar2016] Lazzarotto, A. (2016). A methodology for the
+       calculation of response functions for geothermal fields with
+       arbitrarily oriented boreholes – Part 1, Renewable Energy, 86,
+       1380-1393.
+
+    """
+    if not approximation:
+        # Integrand of the inclined finite line source solution
+        f = _finite_line_source_inclined_integrand(
+            rb1, x1, y1, H1, D1, tilt1, orientation1,
+            x2, y2, H2, D2, tilt2, orientation2,
+            reaSource, imgSource, M)
+    
+        # Evaluate integral
+        if isinstance(time, (np.floating, float)):
+            # Lower bound of integration
+            a = 1.0 / np.sqrt(4.0*alpha*time)
+            h = 0.5 / H2 * quad_vec(f, a, np.inf, epsabs=1e-4, epsrel=1e-6)[0]
+        else:
+            # The real and image parts are split to avoid overflow in the
+            # integrand function
+            # Lower bound of integration
+            a = 1.0 / np.sqrt(4.0*alpha*time)
+            # Upper bound of integration
+            b = np.concatenate(([np.inf], a[:-1]))
+            h = np.cumsum(
+                np.stack(
+                    [0.5 / H2 * quad_vec(
+                        f, a_i, b_i, epsabs=1e-4, epsrel=1e-6)[0]
+                     for i, (a_i, b_i) in enumerate(zip(a, b))],
+                    axis=-1),
+                axis=-1)
+    else:
+        h = finite_line_source_inclined_approximation(
+            time, alpha, rb1, x1, y1, H1, D1, tilt1, orientation1,
+            x2, y2, H2, D2, tilt2, orientation2,
+            reaSource=reaSource, imgSource=imgSource, M=M, N=N)
+    return h
+
+
 def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
     """
     Integrand of the finite line source solution.
@@ -596,6 +1147,117 @@ def _finite_line_source_integrand(dis, H1, D1, H2, D2, reaSource, imgSource):
     return f
 
 
+def _finite_line_source_inclined_integrand(
+        rb1, x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2, orientation2,
+        reaSource, imgSource, M):
+    """
+    Integrand of the inclined Finite Line Source (FLS) solution.
+
+    Parameters
+    ----------
+    rb1 : array
+        Radii of the emitting heat sources.
+    x1 : float or array
+        x-Positions of the emitting heat sources.
+    y1 : float or array
+        y-Positions of the emitting heat sources.
+    H1 : float or array
+        Lengths of the emitting heat sources.
+    D1 : float or array
+        Buried depths of the emitting heat sources.
+    tilt1 : float or array
+        Angles (in radians) from vertical of the emitting heat sources.
+    orientation1 : float or array
+        Directions (in radians) of the tilt the emitting heat sources.
+    x2 : array
+        x-Positions of the receiving heat sources.
+    y2 : array
+        y-Positions of the receiving heat sources.
+    H2 : float or array
+        Lengths of the receiving heat sources.
+    D2 : float or array
+        Buried depths of the receiving heat sources.
+    tilt2 : float or array
+        Angles (in radians) from vertical of the receiving heat sources.
+    orientation2 : float or array
+        Directions (in radians) of the tilt the receiving heat sources.
+    reaSource : bool
+        True if the real part of the FLS solution is to be included.
+        Default is True.
+    imgSource : bool
+        True if the image part of the FLS solution is to be included.
+    M : int
+        Number of points for the Gauss-Legendre quadrature rule along the
+        receiving heat sources.
+
+    Returns
+    -------
+    f : callable
+        Integrand of the finite line source solution. Can be vector-valued.
+
+    Notes
+    -----
+    All arrays (x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2,
+    orientation2) must follow numpy array broadcasting rules.
+
+    """
+    output_shape = np.broadcast_shapes(
+            *[np.shape(arg) for arg in (
+                rb1, x1, y1, H1, D1, tilt1, orientation1,
+                x2, y2, H2, D2, tilt2, orientation2)])
+    # Roots
+    x, w = roots_legendre(M)
+    u = (0.5 * x + 0.5).reshape((-1,) + (1,) * len(output_shape))
+    w = w / 2
+    # Params
+    sb1 = np.sin(tilt1)
+    sb2 = np.sin(tilt2)
+    cb1 = np.cos(tilt1)
+    cb2 = np.cos(tilt2)
+    dx = x1 - x2
+    dy = y1 - y2
+    if reaSource and imgSource:
+        # Full (real + image) FLS solution
+        dzRea = D1 - D2
+        dzImg = D1 + D2
+        rr = np.maximum(dx**2 + dy**2, rb1**2)
+        kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2
+        kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2
+        kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea
+        kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg
+        kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea
+        kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg
+        f = lambda s: \
+            ((H1 / s * (
+                np.exp(-(rr + dzRea**2) * s**2 + s**2 * (u**2 * H1**2 * (kRea_0**2 - 1) + 2 * u * H1 * (kRea_0 * kRea_2 - kRea_1) + kRea_2**2)) \
+                    * (erf((u *  H1 * kRea_0 + kRea_2) * s) - erf((u * H1 * kRea_0 + kRea_2 - H2) * s))
+                - np.exp(-(rr + dzImg**2) * s**2 + s**2 * (u**2 * H1**2 * (kImg_0**2 - 1) + 2 * u * H1 * (kImg_0 * kImg_2 - kImg_1) + kImg_2**2)) \
+                    * (erf((u *  H1 * kImg_0 + kImg_2) * s) - erf((u * H1 * kImg_0 + kImg_2 - H2) * s)))).T @ w).T
+    elif reaSource:
+        # Real FLS solution
+        dzRea = D1 - D2
+        rr = np.maximum(dx**2 + dy**2, rb1**2)
+        kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2
+        kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea
+        kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea
+        f = lambda s: \
+            ((H1 / s * np.exp(-(rr + dzRea**2) * s**2 + s**2 * (u**2 * H1**2 * (kRea_0**2 - 1) + 2 * u * H1 * (kRea_0 * kRea_2 - kRea_1) + kRea_2**2)) \
+                * (erf((u *  H1 * kRea_0 + kRea_2) * s) - erf((u * H1 * kRea_0 + kRea_2 - H2) * s))).T @ w).T
+    elif imgSource:
+        # Image FLS solution
+        dzImg = D1 + D2
+        kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2
+        kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg
+        kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg
+        f = lambda s: \
+            -((H1 / s * np.exp(-(dx**2 + dy**2 + dzImg**2) * s**2 + s**2 * (u**2 * H1**2 * (kImg_0**2 - 1) + 2 * u * H1 * (kImg_0 * kImg_2 - kImg_1) + kImg_2**2)) \
+                * (erf((u *  H1 * kImg_0 + kImg_2) * s) - erf((u * H1 * kImg_0 + kImg_2 - H2) * s))).T @ w).T
+    else:
+        # No heat source
+        f = lambda s: np.zeros(output_shape)
+    return f
+
+
 def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2, N2, reaSource, imgSource):
     """
     Integrand of the finite line source solution.
@@ -672,7 +1334,7 @@ def _finite_line_source_equivalent_boreholes_integrand(dis, wDis, H1, D1, H2, D2
 
 def _finite_line_source_steady_state(dis, H1, D1, H2, D2, reaSource, imgSource):
     """
-    Integrand of the finite line source solution.
+    Steady-state finite line source solution.
 
     Parameters
     ----------
@@ -741,3 +1403,113 @@ def _finite_line_source_steady_state(dis, H1, D1, H2, D2, reaSource, imgSource):
         h = np.zeros(np.broadcast_shapes(
             *[np.shape(arg) for arg in (dis, H1, D1, H2, D2)]))
     return h
+
+
+def _finite_line_source_inclined_steady_state(
+        rb1, x1, y1, H1, D1, tilt1, orientation1, x2, y2, H2, D2, tilt2,
+        orientation2, reaSource, imgSource, M=11):
+    """
+    Steady-state inclined Finite Line Source (FLS) solution.
+
+    Parameters
+    ----------
+    rb1 : array
+        Radii of the emitting heat sources.
+    x1 : float or array
+        x-Positions of the emitting heat sources.
+    y1 : float or array
+        y-Positions of the emitting heat sources.
+    H1 : float or array
+        Lengths of the emitting heat sources.
+    D1 : float or array
+        Buried depths of the emitting heat sources.
+    tilt1 : float or array
+        Angles (in radians) from vertical of the emitting heat sources.
+    orientation1 : float or array
+        Directions (in radians) of the tilt the emitting heat sources.
+    x2 : array
+        x-Positions of the receiving heat sources.
+    y2 : array
+        y-Positions of the receiving heat sources.
+    H2 : float or array
+        Lengths of the receiving heat sources.
+    D2 : float or array
+        Buried depths of the receiving heat sources.
+    tilt2 : float or array
+        Angles (in radians) from vertical of the receiving heat sources.
+    orientation2 : float or array
+        Directions (in radians) of the tilt the receiving heat sources.
+    reaSource : bool
+        True if the real part of the FLS solution is to be included.
+        Default is True.
+    imgSource : bool
+        True if the image part of the FLS solution is to be included.
+    M : int
+        Number of points for the Gauss-Legendre quadrature rule along the
+        receiving heat sources.
+        Default is 11.
+
+    Returns
+    -------
+    h : Steady-state inclined finite line source solution.
+
+    Notes
+    -----
+    All arrays (dis, H1, D1, H2, D2) must follow numpy array broadcasting
+    rules.
+
+    """
+    output_shape = np.broadcast_shapes(
+            *[np.shape(arg) for arg in (
+                rb1, x1, y1, H1, D1, tilt1, orientation1,
+                x2, y2, H2, D2, tilt2, orientation2)])
+    # Roots
+    x, w = roots_legendre(M)
+    u = (0.5 * x + 0.5).reshape((-1,) + (1,) * len(output_shape))
+    w = w / 2
+    # Params
+    sb1 = np.sin(tilt1)
+    sb2 = np.sin(tilt2)
+    cb1 = np.cos(tilt1)
+    cb2 = np.cos(tilt2)
+    dx = x1 - x2
+    dy = y1 - y2
+    # Steady-state solution
+    if reaSource and imgSource:
+        # Full (real + image) FLS solution
+        dzRea = D1 - D2
+        dzImg = D1 + D2
+        rr = np.maximum(dx**2 + dy**2, rb1**2)
+        kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2
+        kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2
+        kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea
+        kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg
+        kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea
+        kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg
+        h = 0.5 * H1 / H2 * ((
+            np.log((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + H2**2 + 2*u*H1*kRea_1 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) + 2*H2**2 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) / ((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + 2*u*H1*kRea_1) - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0)))
+            - np.log((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + H2**2 + 2*u*H1*kImg_1 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) + 2*H2**2 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) / ((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + 2*u*H1*kImg_1) - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0)))
+            ).T @ w).T
+    elif reaSource:
+        # Real FLS solution
+        dzRea = D1 - D2
+        rr = np.maximum(dx**2 + dy**2, rb1**2)
+        kRea_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) + cb1 * cb2
+        kRea_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzRea
+        kRea_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) + cb2 * dzRea
+        h = 0.5 * H1 / H2 * ((
+            np.log((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + H2**2 + 2*u*H1*kRea_1 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) + 2*H2**2 - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0) / ((2 * H2 * np.sqrt(rr + dzRea**2 + u**2*H1**2 + 2*u*H1*kRea_1) - 2*H2*kRea_2 - 2*u*H1*H2*kRea_0)))
+            ).T @ w).T
+    elif imgSource:
+        # Image FLS solution
+        dzImg = D1 + D2
+        kImg_0 = sb1 * sb2 * np.cos(orientation1 - orientation2) - cb1 * cb2
+        kImg_1 = sb1 * (np.cos(orientation1) * dx + np.sin(orientation1) * dy) + cb1 * dzImg
+        kImg_2 = sb2 * (np.cos(orientation2) * dx + np.sin(orientation2) * dy) - cb2 * dzImg
+        h = 0.5 * H1 / H2 * ((
+            - np.log((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + H2**2 + 2*u*H1*kImg_1 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) + 2*H2**2 - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0) / ((2 * H2 * np.sqrt(dx**2 + dy**2 + dzImg**2 + u**2*H1**2 + 2*u*H1*kImg_1) - 2*H2*kImg_2 - 2*u*H1*H2*kImg_0)))
+            ).T @ w).T
+    else:
+        # No heat source
+        h = np.zeros(output_shape)
+    return h
diff --git a/pygfunction/utilities.py b/pygfunction/utilities.py
index b6e6d922..0cf6690a 100644
--- a/pygfunction/utilities.py
+++ b/pygfunction/utilities.py
@@ -6,6 +6,22 @@
 import warnings
 
 
+def cardinal_point(direction):
+    accepted_inputs = ['E', 'ENE', 'NE', 'NNE',
+                       'N', 'NNW', 'NW', 'WNW',
+                       'W', 'WSW', 'SW', 'SSW',
+                       'S', 'SSE', 'SE', 'ESE']
+    direction = direction.upper()
+    if direction not in accepted_inputs:
+        raise ValueError("The indicated direction {} is not encompassed in the "
+                         "cardinal point function.".format(direction))
+    # Create the cardinal compass
+    compass = {}
+    for i, pointer in enumerate(accepted_inputs):
+        compass[pointer] = i * np.pi / 8.
+    return compass[direction]
+
+
 def erfint(x):
     """
     Integral of the error function.


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index ea33edef..a1d08193 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -27,6 +27,19 @@ def single_borehole_short():
     return [gt.boreholes.Borehole(H, D, r_b, x, y)]
 
 
+@pytest.fixture
+def single_borehole_inclined():
+    H = 150.        # Borehole length [m]
+    D = 4.          # Borehole buried depth [m]
+    r_b = 0.075     # Borehole radius [m]
+    x = 0.          # Borehole x-position [m]
+    y = 0.          # Borehole y-position [m]
+    tilt = np.pi/6          # Borehole inclination [rad]
+    orientation = np.pi/3   # Borehole orientation [rad]
+    return [gt.boreholes.Borehole(
+        H, D, r_b, x, y, tilt=tilt, orientation=orientation)]
+
+
 @pytest.fixture
 def ten_boreholes_rectangular():
     H = 150.            # Borehole length [m]
@@ -43,6 +56,12 @@ def three_boreholes_unequal():
             gt.boreholes.Borehole(177., 5., 0.085, -1., 7.),]
 
 
+@pytest.fixture
+def two_boreholes_inclined():
+    return [gt.boreholes.Borehole(150., 4., 0.075, 0., 0., tilt=np.radians(20.), orientation=np.pi/2),
+            gt.boreholes.Borehole(150., 4., 0.075, 15., 0., tilt=np.radians(20.), orientation=3*np.pi/2),]
+
+
 # =============================================================================
 # gfunction fixtures
 # =============================================================================
@@ -50,6 +69,7 @@ def three_boreholes_unequal():
 def uniform_segments():
     options = {'nSegments': 12,
                'segment_ratios': None,
+               'mQuad': 21,
                'approximate_FLS': False}
     return options
 
@@ -58,6 +78,7 @@ def uniform_segments():
 def uniform_segments_approx():
     options = {'nSegments': 12,
                'segment_ratios': None,
+               'mQuad': 21,
                'approximate_FLS': True,
                'nFLS': 10}
     return options
@@ -68,6 +89,7 @@ def unequal_segments():
     segment_ratios = gt.utilities.segment_ratios(8, end_length_ratio=0.02)
     options = {'nSegments': 8,
                'segment_ratios': segment_ratios,
+               'mQuad': 21,
                'approximate_FLS': False}
     return options
 
@@ -77,6 +99,7 @@ def unequal_segments_approx():
     segment_ratios = gt.utilities.segment_ratios(8, end_length_ratio=0.02)
     options = {'nSegments': 8,
                'segment_ratios': segment_ratios,
+               'mQuad': 21,
                'approximate_FLS': True,
                'nFLS': 10}
     return options
diff --git a/tests/gfunction_test.py b/tests/gfunction_test.py
index a525b569..2cfab47f 100644
--- a/tests/gfunction_test.py
+++ b/tests/gfunction_test.py
@@ -7,6 +7,9 @@
 import pygfunction as gt
 
 
+# =============================================================================
+# Test gFunction (vertical boreholes)
+# =============================================================================
 # Test 'UBWT' g-functions for different bore fields using all solvers,
 # unequal/uniform segments, and with/without the FLS approximation
 @pytest.mark.parametrize("field, method, opts, expected", [
@@ -227,3 +230,45 @@ def test_gfunctions_MIFT(field, method, opts, expected, request):
         network, alpha, time=time, method=method, options=options,
         boundary_condition='MIFT')
     assert np.allclose(gFunc.gFunc, expected)
+
+
+# =============================================================================
+# Test gFunction (inclined boreholes)
+# =============================================================================
+# Test 'UBWT' g-functions for a field of inclined boreholes using all solvers,
+# unequal/uniform segments, and with/without the FLS approximation
+@pytest.mark.parametrize("method, opts, expected", [
+        # 'similarities' solver - unequal segments
+        ('similarities', 'unequal_segments', np.array([5.67249989, 6.72866814, 7.15134705])),
+        # 'similarities' solver - uniform segments
+        ('similarities', 'uniform_segments', np.array([5.68324619, 6.74356205, 7.16738741])),
+        # 'similarities' solver - unequal segments, FLS approximation
+        ('similarities', 'unequal_segments_approx', np.array([5.66984803, 6.72564218, 7.14826009])),
+        # 'similarities' solver - uniform segments, FLS approximation
+        ('similarities', 'uniform_segments_approx', np.array([5.67916493, 6.7395222 , 7.16339216])),
+        # 'detailed' solver - unequal segments
+        ('detailed', 'unequal_segments', np.array([5.67249989, 6.72866814, 7.15134705])),
+        # 'detailed' solver - uniform segments
+        ('detailed', 'uniform_segments', np.array([5.68324619, 6.74356205, 7.16738741])),
+        # 'detailed' solver - unequal segments, FLS approximation
+        ('detailed', 'unequal_segments_approx', np.array([5.66984803, 6.72564218, 7.14826009])),
+        # 'detailed' solver - uniform segments, FLS approximation
+        ('detailed', 'uniform_segments_approx', np.array([5.67916493, 6.7395222 , 7.16339216])),
+    ])
+def test_gfunctions_UBWT(two_boreholes_inclined, method, opts, expected, request):
+    # Extract the bore field from the fixture
+    boreholes = two_boreholes_inclined
+    # Extract the g-function options from the fixture
+    options = request.getfixturevalue(opts)
+    # Mean borehole length [m]
+    H_mean = np.mean([b.H for b in boreholes])
+    alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+    # Bore field characteristic time [s]
+    ts = H_mean**2 / (9 * alpha)
+    # Times for the g-function [s]
+    time = np.array([0.1, 1., 10.]) * ts
+    # g-Function
+    gFunc = gt.gfunction.gFunction(
+        boreholes, alpha, time=time, method=method, options=options,
+        boundary_condition='UBWT')
+    assert np.allclose(gFunc.gFunc, expected)
diff --git a/tests/heat_transfer_test.py b/tests/heat_transfer_test.py
index 6fa36b2b..bc0a263a 100644
--- a/tests/heat_transfer_test.py
+++ b/tests/heat_transfer_test.py
@@ -9,7 +9,7 @@
 
 
 # =============================================================================
-# Test finite_line_source
+# Test finite_line_source (vertical boreholes)
 # =============================================================================
 # Test finite_line_source for single borehole to single borehole solution over
 # a single time step, using different combinations of reaSource, imgSource, and
@@ -36,7 +36,7 @@
         ('single_borehole', 'single_borehole_short', False, True, True, 10, -0.36968273),
         ('single_borehole', 'single_borehole_short', False, False, True, 10, 0.),
     ])
-def test_finite_line_source_one_to_one_single_time_step(
+def test_finite_line_source_vertical_to_vertical_single_time_step(
         borehole1, borehole2, reaSource, imgSource, approximation, N, expected,
         request):
         # Extract boreholes from fixtures
@@ -54,7 +54,7 @@ def test_finite_line_source_one_to_one_single_time_step(
 
 
 # Test finite_line_source for single borehole to single borehole solution over
-# a multiple time steps, using different combinations of reaSource, imgSource,
+# multiple time steps, using different combinations of reaSource, imgSource,
 # and with/without the FLS approximation
 @pytest.mark.parametrize("borehole1, borehole2, reaSource, imgSource, approximation, N, expected", [
         # Same borehole
@@ -78,7 +78,7 @@ def test_finite_line_source_one_to_one_single_time_step(
         ('single_borehole', 'single_borehole_short', False, True, True, 10, np.array([-0.07114505, -0.36968273, -0.74142343])),
         ('single_borehole', 'single_borehole_short', False, False, True, 10, np.array([0., 0., 0.])),
     ])
-def test_finite_line_source_one_to_one_multiple_time_steps(
+def test_finite_line_source_vertical_to_vertical_multiple_time_steps(
         borehole1, borehole2, reaSource, imgSource, approximation, N, expected,
         request):
         # Extract boreholes from fixtures
@@ -109,7 +109,7 @@ def test_finite_line_source_one_to_one_multiple_time_steps(
         ('single_borehole', 'single_borehole_short', False, True, -0.9988295508738515),
         ('single_borehole', 'single_borehole_short', False, False, 0.),
     ])
-def test_finite_line_source_one_to_one_steady_state(
+def test_finite_line_source_vertical_to_vertical_steady_state(
         borehole1, borehole2, reaSource, imgSource, expected, request):
         b1 = request.getfixturevalue(borehole1)[0]
         b2 = request.getfixturevalue(borehole2)[0]
@@ -155,7 +155,7 @@ def test_finite_line_source_one_to_one_steady_state(
              [-0.08137532, -0.08865631, -0.08193464]])),
         (False, False, True, 10, np.zeros((3, 3))),
     ])
-def test_finite_line_source_multiple_boreholes_single_time_step(
+def test_finite_line_source_multiple_vertical_boreholes_single_time_step(
         three_boreholes_unequal, reaSource, imgSource, approximation, N,
         expected):
         # Extract boreholes from fixture
@@ -241,7 +241,7 @@ def test_finite_line_source_multiple_boreholes_single_time_step(
               [-0.01030019, -0.08193464, -0.30431902]]])),
         (False, False, True, 10, np.zeros((3, 3, 3))),
     ])
-def test_finite_line_source_multiple_boreholes_multiple_time_steps(
+def test_finite_line_source_multiple_vertical_boreholes_multiple_time_steps(
         three_boreholes_unequal, reaSource, imgSource, approximation, N,
         expected):
         # Extract boreholes from fixture
@@ -275,7 +275,7 @@ def test_finite_line_source_multiple_boreholes_multiple_time_steps(
              [-0.54949639, -0.40882024, -0.60213297]])),
         (False, False, np.zeros((3, 3))),
     ])
-def test_finite_line_source_multiple_boreholes_steady_state(
+def test_finite_line_source_multiple_vertical_boreholes_steady_state(
         three_boreholes_unequal, reaSource, imgSource, expected):
         # Extract boreholes from fixture
         boreholes = three_boreholes_unequal
@@ -287,3 +287,241 @@ def test_finite_line_source_multiple_boreholes_steady_state(
             time, alpha, boreholes, boreholes, reaSource=reaSource,
             imgSource=imgSource)
         assert np.allclose(h, expected)
+
+
+# =============================================================================
+# Test finite_line_source (inclined boreholes)
+# =============================================================================
+# Test finite_line_source for the test case of Lazzarotto (2016) over a single
+# and multiple time steps, using different combinations of reaSource,
+# imgSource, and with/without the FLS approximation
+@pytest.mark.parametrize("tts, reaSource, imgSource, approximation, M, N, expected", [
+        # One time step
+        (1., True, True, False, 21, 10, 0.055792623436062005),
+        (1., True, False, False, 21, 10, 0.055798508426994645),
+        (1., False, True, False, 21, 10, -5.884990932696544e-06),
+        (1., False, False, False, 21, 10, 0.),
+        # One time step - FLS approximation
+        (1., True, True, True, 21, 10, 0.055807785298494346),
+        (1., True, False, True, 21, 10, 0.05581324608019175),
+        (1., False, True, True, 21, 10, -5.4607816973996686e-06),
+        (1., False, False, True, 21, 10, 0.),
+        # Multiple time steps
+        (np.array([0.1, 1., 10., 100.]), True, True, False, 21, 10, np.array(
+            [2.03508637e-04, 5.57926234e-02, 2.65993224e-01, 3.25925673e-01])),
+        (np.array([0.1, 1., 10., 100.]), True, False, False, 21, 10, np.array(
+            [2.03508637e-04, 5.57985084e-02, 2.89690224e-01, 4.50763618e-01])),
+        (np.array([0.1, 1., 10., 100.]), False, True, False, 21, 10, np.array(
+            [0., -5.88499093e-06, -2.36969998e-02, -1.24837945e-01])),
+        (np.array([0.1, 1., 10., 100.]), False, False, False, 21, 10, np.array(
+            [0., 0., 0., 0.])),
+        # Multiple time steps - FLS approximation
+        (np.array([0.1, 1., 10., 100.]), True, True, True, 21, 10, np.array(
+            [2.03364838e-04, 5.58077853e-02, 2.65949066e-01, 3.26037282e-01])),
+        (np.array([0.1, 1., 10., 100.]), True, False, True, 21, 10, np.array(
+            [2.03364838e-04, 5.58132461e-02, 2.89638257e-01, 4.50797447e-01])),
+        (np.array([0.1, 1., 10., 100.]), False, True, True, 21, 10, np.array(
+            [0., -5.46078170e-06, -2.36891904e-02, -1.24760165e-01])),
+        (np.array([0.1, 1., 10., 100.]), False, False, True, 21, 10, np.array(
+            [0., 0., 0., 0.])),
+    ])
+def test_finite_line_source_Lazzarotto(
+        tts, reaSource, imgSource, approximation, M, N, expected):
+        # Extract boreholes from fixtures
+        b1 = gt.boreholes.Borehole(
+            20., 5., 0.075, 3., 1.5, tilt=np.pi/15, orientation=0.)
+        b2 = gt.boreholes.Borehole(
+            20., 25., 0.075, 0., 5., tilt=np.pi/15, orientation=4*np.pi/3)
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        ts = b1.H**2 / (9 * alpha)  # Borehole characteristic time [s]
+        # Time for FLS calculation [s]
+        time = ts * tts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b2, reaSource=reaSource, imgSource=imgSource,
+            approximation=approximation, M=M, N=N)
+        assert np.allclose(h, expected)
+
+
+# Test finite_line_source for the test case of Lazzarotto (2016) at
+# steady-state, using different combinations of reaSource and imgSource
+@pytest.mark.parametrize("reaSource, imgSource, M, expected", [
+        (True, True, 21, 0.32901256907275805),
+        (True, False, 21, 0.5345970978237233),
+        (False, True, 21, -0.2055845287509652),
+        (False, False, 21, 0.),
+    ])
+def test_finite_line_source_Lazzarotto_steady_state(
+        reaSource, imgSource, M, expected):
+        # Configure boreholes
+        b1 = gt.boreholes.Borehole(
+            20., 5., 0.075, 3., 1.5, tilt=np.pi/15, orientation=0.)
+        b2 = gt.boreholes.Borehole(
+            20., 25., 0.075, 0., 5., tilt=np.pi/15, orientation=4*np.pi/3)
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        # Time for FLS calculation [s]
+        time = np.inf
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b2, reaSource=reaSource, imgSource=imgSource, M=M)
+        assert np.isclose(h, expected)
+
+
+# Test finite_line_source for an inclined borehole on itself for a single
+# and multiple time steps, using different combinations of reaSource,
+# imgSource, and with/without the FLS approximation
+@pytest.mark.parametrize("tts, reaSource, imgSource, approximation, M, N, expected", [
+        # One time step
+        (1., True, True, False, 21, 10, 6.387924910257989),
+        (1., True, False, False, 21, 10, 6.5328321472166095),
+        (1., False, True, False, 21, 10, -0.14490784205005186),
+        (1., False, False, False, 21, 10, 0.),
+        # One time step - FLS approximation
+        (1., True, True, True, 21, 10, 6.38795301165354),
+        (1., True, False, True, 21, 10, 6.532855370434292),
+        (1., False, True, True, 21, 10, -0.1449023587807526),
+        (1., False, False, True, 21, 10, 0.),
+        # Multiple time steps
+        (np.array([0.1, 1., 10.]), True, True, False, 21, 10, np.array(
+            [5.61499596, 6.38792491, 6.61329376])),
+        (np.array([0.1, 1., 10.]), True, False, False, 21, 10, np.array(
+            [5.63802126, 6.53283215, 7.03111215])),
+        (np.array([0.1, 1., 10.]), False, True, False, 21, 10, np.array(
+            [-0.0230257 , -0.14490784, -0.41781906])),
+        (np.array([0.1, 1., 10.]), False, False, False, 21, 10, np.array(
+            [0., 0., 0.])),
+        # Multiple time steps - FLS approximation
+        (np.array([0.1, 1., 10.]), True, True, True, 21, 10, np.array(
+            [5.61500526, 6.38795301, 6.61341841])),
+        (np.array([0.1, 1., 10.]), True, False, True, 21, 10, np.array(
+            [5.63803137, 6.53285537, 7.03118607])),
+        (np.array([0.1, 1., 10.]), False, True, True, 21, 10, np.array(
+            [-0.02302611, -0.14490236, -0.41776766])),
+        (np.array([0.1, 1., 10.]), False, False, True, 21, 10, np.array(
+            [0., 0., 0.])),
+    ])
+def test_finite_line_source_inclined_to_self(
+        single_borehole_inclined, tts, reaSource, imgSource, approximation,
+        M, N, expected):
+        # Extract borehole from fixtures
+        b1 = single_borehole_inclined[0]
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        ts = b1.H**2 / (9 * alpha)  # Borehole characteristic time [s]
+        # Time for FLS calculation [s]
+        time = ts * tts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b1, reaSource=reaSource, imgSource=imgSource,
+            approximation=approximation, M=M, N=N)
+        assert np.allclose(h, expected)
+
+
+# Test finite_line_source for an inclined borehole on itself at
+# steady-state, using different combinations of reaSource and imgSource
+@pytest.mark.parametrize("reaSource, imgSource, M, expected", [
+        (True, True, 21, 6.628532759512739),
+        (True, False, 21, 7.295473329568198),
+        (False, True, 21, -0.6669405700554597),
+        (False, False, 21, 0.),
+    ])
+def test_finite_line_source_inclined_to_self_steady_state(
+        single_borehole_inclined, reaSource, imgSource, M, expected):
+        # Extract borehole from fixtures
+        b1 = single_borehole_inclined[0]
+        alpha = 1e-6                # Ground thermal diffusivity [m2/s]
+        # Time for FLS calculation [s]
+        time = np.inf
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, b1, b1, reaSource=reaSource, imgSource=imgSource, M=M)
+        assert np.isclose(h, expected)
+
+
+# Test finite_line_source for multiple inclined boreholes to multiple inclined
+# boreholes solution over a single or multiple time steps, using different
+# combinations of reaSource, imgSource, and with/without the FLS approximation
+@pytest.mark.parametrize("tts, reaSource, imgSource, approximation, M, N, expected", [
+        # One time step
+        (1., True, True, False, 21, 10, np.array(
+            [[6.40344545, 0.38333631],
+             [0.38333631, 6.40344545]])),
+        (1., True, False, False, 21, 10, np.array(
+            [[6.53283215, 0.48927282],
+             [0.48927282, 6.53283215]])),
+        (1., False, True, False, 21, 10, np.array(
+            [[-0.1293867 , -0.10593651],
+             [-0.10593651, -0.1293867]])),
+        (1., False, False, False, 21, 10, np.zeros((2, 2))),
+        # One time step - FLS approximation
+        (1., True, True, True, 21, 10, np.array(
+            [[6.40347323, 0.3833017],
+             [0.3833017 , 6.40347323]])),
+        (1., True, False, True, 21, 10, np.array(
+            [[6.53285537, 0.48923122],
+             [0.48923122, 6.53285537]])),
+        (1., False, True, True, 21, 10, np.array(
+            [[-0.12938214, -0.10592953],
+             [-0.10592953, -0.12938214]])),
+        (1., False, False, True, 21, 10, np.zeros((2, 2))),
+        # Multiple time steps
+        (np.array([0.1, 1., 10.]), True, True, False, 21, 10, np.array(
+            [[[5.61764541, 0.0739292],
+              [0.0739292, 5.61764541]],
+             [[6.40344545, 0.38333631],
+              [0.38333631, 6.40344545]],
+             [[6.64617823, 0.5812878],
+              [0.5812878, 6.64617823]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), True, False, False, 21, 10, np.array(
+            [[[5.63802126, 0.08504218],
+              [0.08504218, 5.63802126]],
+             [[6.53283215, 0.48927282],
+              [0.48927282, 6.53283215]],
+             [[7.03111215, 0.9284482],
+              [0.9284482, 7.03111215]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), False, True, False, 21, 10, np.array(
+            [[[-0.02037585, -0.01111298],
+              [-0.01111298, -0.02037585]],
+             [[-0.1293867, -0.10593651],
+              [-0.10593651, -0.1293867]],
+             [[-0.38493416, -0.34716039],
+              [-0.34716039, -0.38493416]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), False, False, False, 21, 10, np.zeros((2, 2, 3))),
+        # Multiple time steps - FLS approximation
+        (np.array([0.1, 1., 10.]), True, True, True, 21, 10, np.array(
+            [[[5.61765526, 0.07394688],
+              [0.07394688, 5.61765526]],
+             [[6.40347323, 0.3833017],
+              [0.3833017, 6.40347323]],
+             [[6.64629192, 0.58122997],
+              [0.58122997, 6.64629192]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), True, False, True, 21, 10, np.array(
+            [[[5.63803137, 0.08506054],
+              [0.08506054, 5.63803137]],
+             [[6.53285537, 0.48923122],
+              [0.48923122, 6.53285537]],
+             [[7.03118607, 0.9283558],
+              [0.9283558, 7.03118607]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), False, True, True, 21, 10, np.array(
+            [[[-0.02037611, -0.01111366],
+              [-0.01111366, -0.02037611]],
+             [[-0.12938214, -0.10592953],
+              [-0.10592953, -0.12938214]],
+             [[-0.38489415, -0.34712584],
+              [-0.34712584, -0.38489415]]]).transpose((1, 2, 0))),
+        (np.array([0.1, 1., 10.]), False, False, True, 21, 10, np.zeros((2, 2, 3))),
+    ])
+def test_finite_line_source_multiple_inclined_to_multiple_inclined(
+        two_boreholes_inclined, tts, reaSource, imgSource, approximation, M, N,
+        expected):
+        # Extract boreholes from fixture
+        boreholes = two_boreholes_inclined
+        alpha = 1e-6    # Ground thermal diffusivity [m2/s]
+        # Bore field characteristic time [s]
+        ts = np.mean([b.H for b in boreholes])**2 / (9 * alpha)
+        # Times for FLS calculation [s]
+        time = ts * tts
+        # Evaluate FLS
+        h = gt.heat_transfer.finite_line_source(
+            time, alpha, boreholes, boreholes, reaSource=reaSource,
+            imgSource=imgSource, approximation=approximation, M=M, N=N)
+        assert np.allclose(h, expected)