diff --git a/README.md b/README.md
index 15f18ef6347b42419b551268a48d3f7ef1882be8..0fb8d7c5b287b6b6093cafd8dd9b8762fb053d56 100644
--- a/README.md
+++ b/README.md
@@ -10,4 +10,6 @@ The so far implemented features are:
 *  calculates beam sizes and position for 2 energies
 *  calculates sample position (membrane size and pitch) as well as etched facets
 *  calculates position of the DSSC module 15 as well as filter bar position
-*  calculates focal length change due to the KBS  
\ No newline at end of file
+*  calculates focal length change due to the KBS
+*  calculates beam focusing by KBS for the zone plate 0th order
+*  calculates beam and sample shape for sample at non-normal incidence angle
diff --git a/TZPGcalc.py b/TZPGcalc.py
index 4e5129ba606d7b2d3e3ae0774172cffe2049aa3b..e934626a323deedacc5457224b61a1bfb138b627 100644
--- a/TZPGcalc.py
+++ b/TZPGcalc.py
@@ -9,7 +9,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from matplotlib.patches import Rectangle
+from matplotlib.patches import Rectangle, Polygon
 from matplotlib.colors import hsv_to_rgb
 
 import ipywidgets as widgets
@@ -26,6 +26,7 @@ Z0 = 230*1e-3 # [m]
 d = 3.3 - Z0 # distance between HFM and TZPG
 f1 = 7.3 # HFM focus 2 m behind second interaction point
 F = F*(d-f1)/(d-f1-F)
+KBS_F = f1 - d # KBS focus distance from TZPG
 
 # number of membrane to show
 SampleN = 7
@@ -63,50 +64,86 @@ class TZPGcalc():
         c_gr = hsv_to_rgb([95/360, 60/100, 100/100])
         c_gb = hsv_to_rgb([145/360, 60/100, 100/100])
 
-        self.samBeamsL = {'F0G0': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G-1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F1G0': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_rr, alpha=0.7, lw=None)),
-                'F1G1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gr, alpha=0.7, lw=None)),
-                'F1G-1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gr, alpha=0.7, lw=None))
-               }
-
-        self.detBeamsL = {'F0G0': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G-1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F1G0': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_rr, alpha=0.7, lw=None)),
-                'F1G1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gr, alpha=0.7, lw=None)),
-                'F1G-1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gr, alpha=0.7, lw=None))
-               }
-
-        self.samBeamsH = {'F0G0': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G-1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F1G0': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_rb, alpha=0.7, lw=None)),
-                'F1G1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gb, alpha=0.7, lw=None)),
-                'F1G-1': self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gb, alpha=0.7, lw=None))
-               }
-
-        self.detBeamsH = {'F0G0': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F0G-1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="black", alpha=0.4, lw=None)),
-                'F1G0': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_rb, alpha=0.7, lw=None)),
-                'F1G1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gb, alpha=0.7, lw=None)),
-                'F1G-1': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor=c_gb, alpha=0.7, lw=None))
-               }
-
-        self.detLines = {'module': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, fill=False, facecolor='k')),
-                'Vfilter': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="blue", alpha=0.4)),
-                'Hfilter': self.ax_det.add_patch(Rectangle((0, 0), 1, 1, facecolor="blue", alpha=0.4)),
-                'diamond': self.ax_det.add_patch(Rectangle((-8, -8), 16, 16, facecolor="blue", alpha=0.4, angle=45))
+        self.samBeamsL = {
+            'F0G0': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G-1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F1G0': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_rr, alpha=0.7, lw=None)),
+            'F1G1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_gr, alpha=0.7, lw=None)),
+            'F1G-1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_gr, alpha=0.7, lw=None))
+            }
+
+        self.detBeamsL = {
+            'F0G0': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G-1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F1G0': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_rr, alpha=0.7, lw=None)),
+            'F1G1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_gr, alpha=0.7, lw=None)),
+            'F1G-1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_gr, alpha=0.7, lw=None))
+            }
+
+        self.samBeamsH = {
+            'F0G0': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G-1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F1G0': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_rb, alpha=0.7, lw=None)),
+            'F1G1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_gb, alpha=0.7, lw=None)),
+            'F1G-1': self.ax_sam.add_patch(
+                Polygon([(0, 0)], facecolor=c_gb, alpha=0.7, lw=None))
+            }
+
+        self.detBeamsH = {
+            'F0G0': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F0G-1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor="black", alpha=0.4, lw=None)),
+            'F1G0': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_rb, alpha=0.7, lw=None)),
+            'F1G1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_gb, alpha=0.7, lw=None)),
+            'F1G-1': self.ax_det.add_patch(
+                Polygon([(0, 0)], facecolor=c_gb, alpha=0.7, lw=None))
+            }
+
+        self.detLines = {
+            'module': self.ax_det.add_patch(
+                Rectangle((0, 0), 1, 1, fill=False, facecolor='k')),
+            'Vfilter': self.ax_det.add_patch(
+                Rectangle((0, 0), 1, 1, facecolor="blue", alpha=0.4)),
+            'Hfilter': self.ax_det.add_patch(
+                Rectangle((0, 0), 1, 1, facecolor="blue", alpha=0.4)),
+            'diamond': self.ax_det.add_patch(
+                Rectangle((-8, -8), 16, 16, facecolor="blue", alpha=0.4, angle=45))
                }
 
         # 5x5 membranes
         self.sampleLines = {}
         self.etchLines = {}
         for k in range(SampleN*SampleN):
