diff --git a/src/TZPGcalc/GeoBeams.py b/src/TZPGcalc/GeoBeams.py
index f0fc76a659bf3eda1a90a651af9dd021a68badd6..e1e975256ce4c4062bf1632843cbe1644e132205 100644
--- a/src/TZPGcalc/GeoBeams.py
+++ b/src/TZPGcalc/GeoBeams.py
@@ -30,7 +30,8 @@ class GeoBeams:
             'dVFM': -2,  # Vertical Focusing Mirror position
             'fVFM': 5.05,  # VFM focal length
             'dBOZ': -0.23,  # Beam splitting Off axis Zone plate position
-            'fBOZ': 0.25,  # BOZ focal length
+            'fBOZ_x': 0.25,  # BOZ horizontal focal length
+            'fBOZ_y': 0.25,  # BOZ vertical focal distance
             'theta_grating': 0,  # grating deflection angle in rad
             'dSAMZ': 0.0,  # Sample position
             'WH': 0.8e-3,  # BOZ horizontal width
@@ -55,7 +56,7 @@ class GeoBeams:
             position x2, given the source position at x1.
         """
         dEX, dIHF, dVFM, dHFM, dBOZ = symbols('dEX, dIHF, dVFM, dHFM, dBOZ')
-        fVFM, fHFM, fBOZ = symbols('fVFM, fHFM, fBOZ')
+        fVFM, fHFM, fBOZ_x, fBOZ_y = symbols('fVFM, fHFM, fBOZ_x, fBOZ_y')
         x1, theta1_x, x2 = symbols('x1, theta1_x, x2')
         y1, theta1_y, y2 = symbols('y1, theta1_y, y2')
 
@@ -71,8 +72,8 @@ class GeoBeams:
         self.VFM_f = lambdify(self.elems.keys(), self.VFM, ['numpy'])
 
         # beams for the zone plate first order (0th order is unfocused)
-        self.HBOZ = simplify(ThinLens(fBOZ)*self.HFM)
-        self.VBOZ = simplify(ThinLens(fBOZ)*self.VFM)
+        self.HBOZ = simplify(ThinLens(fBOZ_x)*self.HFM)
+        self.VBOZ = simplify(ThinLens(fBOZ_y)*self.VFM)
         self.HBOZ_f = lambdify(self.elems.keys(), self.HBOZ, ['numpy'])
         self.VBOZ_f = lambdify(self.elems.keys(), self.VBOZ, ['numpy'])
 
@@ -105,7 +106,7 @@ class GeoBeams:
         x += [ray[0].evalf()]
 
         # BOZ
-        ray = (ThinLens(self.elems['fBOZ']) *
+        ray = (ThinLens(self.elems['fBOZ_x']) *
                FreeSpace(-self.elems['dHFM'] + self.elems['dBOZ'])*ray)
         z_x += [self.elems['dBOZ']]
         x += [ray[0].evalf()]
@@ -132,7 +133,7 @@ class GeoBeams:
         y += [ray[0].evalf()]
 
         # BOZ
-        ray = (ThinLens(self.elems['fBOZ']) *
+        ray = (ThinLens(self.elems['fBOZ_y']) *
                FreeSpace(-self.elems['dVFM']+self.elems['dBOZ'])*ray)
         z_y += [self.elems['dBOZ']]
         y += [ray[0].evalf()]
diff --git a/src/TZPGcalc/TZPGcalc.py b/src/TZPGcalc/TZPGcalc.py
index ff0fb2c6885fef26a8fb39c7216428825a2e368c..cf6a17d1ed24101938849d350e74e156864cbe4d 100644
--- a/src/TZPGcalc/TZPGcalc.py
+++ b/src/TZPGcalc/TZPGcalc.py
@@ -10,8 +10,9 @@
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
-from matplotlib.patches import Rectangle, Polygon
+from matplotlib.patches import PathPatch, Rectangle, Polygon
 from matplotlib.colors import hsv_to_rgb
+from matplotlib.path import Path
 
 import ipywidgets as widgets
 from ipywidgets import HBox, VBox
@@ -31,7 +32,9 @@ TZPG_db = {
         'TZPGwV': 1,
         'TZPGoffaxis': 0.75,
         'grating': 3.8,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'O': {
         'design_nrj': 530,
@@ -39,7 +42,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'Fe': {
         'design_nrj': 715,
@@ -47,7 +52,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'Co': {
         'design_nrj': 785,
@@ -55,7 +62,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'Ni': {
         'design_nrj': 860,
@@ -63,7 +72,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'Cu': {
         'design_nrj': 927,
@@ -71,7 +82,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     },
     'Gd': {
         'design_nrj': 1210,
@@ -79,7 +92,9 @@ TZPG_db = {
         'TZPGwV': 0.8,
         'TZPGoffaxis': 0.55,
         'grating': 3.1,
-        'F': 0.25
+        'F_x': 0.25,
+        'F_y': 0.25,
+        '3beams': True
     }
 }
 
@@ -87,14 +102,9 @@ TZPG_db = {
 class TZPGcalc():
     def __init__(self):
         self.geo_beams = GeoBeams()
-        self.initFig()
         self.initWidgets()
-        for v in ['fVFM', 'fHFM', 'EXw', 'IHFw']:
-            if v in ['fVFM', 'fHFM']:
-                scale = 1.0
-            else:
-                scale = 1e6
-            self.widgets[v].value = scale*self.geo_beams.elems[v]
+        self.initFig()
+        self.init_beam_transport()
 
         # spot sizes of all beams
         self.SpotSizes = {}
@@ -106,6 +116,16 @@ class TZPGcalc():
         self.UpdateFig()
         display(self.control)
 
+    def init_beam_transport(self):
+        temp = GeoBeams()
+        # set default value for beam transport
+        for v in ['fVFM', 'fHFM', 'EXw', 'IHFw']:
+            if v in ['fVFM', 'fHFM']:
+                scale = 1.0
+            else:
+                scale = 1e6
+            self.widgets[v].value = scale*temp.elems[v]
+
     def initFig(self):
         "Creates a figure for the sample plane and detector plane images."
 
@@ -204,7 +224,7 @@ class TZPGcalc():
                           angle=45))
                }
 
-        # 5x5 membranes
+        # SampleNxSampleN membranes
         self.sampleLines = {}
         self.etchLines = {}
         for k in range(SampleN*SampleN):
@@ -214,6 +234,63 @@ class TZPGcalc():
                 Rectangle((0, 0), 1, 1, fill=False, facecolor='k',
                           alpha=0.4, ls='--'))
 
+        # Flat Liquid Jet
+        self.FLJ_lines = {}
+        self.FlatLiquidJet()
+
+    def FlatLiquidJet(self):
+        """Draw a Flat Liquid jet.
+        """
+        sw = self.widgets
+        HW = 0.5*sw['FLJ_W'].value  # [mm]
+        L = sw['FLJ_L'].value  # [mm]
+        mf = sw['FLJ_mf'].value  # [mm]
+        incidence = np.deg2rad(sw['samIncidence'].value)  # [rad]
+        ox = sw['samX'].value  # [mm]
+        oy = -L/2 + sw['samY'].value  # [mm]
+
+        # incidence angle squeezes sample
+        ci = np.cos(incidence)
+
+        verts = [
+            [(0*ci, L + oy),   # P0
+             (HW/2*ci, L + oy),  # P1
+             (HW*ci, 0.5*(1+mf)*L + oy),  # P2
+             (HW*ci, mf*L + oy)],  # P3
+            [(HW*ci, mf*L + oy),  # P0
+             (HW*ci, 0.5*mf*L + oy),  # P1
+             (0, 0 + oy),  # P2
+             (0, 0 + oy)]  # P3
+        ]
+
+        # second half image of liquid jet
+        mirror = []
+        for v in verts:
+            mirror.append([(-p[0], p[1]) for p in v])
+
+        verts += mirror
+
+        # apply offsets
+        fverts = []
+        for v in verts:
+            fverts.append([(p[0] + ox, p[1]) for p in v])
+
+        codes = [
+            Path.MOVETO,  # P0
+            Path.CURVE4,  # P1
+            Path.CURVE4,  # P2
+            Path.CURVE4  # P3
+        ]
+
+        for k, v in enumerate(fverts):
+            path = Path(v, codes)
+            if k in self.FLJ_lines:
+                self.FLJ_lines[k].remove()
+
+            self.FLJ_lines[k] = self.ax_sam.add_patch(
+                PathPatch(path, alpha=0.4, facecolor='none', lw=2)
+                )
+
     def RectUpdate(self, rect, xLeft, yBottom, xRight, yTop):
         """Updates the position and size of the given Rectangle.
 
