diff --git a/src/Model/Geometry/Reach.py b/src/Model/Geometry/Reach.py
index ebff2dc7..b2e35f7b 100644
--- a/src/Model/Geometry/Reach.py
+++ b/src/Model/Geometry/Reach.py
@@ -455,3 +455,121 @@ class Reach:
 
                 file_st.write("    999.9990     999.9990     999.9990")
                 file_st.write("\n")
+
+    # Interpolate
+
+    @timer
+    def interpolate_step(self, lim: list[int], guideline1: str, guideline2: str, step: float, longmi: int, linear: bool)
+        """interpolate with a given step between profile lim[0] and profile lim[1]
+        Args:
+            lim: index of profiles for interpolation
+            guideline1: first guideline for distance
+            guideline2: second guideline for distance
+            step: interpolation step
+            longmi: 1 if distance computed on min(dist(guideline1),dist(guideline2))
+                    2 if distance computed on max(dist(guideline1),dist(guideline2))
+                    else: distance computed on (dist(guideline1)+dist(guideline2))/2
+            linear: true if linear interpolation, false if spline interpolation
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        if len(lim) != 2:
+            print('TODO: CRITICAL ERROR!') # error for now, only one interval
+        if lim[0] < 0 .or. lim[0] >= lim[1] .or. lim[1] >= len(self._profiles):
+            print('TODO: CRITICAL ERROR!')
+        if lim[1] < 0 .or. lim[0] > lim[1]:
+            print('TODO: CRITICAL ERROR!')
+        comp, inconmp = self.compute_guidelines()
+        if not guideline1 in comp:
+            print('TODO: CRITICAL ERROR!')
+        if not guideline2 in comp:
+            print('TODO: CRITICAL ERROR!')
+
+        nbsti = [] # number of profiles to add between self._profile[i] and self._profile[i+1]
+        for i in range(lim[0])
+            nbsti.append(0)
+        for i in range(lim[0],lim[1]-1):
+                points = self.profile(i).points
+                point11 = next(filter(lambda p: p.name == guideline1, points))
+                point12 = next(filter(lambda p: p.name == guideline2, points))
+                points = self.profile(i+1).points
+                point21 = next(filter(lambda p: p.name == guideline1, points))
+                point22 = next(filter(lambda p: p.name == guideline2, points))
+                d1 = point11.dist(point21)
+                d2 = point12.dist(point22) # warning: 3D distance
+
+            if (longmi == 1) then
+                d = max(d1, d2)
+            else if (longmi == 2) then
+                d = min(d1, d2)
+            else
+                d = 0.5*(d1 + d2)
+
+            nbsti[i] = nint(d / pas) - 1
+            if(nbsti[i].lt.0) nbsti[i] = 0
+
+        for i in range(lim[1]-1,len(self._profiles)-1)
+            nbsti.append(0)
+
+        return self.interpolate_nprof(nbsti, lineaire)
+
+    @timer
+    def interpolate_nprof(self, nbsti: List[int], linear: bool):
+        """add nbsti[i] profiles between self._profile[i] and self._profile[i+1]
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+            linear: true if linear interpolation, false if spline interpolation
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+
+        if(lineaire):
+            return self.interpolate_linear(nbsti)
+        else:
+            # we might add the spline interpolation later
+            print('TODO: CRITICAL ERROR!')
+
+    @timer
+    def interpolate_linear(self, nbsti: List[int]):
+        """add nbsti[i] profiles between self._profile[i] and self._profile[i+1]
+
+        Args:
+            nbsti: number of profiles to add between self._profile[i] and self._profile[i+1]
+
+        Returns:
+            List of profiles including old profiles and interpolated profiles.
+        """
+        if len(nbsti) != self.number_profiles()-1:
+            print('TODO: CRITICAL ERROR!') # error for now, only one interval
+        sections_out = []
+        kps = self.get_kp()
+        for isect_in in range(len(nbsti)):
+            sections_out.append(self.section(isect_in))
+            dkp = kps[isect_in + 1] - kps[isect_in]
+            sect_in1 = self.section(isect_in)
+            sect_in2 = self.section(isect_in+1)
+            for isect_loc in range(nbsti[isect_in]):
+                # RATIO entre les deux sections initiales
+                d=float(isect_loc)/float(nbsti(isect_in)+1)
+                x1 = sect_in1.x()
+                y1 = sect_in1.y()
+                z1 = sect_in1.z()
+                x2 = sect_in2.x()
+                y2 = sect_in2.y()
+                z2 = sect_in2.z()
+                tags = sect_in1.names()
+                section_out.append(ProfileXYZ(kp = kps[isect_in] + d*dkp,
+                                              reach = self))
+                for i in range(len(x1)): # can do better...
+                    # add the points to the created section
+                    sections_out[-1].add(PointXYZ(x = x1[i] + d*(x2[i] - x1[i]),
+                                                  y = y1[i] + d*(y2[i] - y1[i]),
+                                                  z = z1[i] + d*(z2[i] - z1[i]),
+                                                  name = tags[i]
+                                                 )
+        # don't forget the last one:
+        sections_out.append(self.section(self.number_profiles()-1))
+        return sections_out