-            self.sampleLines[k] = self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, fill=False, facecolor='k'))
-            self.etchLines[k] = self.ax_sam.add_patch(Rectangle((0, 0), 1, 1, fill=False, facecolor='k', alpha=0.4, ls='--'))
+            self.sampleLines[k] = self.ax_sam.add_patch(
+                Rectangle((0, 0), 1, 1, fill=False, facecolor='k'))
+            self.etchLines[k] = self.ax_sam.add_patch(
+                Rectangle((0, 0), 1, 1, fill=False, facecolor='k', alpha=0.4, ls='--'))
+
 
     def RectUpdate(self, rect, xLeft, yBottom, xRight, yTop):
         """ Updates the position and size of the given Rectangle.
@@ -126,13 +163,54 @@ class TZPGcalc():
         rect.set_height(self.scale*yw)
         rect.set_width(self.scale*xw)
 
-    def UpdateBeams(self, Beams, Z, conf):
+    def PolyUpdate(self, poly, xLeft, yBottom, xRight, yTop):
+        """ Updates the corner position of a Polygon.
+
+            poly: regular Polygon to update
+            xLeft: x position of the left corner
+            yBottom: y position of the bottom corner
+            xRight: x position of the right corner
+            yTop: y position of the top corner
+        """
+
+        xy = self.scale*np.array([
+            [xLeft, yBottom],
+            [xLeft, yTop],
+            [xRight, yTop],
+            [xRight, yBottom]])
+
+        poly.set_xy(xy)
+
+    def LinePlaneIntersection(self, l1, l2, p0, n):
+        """ Calculate the intersection point in space between a line (beam)
+            passing through 2 points l1 and l2 and a (sample) plane passing by
+            p0 with normal n
+
+            l1: [x,y,z] point on line
+            l2: [x,y,z] point on line
+            p0: [x,y,z] point on plane
+            n: plane normal vector
+        """
+
+        # plane parametrized as (p - p0).n = 0
+        # line parametrized as p = l1 + l12*d with d Real
+        l12 = l2 - l1
+
+        if np.dot(l12,n) == 0:
+            return [0,0,0] # line is either in the plane or outside the plane
+        else:
+            d = np.dot((p0 - l1), n)/np.dot(l12,n)
+            return l1 + l12*d
+
+    def UpdateBeams(self, Beams, Z, incidence, conf):
         """ Update the position and size of the beams.
 
-            Beams: dictionary of f'F{f}G{g}' Rectangles for f = 0 and 1 zone plate order and
-                g = +1, 0 and -1 grating order
+            Beams: dictionary of f'F{f}G{g}' Polygon for f = 0 and 1 zone
+                plate order and g = +1, 0 and -1 grating order
             Z: distance Z between the zone plate and the current imaging plane