@@ -266,7 +343,8 @@ class TZPGcalc():
         # shortcut
         sge = self.geo_beams.elems
 
-        sge['fBOZ'] = conf['F']
+        sge['fBOZ_x'] = conf['F_x']
+        sge['fBOZ_y'] = conf['F_y']
         sge['theta_grating'] = conf['theta_grating']
 
         # imaging plane
@@ -290,6 +368,11 @@ class TZPGcalc():
                 self.SpotSizes[img['type']][conf['Energy']][z, k] = (
                     1e3*(np.max(vs) - np.min(vs)))
 
+        # 3 beams configuration
+        b = self.widgets['3beams'].value
+        Beams['F0G0'].set_visible(b)
+        Beams['F1G0'].set_visible(b)
+
     def DetectorUpdate(self, Xoff, Yoff):
         """Draws DSSC detector module with filter mask.
 
@@ -325,8 +408,28 @@ class TZPGcalc():
         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,
-                     incidence=0, etch_angle=54.74):
+    def SampleUpdate(self):
+        if self.widgets['SampleType'].value == 'Membranes Array':
+            b = True
+        elif self.widgets['SampleType'].value == 'Flat Liquid Jet':
+            b = False
+        else:
+            raise ValueError('Sample type must be either "Membranes Array" or'
+                             '"Flat Liquid Jet"')
+
+        for k in range(SampleN*SampleN):
+            self.sampleLines[k].set_visible(b)
+            self.etchLines[k].set_visible(b)
+
+        for k in range(4):
+            self.FLJ_lines[k].set_visible(not(b))
+
+        if b:
+            self.MembraneSampleUpdate()
+        else:
+            self.FlatLiquidJet()
+
+    def MembraneSampleUpdate(self):
         """Draws the sample.
 
         Inputs
