diff --git a/zeroheliumkit/src/__init__.py b/zeroheliumkit/src/__init__.py
index b29618c..471cdc0 100755
--- a/zeroheliumkit/src/__init__.py
+++ b/zeroheliumkit/src/__init__.py
@@ -1,5 +1,5 @@
 from .anchors import Anchor, MultiAnchor, Skeletone, Layer
-from .core import Entity, Structure, GeomCollection
+from .core import Entity, Structure, GeomCollection,ReferenceStructure
 from .supercore import SuperStructure, ContinuousLineBuilder, RoutingConfig, ObjsAlongConfig
 
 from .geometries import (StraightLine, ArbitraryLine, Taper, Fillet,
@@ -17,6 +17,7 @@
     "Structure",
     "GeomCollection",
     "SuperStructure",
+    "ReferenceStructure",
     "ContinuousLineBuilder",
     "Rectangle", "Square", "Circle", "RegularPolygon", "ArcLine",
     "Meander", "MeanderHalf", "PinchGate",
diff --git a/zeroheliumkit/src/core.py b/zeroheliumkit/src/core.py
index 4be9740..6f83bb6 100755
--- a/zeroheliumkit/src/core.py
+++ b/zeroheliumkit/src/core.py
@@ -10,15 +10,17 @@
 """
 
 import copy
+import math
 import matplotlib.pyplot as plt
 from warnings import warn
 
 from shapely import (Point, MultiPoint, LineString, MultiLineString,
                      Polygon, MultiPolygon, GeometryCollection)
+import gdstk
 
 from .plotting import interactive_widget_handler, listify_colors, ColorHandler
-from .importing import Exporter_DXF, Exporter_GDS, Exporter_Pickle
-from .settings import SIZE, SIZE_L, SIZE_S, RED, DARKGRAY
+from .importing import Exporter_DXF, Exporter_GDS, Exporter_Pickle,Reader_GDS
+from .settings import SIZE, SIZE_L, SIZE_S, RED, DARKGRAY, BLACK
 from .anchors import Anchor, MultiAnchor, Skeletone, Layer, get_dxdy
 from .errors import hard_deprecated
 
@@ -38,6 +40,7 @@ class Entity():
         skeletone (Skeletone): Represents a collection of lines linked to the Entity,
             which is an instance of the Skeletone class.
         anchors (MultiAnchor): Represents the anchor points of the Entity, which are instances of MultiAnchor.
+        library: gdstk library containing information about nested cells
 
     """
     layers = []
@@ -369,7 +372,7 @@ def export_pickle(self, filename: str) -> None:
         exp.save()
 
 
-    def export_gds(self, filename: str, layer_cfg: dict) -> None:
+    def export_gds(self, filename: str, layer_cfg: dict,cellname:str='toplevel') -> None:
         """
         Exports all layers as a GDS file.
 
@@ -379,7 +382,7 @@ def export_gds(self, filename: str, layer_cfg: dict) -> None:
                 See `gdspy docs <https://gdspy.readthedocs.io/en/stable/gettingstarted.html#layer-and-datatype>`_ for 'datatype' details.
         """
         zhkdict = self.export_dict(remove_holes=True)
-        exp = Exporter_GDS(filename, zhkdict, layer_cfg)
+        exp = Exporter_GDS(filename, zhkdict, layer_cfg,cellname,library=None)
         exp.save()
 
 
@@ -395,7 +398,6 @@ def export_dxf(self, filename: str, layer_cfg: list) -> None:
         exp = Exporter_DXF(filename, zhkdict, layer_cfg)
         exp.save()
 
-
     #############################
     #### Plotting operations ####
     #############################
@@ -625,3 +627,374 @@ def __init__(self, layers: dict=None):
 
         if self.colors.is_empty:
             self.colors.update_colors(self.layers)
+
+
+class ReferenceStructure(Structure):
+    """
+    Represents a structure that contains layers with a collection of geometries (Points, LineStrings, Polygons, etc.).
+    The ReferenceStructure class provides methods to work with nested gdstk libraries,cells, and references. 
+    Inherits from the Entity class.
+    """
+
+    library = gdstk.Library()
+
+    def __init__(self,cellname='toplevel'):
+        super().__init__()
+        self.topCellName = cellname
+        self.library = gdstk.Library()
+        self.topcell = self.library.new_cell(cellname)
+        self.cellNames = [cell.name for cell in self.library.cells]
+
+        
+
+    ##############################
+    #### reference operations ####
+    ##############################
+
+    def add_reference(self,referenceCell,Coord:tuple=(0,0),rotation:float=0):
+        '''
+        Adds a 2d array of references to the structures library (visibile upon export)
+        
+        Args:
+            referenceCell (gdstk.cell or gds path): the cell to be referenced and arrayed
+            Coord (tuple): (x,y) pair for the center of the first instance of the array
+            rotation (float): rotation of the reference in degrees
+        '''
+
+        if isinstance(referenceCell,str):
+            A = Reader_GDS(referenceCell,verbose=False)
+            topCellName = A.gdsii.top_level()[0].name
+            topcell = A.gdsii[topCellName]
+            referenceCell = topcell
+
+        self.update_dependencies(referenceCell)
+
+        self.topcell.add(gdstk.Reference(referenceCell,Coord,rotation=rotation*3.1415926535/180))
+
+    def add_reference_array(self,referenceCell,initCoord:tuple=(0,0),columns:int=1,rows:int=1,spacing:tuple=(0,0),rotation:float=0):
+        '''
+        Adds a 2d array of references to the structures library (visibile upon export)
+        
+        Args:
+            referenceCell (gdstk.cell or gds path): the cell to be referenced and arrayed
+            initCoord (tuple): (x,y) pair for the center of the first instance of the array
+            columns (int): number of columns of the array
+            rows (int): number of rows of the array
+            spacing (tuple): (dx,dx) spacing vector bewteen the columns and rows centerpoints
+            rotation (float): rotation of the references in degrees
+        '''
+        if isinstance(referenceCell,str):
+            A = Reader_GDS(referenceCell,verbose=False)
+            topCellName = A.gdsii.top_level()[0].name
+            topcell = A.gdsii[topCellName]
+            referenceCell = topcell
+
+        self.update_dependencies(referenceCell)
+
+        self.topcell.add(gdstk.Reference(referenceCell,initCoord,columns=columns,rows=rows,spacing=spacing,rotation=rotation*3.1415926535/180))
+
+    def update_dependencies(self,referenceCell):
+        # brings along any dependent cells that the new reference cell references
+        if referenceCell.name not in self.cellNames:
+            self.library.add(referenceCell)
+            self.cellNames.append(referenceCell.name)
+            for dependency in referenceCell.dependencies(True):
+                if dependency.name not in self.cellNames:
+                    self.library.add(dependency)
+                    self.cellNames.append(dependency.name)
+
+    def get_position_dict(self):
+        bounds = []
+
+        for lname in self.layers:
+            layer = getattr(self, lname, None)
+            if layer is None or layer.is_empty:
+                continue
+            bounds.append(layer.polygons.bounds)  # (xmin, ymin, xmax, ymax)
+
+        bb = self.library[self.topCellName].bounding_box()
+        if bb is not None:
+            (xmin, ymin), (xmax, ymax) = bb
+            bounds.append((xmin, ymin, xmax, ymax))
+
+        if not bounds:
+            raise ValueError(
+                "get_position_dict: ReferenceStructure has no polygons or "
+                "references to compute a bounding box from."
+            )
+
+        xmin = min(b[0] for b in bounds)
+        ymin = min(b[1] for b in bounds)
+        xmax = max(b[2] for b in bounds)
+        ymax = max(b[3] for b in bounds)
+
+        p_dict = {}
+        p_dict['xmin'] = xmin
+        p_dict['xmax'] = xmax
+        p_dict['ymin'] = ymin
+        p_dict['ymax'] = ymax
+        p_dict['dx'] = p_dict['xmax']-p_dict['xmin']
+        p_dict['dy'] = p_dict['ymax']-p_dict['ymin']
+        p_dict['x0'] = (p_dict['xmax']+p_dict['xmin'])/2
+        p_dict['y0'] = (p_dict['ymax']+p_dict['ymin'])/2
+        return p_dict
+
+    def move(self, dx: float, dy: float):
+        super().move(dx, dy)
+        for ref in self.topcell.references:
+            ox, oy = ref.origin
+            ref.origin = (ox + dx, oy + dy)
+        return self
+
+    def rotate(self, angle: float = 0, origin=(0, 0)):
+        super().rotate(angle, origin)
+        theta = math.radians(angle)
+        ox0, oy0 = origin
+        cos_t, sin_t = math.cos(theta), math.sin(theta)
+
+        def rotate_vec(v):
+            x, y = v
+            return (x * cos_t - y * sin_t, x * sin_t + y * cos_t)
+
+        for ref in self.topcell.references:
+            x, y = ref.origin
+            dx, dy = x - ox0, y - oy0
+            ref.origin = (ox0 + dx * cos_t - dy * sin_t, oy0 + dx * sin_t + dy * cos_t)
+            ref.rotation += theta
+
+            rep = ref.repetition
+            if rep is None:
+                continue
+
+            if rep.spacing is not None:
+                v1_local, v2_local = (rep.spacing[0], 0.0), (0.0, rep.spacing[1])
+            elif rep.v1 is not None and rep.v2 is not None:
+                v1_local, v2_local = rep.v1, rep.v2
+            else:
+                local_offsets = rep.get_offsets()
+                rotated = [rotate_vec(tuple(o)) for o in local_offsets]
+                ref.repetition = gdstk.Repetition(offsets=rotated[1:])
+                continue
+
+            ref.repetition = gdstk.Repetition(
+                columns=rep.columns, rows=rep.rows,
+                v1=rotate_vec(v1_local), v2=rotate_vec(v2_local),
+            )
+        return self
+
+    ##############################
+    #### Plotting operations  ####
+    ##############################
+
+    def quickplot_with_references(
+            self,
+            size="large",
+            color_config: dict=None,
+            zoom: tuple=None,
+            show_idx: bool=False,
+            off: list=[],
+            labels: bool=False,
+            draw_anchor_dir: bool=True,
+            ax=None,
+            library_mode: str="bbox",
+            library_color=None,
+            library_alpha: float=0.3,
+            library_labels: bool=True,
+            library_fontsize: float=8,
+            export_config: dict=None,
+            **kwargs
+            ) -> None:
+        """
+        ... (same docstring as before, plus:)
+
+        Args:
+            export_config (dict): the same {zhk_layer_name: {"layer": int,
+                "datatype": int}} dict you'd pass to export_gds. When given in
+                "polygons" mode, it's used to map each referenced polygon's GDS
+                (layer, datatype) back to a zhk layer name so it can be colored
+                from `color_config` exactly like quickplot() colors this
+                object's own layers. Ignored in "bbox" mode.
+        """
+        ax = self.quickplot(size=size, color_config=color_config, zoom=None,
+                            show_idx=show_idx, off=off, labels=labels,
+                            draw_anchor_dir=draw_anchor_dir, ax=ax, **kwargs)
+
+        if library_mode == "bbox":
+            self._plot_library_bboxes(ax, color=library_color or DARKGRAY,
+                                    alpha=library_alpha, labels=library_labels,
+                                    fontsize=library_fontsize)
+        elif library_mode == "polygons":
+            self._plot_library_polygons(ax, color=library_color, alpha=library_alpha,
+                                        color_config=color_config, export_config=export_config)
+        elif library_mode not in (None, False):
+            raise ValueError(f"Unknown library_mode {library_mode!r}, expected 'bbox', 'polygons', or None")
+
+        if zoom is not None:
+            xmin, xmax = ax.get_xlim()
+            ymin, ymax = ax.get_ylim()
+            x0, y0 = zoom[0]
+            dx = round((xmax - xmin) / zoom[1] / 2)
+            dy = round((ymax - ymin) / zoom[1] / 2)
+            if len(zoom) > 2:
+                dy = dy / zoom[2]
+            ax.set_xlim(x0 - dx, x0 + dx)
+            ax.set_ylim(y0 - dy, y0 + dy)
+
+        ax.set_aspect('equal')
+        return ax
+
+    def _iter_top_references(self):
+        """
+        Yields (cell, (x, y), rotation, magnification, x_reflection) for every
+        reference sitting directly on self.topcell, expanding array
+        repetitions (add_reference_array) into individual placements.
+        """
+        for ref in self.topcell.references:
+            cell = ref.cell
+            rotation = ref.rotation or 0.0
+            mag = ref.magnification or 1.0
+            x_refl = ref.x_reflection
+            ox, oy = ref.origin
+            rep = ref.repetition
+            if rep is None:
+                yield cell, (ox, oy), rotation, mag, x_refl
+            else:
+                for dx, dy in rep.get_offsets():
+                    yield cell, (ox + dx, oy + dy), rotation, mag, x_refl
+
+    def _plot_library_bboxes(self, ax, color=None, alpha=0.9, labels=True, fontsize=8):
+        """Draws a labeled, axis-aligned bounding box for every reference on the topcell."""
+        from matplotlib.patches import Rectangle
+
+        color = color or "magenta"
+        topcell = self.library[self.topCellName]   # look up by name, not via self.topcell
+
+        for ref in topcell.references:
+            bb = ref.bounding_box()          # gdstk handles rotation/mag/reflection/repetition for us
+            if bb is None:
+                continue  # empty referenced cell, nothing to draw
+            (xmin, ymin), (xmax, ymax) = bb
+
+            ax.add_patch(Rectangle((xmin, ymin), xmax - xmin, ymax - ymin,
+                                    fill=False, edgecolor=color, linestyle="--",
+                                    alpha=alpha, linewidth=1.2))
+
+            if labels:
+                ax.text((xmin + xmax) / 2, (ymin + ymax) / 2, ref.cell.name,
+                        ha="center", va="center", fontsize=fontsize, color=color)
+
+        ax.autoscale_view()
+
+    def _plot_library_polygons(self, ax, color=None, alpha=0.3, color_config=None, export_config=None):
+        """
+        Flattens and draws the actual geometry of every referenced cell (recursively).
+
+        Styling:
+            If `export_config` is given (the same {zhk_layer_name: {"layer": int,
+            "datatype": int}} dict passed to export_gds), each polygon's GDS
+            (layer, datatype) is mapped back to a zhk layer name, and that
+            layer's [color, alpha] is pulled from `color_config` - same dict you'd
+            pass to quickplot's color_config - or from self.colors if
+            color_config is None. This makes the referenced library geometry
+            match the styling of this object's own layers.
+            Any polygon whose (layer, datatype) isn't found in export_config
+            (or export_config isn't given) falls back to `color`/`alpha` if you
+            passed those explicitly, else a tab20 color cycle.
+
+            Every GDS (layer, datatype) that resolves to the same zhk layer name
+            (e.g. several via sub-layers all mapped to "vias") is merged into a
+            single unioned shape before drawing, so it's plotted once, flat -
+            exactly like quickplot() draws one of this object's own layers -
+            instead of as several overlapping, alpha-stacking patches. Layers are
+            then drawn bottom-to-top in `color_config`'s (or self.colors') order,
+            matching quickplot()'s layer ordering, and get the same black outline
+            quickplot() gives semi-transparent layers.
+        """
+        from collections import defaultdict
+        from shapely.ops import unary_union
+        from shapely.plotting import plot_polygon, plot_line
+
+        topcell = self.library[self.topCellName]   # look up by name, not via self.topcell
+
+        polys = topcell.get_polygons(apply_repetitions=True, include_paths=True, depth=None)
+
+        # (gds layer, datatype) -> zhk layer name
+        spec_to_name = {}
+        if export_config is not None:
+            for lname, spec in export_config.items():
+                spec_to_name[(spec["layer"], spec.get("datatype", 0))] = lname
+
+        # group by resolved zhk layer name (falling back to "L{layer}_{datatype}" for
+        # anything export_config doesn't know about) so every sub-layer sharing a name
+        # ends up in one flat shape rather than several stacked collections
+        by_name = defaultdict(list)
+        for p in polys:
+            lname = spec_to_name.get((p.layer, p.datatype), f"L{p.layer}_{p.datatype}")
+            by_name[lname].append(Polygon(p.points))
+
+        # zhk layer name -> [color, alpha], same source quickplot() itself uses
+        named_styles = listify_colors(color_config) if color_config else self.colors.colors
+
+        # bottom-to-top order: quickplot()'s layer order first, then any leftovers
+        ordered_names = [n for n in named_styles if n in by_name]
+        ordered_names += [n for n in by_name if n not in named_styles]
+
+        cmap = plt.get_cmap("tab20")
+        for i, lname in enumerate(ordered_names):
+            merged = unary_union(by_name[lname])
+            style = named_styles.get(lname)
+
+            if style is not None:
+                c, a = style[0], style[1]
+            elif color is not None:
+                c, a = color, alpha
+            else:
+                c, a = cmap(i % 20), alpha
+
+            plot_polygon(merged, ax=ax, color=c, alpha=a, edgecolor=BLACK,
+                        add_points=False, label=lname)
+            if a != 1:
+                plot_line(merged.boundary, ax=ax, color=BLACK, add_points=False, lw=1.5)
+        ax.autoscale_view()
+
+    ##############################
+    #### Exporting operations ####
+    ##############################
+
+    def export_gds(self, filename: str, layer_cfg: dict) -> None:
+        """
+        Exports all layers as a GDS file.
+
+        Args:
+            filename (str): The name of the gds file to be exported.
+            layer_cfg (dict): A dictionary containing the layer configuration.
+                See `gdspy docs <https://gdspy.readthedocs.io/en/stable/gettingstarted.html#layer-and-datatype>`_ for 'datatype' details.
+        """
+        zhkdict = self.export_dict(remove_holes=True)
+        exp = Exporter_GDS(filename, zhkdict, layer_cfg,self.topCellName,self.library)
+        exp.save()
+
+    ##############################
+    #### Importing operations ####
+    ##############################
+
+    def import_gds(self,componentFolder,componentName,export_config,plot_config=None):
+        # open the GDS file
+        A = Reader_GDS(f'{componentFolder}/{componentName}.gds',verbose=False)
+
+        # Add the lib to the struct
+        self.library = A.gdsii
+
+        # Update class attributes
+        self.cellNames = [cell.name for cell in self.library.cells]
+        self.topCellName = A.gdsii.top_level()[0].name
+        self.topcell = self.library[self.topCellName]
+        
+        # import the geometry from the top cell
+        cell = A.cells[self.topCellName]
+        for layerNumber,geom in cell.items():
+            layerName = [k for k,v in export_config.items() if v['layer']==layerNumber][0]
+            if plot_config is None:
+                self.add(Layer(layerName,geom))
+            else:
+                self.add(Layer(layerName,geom,plot_config[layerName]))
diff --git a/zeroheliumkit/src/importing.py b/zeroheliumkit/src/importing.py
index 4f8f9e4..dd65fb3 100755
--- a/zeroheliumkit/src/importing.py
+++ b/zeroheliumkit/src/importing.py
@@ -34,38 +34,48 @@ class Exporter_GDS():
                           Expected keys like {"layer": int, "datatype": int}.
     """
 
-    __slots__ = "name", "zhk_layers", "gdsii", "layer_cfg"
+    __slots__ = "name", "zhk_layers", "gdsii", "layer_cfg","cellname","library"
 
-    def __init__(self, name: str, zhk_layers: dict, layer_cfg: dict) -> None:
+    def __init__(self, name: str, zhk_layers: dict, layer_cfg: dict,cellname:str,library) -> None:
         self.name = name
         self.zhk_layers = zhk_layers
         self.layer_cfg = layer_cfg
+        self.cellname = cellname
+        self.library = library
         self.preapre_gds()
 
     def preapre_gds(self) -> None:
-        """
-        Prepare the GDSII library by creating a top-level cell and adding polygons.
-
-        Notes vs gdspy:
-        - gdstk does not have `exclude_from_current`; cells are not automatically "current".
-        - gdstk polygons use `layer` and `datatype` (same concepts).
-        """
-        self.gdsii = gdstk.Library()
-        cell = gdstk.Cell("toplevel")
-        self.gdsii.add(cell)
-
-        for lname, l_property in self.layer_cfg.items():
+        if self.library is None:
+            self.gdsii = gdstk.Library()
+            topcell = gdstk.Cell(self.cellname)
+            self.gdsii.add(topcell)
+        else:
+            self.gdsii = self.library
+            # Check if cell already exists in the library
+            existing_names = {c.name: c for c in self.gdsii.cells}
+            if self.cellname in existing_names:
+                # Reuse existing cell, just append geometry to it
+                topcell = existing_names[self.cellname]
+            else:
+                # Cell not found, create a new one and add it
+                topcell = gdstk.Cell(self.cellname)
+                self.gdsii.add(topcell)
+
+        # Add structure geometry to topcell (existing or new)
+        for lname in self.zhk_layers.keys():
+            if lname in ['skeletone', 'anchors']:
+                continue
+            l_property = self.layer_cfg[lname]
             polygons = self.zhk_layers[lname].polygons
             for poly in polygons.geoms:
                 points = list(poly.exterior.coords)
-
-                # Optional: shapely exterior repeats the first point at the end.
-                # gdstk is fine either way, but removing the duplicate keeps things tidy.
                 if len(points) > 1 and points[0] == points[-1]:
                     points = points[:-1]
-
                 gds_poly = gdstk.Polygon(points, **l_property)
-                cell.add(gds_poly)
+                topcell.add(gds_poly)
+
+
+        
 
     def save(self):
         """
@@ -89,17 +99,18 @@ class Reader_GDS():
         cells (dict): Dict mapping cellname -> {layer_number -> MultiPolygon}.
     """
 
-    __slots__ = "filename", "geometries", "gdsii", "cells"
+    __slots__ = "filename", "geometries", "gdsii", "cells","references"
 
-    def __init__(self, filename: str, cellname: str = "toplevel"):
+    def __init__(self, filename: str,verbose=True):
         self.filename = filename
         self.geometries = {}
         self.cells = {}
+        self.references = {}
         self.gdsii = gdstk.read_gds(filename)
-        self.extract_geometries()
-        self.prepare_dict(cellname)
+        self.extract_geometries(verbose)
+        self.prepare_dict()
 
-    def extract_geometries(self) -> None:
+    def extract_geometries(self,verbose) -> None:
         cells_out = {}
 
         # gdstk: library.cells is a list of Cell objects
@@ -126,7 +137,8 @@ def extract_geometries(self) -> None:
                 by_layer.setdefault(layer, []).append(shp)
 
             layer_numbers = sorted(by_layer.keys())
-            print(f"{self.filename} // Layers in cell '{name}': {layer_numbers}")
+            if verbose:
+                print(f"{self.filename} // Layers in cell '{name}': {layer_numbers}")
 
             # Build MultiPolygon per layer via unary_union
             layer_map = {}
@@ -139,16 +151,21 @@ def extract_geometries(self) -> None:
 
             cells_out[name] = layer_map
 
+            # collect the references
+            self.references[name] = cell.references
+                
+
+
         self.cells = cells_out
 
-    def prepare_dict(self, cellname: str = "toplevel") -> None:
+    def prepare_dict(self) -> None:
+        cellname = list(self.cells.keys())[0]
         geoms = self.cells[cellname]
         self.geometries = {
             ("L" + str(k) if isinstance(k, numbers.Number) else k): v
             for k, v in geoms.items()
         }
 
-
 class Exporter_DXF():
     """
     Helper class to export zhk dictionary with geometries into .dxf file.