-            conf: dictionnary of distance for calculation {'F', 'TZPGwH', 'TZPGwV', 'TZPGo', 'theta_grating'}
+            incidence: incidence angle of the imaging plane in rad
+            conf: dictionnary of distance for calculation {'F', 'TZPGwH',
+                'TZPGwV', 'TZPGo', 'theta_grating'}
         """
         F = conf['F']
         wH = conf['TZPGwH']
@@ -140,36 +218,34 @@ class TZPGcalc():
         o = conf['TZPGo']
         offaxis = wV/2 + o
 
-        # side view X
-        Fx = F*np.arctan(conf['theta_grating'])
-
-        Xdg0 = np.abs(Z - F)/F*wH/2
-        Xdg1L = Z/F*(Fx - wH/2) + wH/2
-        Xdg1 = Z*np.tan(conf['theta_grating'])
-        Xdg1H = Z/F*(Fx + wH/2) - wH/2
-
-        # top view Y
-        YdfL = (o*(Z - F)/F + offaxis)
-        Ydf = ((o + wV/2)*(Z - F)/F + offaxis)
-        YdfH = ((o + wV)*(Z - F)/F + offaxis)
-
-        if Z < F: # before zone plate focus, low and high beam edges are swapped
-            Xdg1L, Xdg1H = (Xdg1H, Xdg1L)
-            YdfL, YdfH = (YdfH, YdfL)
-
-        self.RectUpdate(Beams['F0G0'],
-                   -wH/2, -wV/2, wH/2, wV/2)
-        self.RectUpdate(Beams['F0G1'],
-                   Xdg1-wH/2, -wV/2, Xdg1+wH/2, wV/2)
-        self.RectUpdate(Beams['F0G-1'],
-                   -Xdg1-wH/2, -wV/2, -Xdg1+wH/2, wV/2)
-
-        self.RectUpdate(Beams['F1G0'],
-                   -Xdg0, YdfL, Xdg0, YdfH)
-        self.RectUpdate(Beams['F1G1'],
-                   Xdg1L, YdfL, Xdg1H, YdfH)
-        self.RectUpdate(Beams['F1G-1'],
-                   -Xdg1H, YdfL, -Xdg1L, YdfH)
+        # imaging plane
+        n = np.array([np.sin(incidence), 0, np.cos(incidence)])
+        p0 = np.array([0,0,Z])
+
+        # zone plate 4 corner points
+        l1_list = [
+            np.array([wH/2,wV/2,0]),
+            np.array([wH/2,-wV/2,0]),
+            np.array([-wH/2,-wV/2,0]),
+            np.array([-wH/2,wV/2,0])
+            ]
+
+        # 6 beam focus point
+        l2_list = {
+            'F0G0': np.array([0,0,KBS_F]),
+            'F0G1': np.array([KBS_F*np.arctan(conf['theta_grating']),0,KBS_F]),
+            'F0G-1': np.array([-KBS_F*np.arctan(conf['theta_grating']),0,KBS_F]),
+            'F1G0': np.array([0,offaxis,F]),
+            'F1G1': np.array([F*np.arctan(conf['theta_grating']),offaxis,F]),
+            'F1G-1': np.array([-F*np.arctan(conf['theta_grating']),offaxis,F])
+            }
+
+        for beam in l2_list.keys():
+            l2 = l2_list[beam]
+            corners = []
+            for l1 in l1_list:
+                corners.append(self.LinePlaneIntersection(l1, l2, p0, n)[:2])
+            Beams[beam].set_xy(self.scale*np.array(corners))
 
     def DetectorUpdate(self, Xoff, Yoff):
         """ Draw DSSC detector module with filter mask.
@@ -199,7 +275,8 @@ class TZPGcalc():
         # moving rotated rectangles is a pain in matplotlib
         self.detLines['diamond'].set_xy((self.scale*Xoff, self.scale*(Yoff - diamondW/2*np.sqrt(2))))
 
-    def SampleUpdate(self, w, p, Xoff, Yoff, thickness=0.525):
+    def SampleUpdate(self, w, p, Xoff, Yoff, thickness=0.525,
+                     incidence=0, etch_angle=54.74):
         """ Draw the sample.
 
             w: membrane width
@@ -207,15 +284,26 @@ class TZPGcalc():
             Xoff: sample x offset
             Yoff: sample y offset
             thickness: sample thickness used to calculate the etched facets
+            incidence: incidence angle in rad
+            etch_angle: etching angle from surface in rad
         """
         # Si etching angle