@@ -339,6 +442,15 @@ class TZPGcalc():
             incidence: incidence angle in rad
             etch_angle: etching angle from surface in rad
         """
+        sw = self.widgets
+        w = sw['samw'].value*1e-3  # [m]
+        p = sw['samp'].value*1e-3  # [m]
+        Xoff = sw['samX'].value*1e-3  # [m]
+        Yoff = sw['samY'].value*1e-3  # [m]
+        thickness = sw['samthickness'].value*1e-6  # [m]
+        etch_angle = np.deg2rad(sw['samEtchAngle'].value)  # [rad]
+        incidence = np.deg2rad(sw['samIncidence'].value)  # [rad]
+
         # Si etching angle
         wp = w + 2*thickness/np.tan(etch_angle)
 
@@ -371,12 +483,13 @@ class TZPGcalc():
         # we calculate the optics for the central wavelength
         nrjL = sw['nrjL'].value  # [eV]
         nrjH = sw['nrjH'].value  # [eV]
+        nrjD = sw['nrjD'].value  # [eV]
         wlL = 1240/nrjL*1e-9
         wlH = 1240/nrjH*1e-9
-        nrjD = sw['nrjD'].value  # [eV]
         wlD = 1240/nrjD*1e-9
 
-        F = sw['F'].value  # [m] hidden nominal focal length
+        F_x = sw['F_x'].value  # [m] Nominal horiz. BOZ focal length
+        F_y = sw['F_y'].value  # [m] Nominal vert. BOZ focal length
         theta_grating = sw['grating'].value*1e-3  # [rad]
         sampleZ = sw['SAMZ'].value*1e-3  # [m]
         samIncidence = np.deg2rad(sw['samIncidence'].value)  # [rad]
@@ -397,21 +510,30 @@ class TZPGcalc():
         # zone plate radius at the point further away from the optical axis
         rn = np.sqrt((sge['WH']/2.0)**2 +
                      (sge['WV']/2.0 + np.abs(sge['offaxis']))**2)
-        dr_nominal = wlD * F / (2*rn)
-        sw['dr_label'].value = (
-            f'Outer Zone Plate width dr:{int(np.round(dr_nominal*1e9))} nm')
+        dr_nominal_x = wlD * F_x / (2*rn)
+        dr_nominal_y = wlD * F_y / (2*rn)
+        sw['dr_label_x'].value = (
+            f'Outer Zone Plate width dr for Horiz. '
+            f'focus:{int(np.round(dr_nominal_x*1e9))} nm')
+        sw['dr_label_y'].value = (
+            f'for Vert. focus:{int(np.round(dr_nominal_y*1e9))} nm')
 
         # Optics properties (focal length and grating angle) for the
         # low energy and high energy photon
-        F_L = 2*rn*dr_nominal/wlL
+        F_L_x = 2*rn*dr_nominal_x/wlL
+        F_L_y = 2*rn*dr_nominal_y/wlL
         G_L = np.arcsin(wlL/d_nominal)
         confL = {'Energy': 'L',
-                 'F': F_L,
+                 'F_x': F_L_x,
+                 'F_y': F_L_y,
                  'theta_grating': G_L}
-        F_H = 2*rn*dr_nominal/wlH
+
+        F_H_x = 2*rn*dr_nominal_x/wlH
+        F_H_y = 2*rn*dr_nominal_y/wlH
         G_H = np.arcsin(wlH/d_nominal)
         confH = {'Energy': 'H',
-                 'F': F_H,
+                 'F_x': F_H_x,
+                 'F_y': F_H_y,
                  'theta_grating': G_H}
 
         # Configuration for imaging plane
@@ -439,14 +561,7 @@ class TZPGcalc():
         self.DetectorUpdate(detXoff, detYoff)
 
         # update the sample
-        samw = self.widgets['samw'].value*1e-3  # [m]
-        samp = self.widgets['samp'].value*1e-3  # [m]
-        samXoff = self.widgets['samX'].value*1e-3  # [m]
-        samYoff = self.widgets['samY'].value*1e-3  # [m]
-        samthickness = self.widgets['samthickness'].value*1e-6  # [m]
-        samEtchAngle = np.deg2rad(self.widgets['samEtchAngle'].value)  # [rad]
-        self.SampleUpdate(samw, samp, samXoff, samYoff, samthickness,
-                          samIncidence, samEtchAngle)
+        self.SampleUpdate()
 
     def initWidgets(self):
         """ Creates the necessary interactive widget controls.
