diff --git a/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py b/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py
index 0c755d7..e479169 100644
--- a/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py
+++ b/openfast_toolbox/fastfarm/FASTFarmCaseCreation.py
@@ -531,6 +531,7 @@ def _checkInputs(self):
         self.hasSrvD                   = False
         self.hasHD                     = False
         self.hasMD                     = False
+        self.hasBath                   = False
         self.hasSS                     = False
         self.hasSubD                   = False
         self.hasBD                     = False
@@ -824,6 +825,9 @@ def copyTurbineFilesForEachCase(self, writeFiles=True):
                             self.MoorDynFile.write(os.path.join(currPath,f'{self.MDfilename}{t+1}_mod.dat'))
                         else:
                             if t==0: shutilcopy2_untilSuccessful(self.MDfilepath, os.path.join(currPath, self.MDfilename))
+                    
+                    if writeFiles and self.hasBath:
+                        if t==0: shutilcopy2_untilSuccessful(self.bathfilepath, os.path.join(currPath, self.bathfilename))
 
                     # Update each turbine's OpenFAST input
                     self.turbineFile['TMax']         = self.tmax
@@ -1043,6 +1047,7 @@ def setTemplateFilename(self, templatePath=None, templateFiles=None):
             'SEDfilename'             : None,
             'HDfilename'              : 'HydroDyn.dat', # ending with .T for per-turbine HD, .dat for holisitc
             'MDfilename'              : 'MoorDyn.T',    # ending with .T for per-turbine MD, .dat for holistic
+            'bathfilename'            : 'bathymetry.txt',   
             'SSfilename'              : 'SeaState.dat',
             'SrvDfilename'            : 'ServoDyn.T',
             'ADfilename'              : 'AeroDyn.dat',
@@ -1072,6 +1077,7 @@ def setTemplateFilename(self, templatePath=None, templateFiles=None):
         self.SEDfilename     = "unused";  self.SEDfilepath     = "unused"
         self.HDfilename      = "unused";  self.HDfilepath      = "unused"
         self.MDfilename      = "unused";  self.MDfilepath      = "unused"
+        self.bathfilename    = "unused";  self.bathfilepath    = "unused"
         self.SSfilename      = "unused";  self.SSfilepath      = "unused"
         self.SrvDfilename    = "unused";  self.SrvDfilepath    = "unused"
         self.ADfilename      = "unused";  self.ADfilepath      = "unused"
@@ -1102,7 +1108,7 @@ def setTemplateFilename(self, templatePath=None, templateFiles=None):
         self.templatePathabs = os.path.abspath(self.templatePath)
 
         # Check and set the templateFiles
-        valid_keys = {'EDfilename', 'SEDfilename', 'HDfilename', 'MDfilename', 'SSfilename',
+        valid_keys = {'EDfilename', 'SEDfilename', 'HDfilename', 'MDfilename', 'bathfilename', 'SSfilename',
                       'SrvDfilename', 'ADfilename', 'ADskfilename', 'SubDfilename', 'IWfilename', 'BDfilename',
                       'BDbladefilename', 'EDbladefilename', 'EDtowerfilename', 'ADbladefilename', 'turbfilename',
                       'libdisconfilepath', 'controllerInputfilename', 'coeffTablefilename', 'hydroDatapath',
@@ -1168,6 +1174,11 @@ def checkIfExists(f):
                 self.MDfilename = value
                 self.hasMD = True
 
+            elif key == 'bathfilename':
+                self.bathfilepath = os.path.join(self.templatePath, value)
+                checkIfExists(self.bathfilepath)
+                self.bathfilename = value
+                self.hasBath = True
             elif key == 'SSfilename':
                 if not value.endswith('.dat'):
                     raise ValueError(f'The SeaState filename should end in `.dat`.')
@@ -1811,7 +1822,7 @@ def getDomainParameters(self):
         
         # Determine offsets from turbines coordinate frame to TurbSim coordinate frame
         self.yoffset_turbsOrigin2TSOrigin = -( (self.TSlowbox.ymax - self.TSlowbox.ymin)/2 + self.TSlowbox.ymin )
-        self.xoffset_turbsOrigin2TSOrigin = - self.extent_low[0]*self.D
+        self.xoffset_turbsOrigin2TSOrigin = -self.extent_low[0]*self.D
         
         if self.verbose>0:
             print(f"    The y offset between the turbine ref frame and turbsim is {self.yoffset_turbsOrigin2TSOrigin}")
@@ -2138,8 +2149,8 @@ def FF_setup(self, outlistFF=None, **kwargs):
         
         offset=10
         planes_xy = [self.zhub+self.zbot]
-        planes_yz = np.unique(xWT+offset)
-        planes_xz = np.unique(yWT)
+        planes_yz = np.unique(np.round(xWT+offset, 2))
+        planes_xz = np.unique(np.round(yWT, 2))
         
         # Number of planes must be at most 9
         self.planes_xy = planes_xy[0:9]
@@ -2750,11 +2761,12 @@ def plot(self, figsize=(14,7), fontsize=13, saveFig=True, returnFig=False, figFo
 
                 # plot turbine disk accoding to all yaws in current wdir
                 allyaw_currwdir = self.allCases.where(self.allCases['inflow_deg']==inflow,drop=True).sel(turbine=currTurbine)['yaw']
+                phi             = self.allCases.where(self.allCases['inflow_deg']==inflow,drop=True).sel(turbine=currTurbine)['phi']
                 _, ind = np.unique(allyaw_currwdir, axis=0, return_index=True)
                 yaw_currwdir = allyaw_currwdir[np.sort(ind)].values # duplicates removed, same order as original array
                 for yaw in yaw_currwdir:
-                    ax.plot([dst.x.values-(dst.D.values/2)*sind(yaw), dst.x.values+(dst.D.values/2)*sind(yaw)],
-                            [dst.y.values-(dst.D.values/2)*cosd(yaw), dst.y.values+(dst.D.values/2)*cosd(yaw)], c=color, alpha=alphas[j])
+                    ax.plot([dst.x.values-(dst.D.values/2)*sind(yaw+phi), dst.x.values+(dst.D.values/2)*sind(yaw+phi)],
+                            [dst.y.values-(dst.D.values/2)*cosd(yaw+phi), dst.y.values+(dst.D.values/2)*cosd(yaw+phi)], c=color, alpha=alphas[j])
 
             # plot convex hull of farm (or line) for given inflow
             turbs = self.wts_rot_ds.sel(inflow_deg=inflow)[['x','y']].to_array().transpose()
@@ -2772,7 +2784,7 @@ def plot(self, figsize=(14,7), fontsize=13, saveFig=True, returnFig=False, figFo
         handles, labels = plt.gca().get_legend_handles_labels()
         by_label = dict(zip(labels, handles))
         if showLegend:
-            plt.legend(by_label.values(), by_label.keys(), loc='upper left', bbox_to_anchor=(1.02,1.015), fontsize=fontsize, ncols=int(self.nTurbines/25))
+            plt.legend(by_label.values(), by_label.keys(), loc='upper left', bbox_to_anchor=(1.02,1.015), fontsize=fontsize, ncols=np.ceil(self.nTurbines/25))
 
         ax.set_xlabel("x [m]", fontsize=fontsize)
         ax.set_ylabel("y [m]", fontsize=fontsize)