-        wp = w +2*thickness/np.tan(np.deg2rad(54.74))
+        wp = w +2*thickness/np.tan(etch_angle)
+
+        # incidence angle squeezes sample and etch lines
+        # and induces an apparent shift off the etch lines
+        ci = np.cos(incidence)
+        thsi = thickness*np.sin(incidence)
 
         j = 0
         for k in range(-(SampleN-1)//2, (SampleN-1)//2+1):
             for l in range(-(SampleN-1)//2, (SampleN-1)//2+1):
-                self.RectUpdate(self.sampleLines[j], k*p - w/2 + Xoff, l*p - w/2 - Yoff, k*p + w/2 + Xoff, l*p + w/2 - Yoff)
-                self.RectUpdate(self.etchLines[j], k*p - wp/2 + Xoff, l*p - wp/2 - Yoff, k*p + wp/2 + Xoff, l*p + wp/2 - Yoff)
+                self.RectUpdate(self.sampleLines[j],
+                    ci*(k*p - w/2 + Xoff), l*p - w/2 - Yoff,
+                    ci*(k*p + w/2 + Xoff), l*p + w/2 - Yoff)
+                self.RectUpdate(self.etchLines[j],
+                    ci*(k*p - wp/2 + Xoff)+thsi, l*p - wp/2 - Yoff,
+                    ci*(k*p + wp/2 + Xoff)+thsi, l*p + wp/2 - Yoff)
                 j+=1
 
     def UpdateFig(self):
@@ -232,6 +320,7 @@ class TZPGcalc():
 
         theta_grating = self.grating_slider.value*1e-3 # [rad]
         sampleZ = self.samz_slider.value*1e-3 # [m]
+        samIncidence = np.deg2rad(self.samIncidence_slider.value) # [rad]
         detectorZ = self.det_slider.value*1e-3 # [m]
         TZPGwH = self.TZPGwH_slider.value*1e-3 #[m]
         TZPGwV = self.TZPGwV_slider.value*1e-3 #[m]
@@ -253,10 +342,10 @@ class TZPGcalc():
                  'TZPGwH':TZPGwH, 'TZPGwV':TZPGwV,  'TZPGo':TZPGo}
 
         # update the beams
-        self.UpdateBeams(self.samBeamsL, Z0 + sampleZ, confL)
-        self.UpdateBeams(self.detBeamsL, Z0 + detectorZ, confL)
-        self.UpdateBeams(self.samBeamsH, Z0 + sampleZ, confH)
-        self.UpdateBeams(self.detBeamsH, Z0 + detectorZ, confH)
+        self.UpdateBeams(self.samBeamsL, Z0 + sampleZ, samIncidence, confL)
+        self.UpdateBeams(self.detBeamsL, Z0 + detectorZ, 0, confL)
+        self.UpdateBeams(self.samBeamsH, Z0 + sampleZ, samIncidence, confH)
+        self.UpdateBeams(self.detBeamsH, Z0 + detectorZ, 0, confH)
 
         # update the detector
         detXoff = self.detX_slider.value*1e-3 #[m]
@@ -269,7 +358,9 @@ class TZPGcalc():
         samXoff = self.samX_slider.value*1e-3 #[m]
         samYoff = self.samY_slider.value*1e-3 #[m]
         samthickness = self.samthickness_slider.value*1e-6 #[m]
-        self.SampleUpdate(samw, samp, samXoff, samYoff, samthickness)
+        samEtchAngle = np.deg2rad(self.samEtchAngle_slider.value) #[rad]
+        self.SampleUpdate(samw, samp, samXoff, samYoff, samthickness,
+            samIncidence, samEtchAngle)
 
     def initWidgets(self):
         """ Creates the necessary interactive widget controls.
@@ -379,12 +470,28 @@ class TZPGcalc():
             step=1,
             readout_format='.0f',
         )
+        self.samIncidence_slider = widgets.FloatSlider(
+            value=0,
+            min=0,
+            max=90,
+            step=1,
+            readout_format='.0f',
+        )
+        self.samEtchAngle_slider = widgets.FloatSlider(
+            value=54.74,
+            min=0,
+            max=90,
+            step=0.01,
+            readout_format='.2f',
+        )
         samTab = VBox(children=[HBox([widgets.Label(value='Sample Z (mm):'), self.samz_slider]),
                                 HBox(children=[HBox([widgets.Label(value='Membrane width (mm):'), self.samw_slider]),
                                                HBox([widgets.Label(value='Membrane pitch (mm):'), self.samp_slider])]),
                                 HBox(children=[HBox([widgets.Label(value='Sample X-Offset (mm):'), self.samX_slider]),
                                                HBox([widgets.Label(value='Sample Y-Offset (mm):'), self.samY_slider])]),
-                                HBox([widgets.Label(value='Substrate thickness (um):'), self.samthickness_slider])
+                                HBox([widgets.Label(value='Substrate thickness (um):'), self.samthickness_slider]),
+                                HBox([HBox([widgets.Label(value='Normal incidence (deg):'), self.samIncidence_slider]),
+                                    HBox([widgets.Label(value='Etch angle from surface (deg):'), self.samEtchAngle_slider])])
                                 ])
 
         #detector tab