@@ -476,6 +591,14 @@ class TZPGcalc():
                          ])
 
         # Source
+        self.Reset = widgets.Button(
+            description='Reset',
+        )
+
+        @self.Reset.on_click
+        def reset_on_click(b):
+            self.init_beam_transport()
+
         self.widgets['EXw'] = widgets.BoundedIntText(
             value=100,
             min=0,
@@ -513,6 +636,7 @@ class TZPGcalc():
             layout=layout
         )
         SourceTab = VBox([
+            self.Reset,
             self.widgets['EXw'],
             self.widgets['IHFw'],
             self.widgets['fVFM'],
@@ -534,17 +658,33 @@ class TZPGcalc():
             self.widgets['TZPGwV'].value = v['TZPGwV']
             self.widgets['TZPGoffaxis'].value = v['TZPGoffaxis']
             self.widgets['grating'].value = v['grating']
-            self.widgets['F'].value = v['F']
+            self.widgets['F_x'].value = v['F_x']
+            self.widgets['F_y'].value = v['F_y']
+            self.widgets['3beams'].value = v['3beams']
+
             # necessary to recompute grating pitch and outer zone plate width
             self.UpdateFig()
 
         self.widgets['Type'].observe(TZPGtype, names='value')
 
         # hidden nominal zone plate focus
-        self.widgets['F'] = widgets.BoundedFloatText(
+        self.widgets['F_x'] = widgets.BoundedFloatText(
             value=0.25,
             min=0,
-            max=1
+            max=1,
+            step=0.01,
+            description='Focal length (m) Horiz.:',
+            style=style,
+            layout=layout
+        )
+        self.widgets['F_y'] = widgets.BoundedFloatText(
+            value=0.25,
+            min=0,
+            max=1,
+            step=0.01,
+            description='Vert.:',
+            style=style,
+            layout=layout
         )
 
         self.widgets['nrjL'] = widgets.BoundedIntText(
@@ -571,7 +711,7 @@ class TZPGcalc():
             max=3200,
             step=1,
             width=4,
-            description='Design:',
+            description='Design energy (eV):',
             style=style,
             layout=layout
         )
@@ -611,15 +751,26 @@ class TZPGcalc():
             style=style,
             layout=layout
         )
-        self.widgets['dr_label'] = widgets.Label(value='dr')
+        self.widgets['3beams'] = widgets.Checkbox(
+            value=True,
+            description='3 beams:',
+            style=style,
+            layout=layout
+        )
+        self.widgets['dr_label_x'] = widgets.Label(value='dr_x')
+        self.widgets['dr_label_y'] = widgets.Label(value='dr_y')
         self.widgets['d_label'] = widgets.Label(value='dr')
         TZPGTab = VBox([
             self.widgets['Type'],
             HBox([self.widgets['nrjD'],
-                  self.widgets['grating'],
-                  self.widgets['TZPGoffaxis']]),
-            HBox([self.widgets['d_label'],
-                  self.widgets['dr_label']]),
+                  self.widgets['F_x'],
+                  self.widgets['F_y']
+                  ]),
+            HBox([self.widgets['grating'],
+                  self.widgets['TZPGoffaxis'],
+                  self.widgets['3beams']]),
+            HBox([self.widgets['dr_label_x'],
+                  self.widgets['dr_label_y']]),
             HBox([widgets.Label(value='Energy range (eV):'),
                   self.widgets['nrjL'], self.widgets['nrjH']]),
             HBox([widgets.Label(value='TZPG (mm):'),
@@ -628,6 +779,13 @@ class TZPGcalc():
             ])
 
         # sample part
+        self.widgets['SampleType'] = widgets.Dropdown(
+            options=['Membranes Array', 'Flat Liquid Jet'],
+            value='Membranes Array',
+            decription='Sample type:',
+            style=style,
+            layout=layout
+        )
         self.widgets['SAMZ'] = widgets.BoundedFloatText(
             value=30.,
             min=-10.,
@@ -637,12 +795,14 @@ class TZPGcalc():
             style=style,
             layout=layout
         )
+
+        # membranes
         self.widgets['samw'] = widgets.BoundedFloatText(
             value=.5,
             min=0.01,
             max=2.0,
             step=.01,
-            description='width:',
+            description='width (mm):',
             style=style,
             layout=layout
         )
@@ -651,10 +811,58 @@ class TZPGcalc():
             min=0.01,
             max=2.0,
             step=.01,
-            description='pitch:',
+            description='pitch (mm):',
+            style=style,
+            layout=layout
+        )
+        self.widgets['samthickness'] = widgets.BoundedFloatText(
+            value=381,
+            min=1,
+            max=1000,
+            step=1,
+            description='Substrate thickness (um):',
+            style=style,
+            layout=layout
+        )
+        self.widgets['samEtchAngle'] = widgets.BoundedFloatText(
+            value=54.74,
+            min=0,
+            max=90,
+            step=0.01,
+            description='etch angle (deg):',
+            style=style,
+            layout=layout
+        )
+
+        # Flat Liquid Jet
+        self.widgets['FLJ_W'] = widgets.BoundedFloatText(
+            value=1.0,
+            min=0,
+            max=5,
+            step=0.1,
+            description='Flat Liquid Jet width (mm):',
+            style=style,
+            layout=layout
+        )
+        self.widgets['FLJ_L'] = widgets.BoundedFloatText(
+            value=4.6,
+            min=0,
+            max=15,
+            step=0.1,
+            description='of length (mm):',
             style=style,
             layout=layout
         )
+        self.widgets['FLJ_mf'] = widgets.BoundedFloatText(
+            value=0.75,
+            min=0,
+            max=1,
+            step=0.01,
+            description='at:',
+            style=style,
+            layout=layout
+        )
+
         self.widgets['samX'] = widgets.BoundedFloatText(
             value=0.,
             min=-10,
@@ -673,15 +881,6 @@ class TZPGcalc():
             style=style,
             layout=layout
         )
-        self.widgets['samthickness'] = widgets.BoundedFloatText(
-            value=381,
-            min=1,
-            max=1000,
-            step=1,
-            description='Substrate thickness (um):',
-            style=style,
-            layout=layout
-        )
         self.widgets['samIncidence'] = widgets.BoundedFloatText(
             value=0,
             min=0,
@@ -691,26 +890,22 @@ class TZPGcalc():
             style=style,
             layout=layout
         )
-        self.widgets['samEtchAngle'] = widgets.BoundedFloatText(
-            value=54.74,
-            min=0,
-            max=90,
-            step=0.01,
-            description='Etch angle from surface (deg):',
-            style=style,
-            layout=layout
-        )
         samTab = VBox([
             self.widgets['SAMZ'],
-            HBox([widgets.Label(value='Membrane (mm), '),
+            self.widgets['SampleType'],
+            HBox([widgets.Label(value='Membranes array, '),
                   self.widgets['samw'],
                   self.widgets['samp']]),
+            HBox([self.widgets['samthickness'],
+                  self.widgets['samEtchAngle']]),
+            HBox([self.widgets['FLJ_W'],
+                  self.widgets['FLJ_mf'],
+                  self.widgets['FLJ_L']]),
             HBox([widgets.Label(value='Sample Offset (mm), '),
                   self.widgets['samX'],
                   self.widgets['samY']]),
-            self.widgets['samthickness'],
             HBox([self.widgets['samIncidence'],
-                  self.widgets['samEtchAngle']])
+                  ])
         ])
 
         # Detector tab
@@ -724,7 +919,7 @@ class TZPGcalc():
             layout=layout
         )
         self.widgets['detX'] = widgets.BoundedFloatText(
-            value=20.,
+            value=34.5,
             min=-50,
             max=50,
             step=0.5,
@@ -733,7 +928,7 @@ class TZPGcalc():
             layout=layout
         )
         self.widgets['detY'] = widgets.BoundedFloatText(
-            value=0.,
+            value=-2.,
             min=-50,
             max=50,
             step=0.5,