diff --git a/README.md b/README.md
index 60cbebfd..0f65b7cb 100644
--- a/README.md
+++ b/README.md
@@ -17,4 +17,4 @@ The purpose of Paramak is to provide geometry for parameter studies and therefor
 
 Paramak produces [CadQuery](https://github.com/CadQuery/CadQuery) objects which can easily be saved to CAD files in several formats (STEP, STL, BRep) and also readily convertable to DAGMC h5m neutronics geometry using tools such as [Cad-to_DAGMC](https://github.com/fusion-energy/cad_to_dagmc).
 
-:point_right: Please see the [Online documentation](https://fusion-energy.github.io/paramak) for installation, usage and API documentation.
+:point_right: Please see the [online documentation](https://fusion-energy.github.io/paramak) for installation, usage and API documentation.
diff --git a/docs/usage.rst b/docs/usage.rst
index d1c72ea3..0fc36347 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -17,30 +17,7 @@ These examples show how to make various reactors with and without different comp
 .. toctree::
    :maxdepth: 3
 
+   usage_vis
    usage_tokamak
    usage_spherical_tokamak
 
-
-Visualization options
-=====================
-
-The reactors (cadquery.Assembly objects) and components (cadquery.Workplane objects) can be visualized in a number of ways.
-First option is to export them to STEP, BREP or STL file then open the file with your favorite CAD software for example `FreeCAD <https://www.freecad.org/>`_.
-See the `CadQuery documentation on saving <https://cadquery.readthedocs.io/en/latest/importexport.html#exporting-step>`_ for more information.
-Other options include the built in visualization tools in CadQuery
-
-.. code-block:: python
-
-    from cadquery.vis import show
-    show(result)  # where result is the returned reactor or component object
-
-or the jupyter_cadquery package which allows for interactive 3D visualization in a web browser.
-or the jupyter_cadquery package which allows for interactive 3D visualization in a web browser.
-
-.. code-block:: python
-
-    # pip install jupyter_cadquery
-    # might needed to downgrade pip with ... python -m pip  install pip==24.0
-    from jupyter_cadquery import show
-    view = show(result)  # where result is the returned reactor or component object
-    view.export_html("3d.html")
diff --git a/docs/usage_spherical_tokamak.rst b/docs/usage_spherical_tokamak.rst
index 90c68d2b..dde8cfd8 100644
--- a/docs/usage_spherical_tokamak.rst
+++ b/docs/usage_spherical_tokamak.rst
@@ -1,61 +1,7 @@
+Spherical Tokamak
+=================
 
-Spherical tokamak from plasma
------------------------------
-
-- The spherical_tokamak_from_plasma function provides a parametric tokamak shaped reactor.
-- This is characterized by a blanket that only goes around the outboard sides of the plasma.
-- This reactor requires few arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build.
-
-
-.. cadquery::
-    :gridsize: 0
-    :select: result
-    :color: #00cd00
-    :width: 100%
-    :height: 600px
-
-    import paramak
-    result = paramak.spherical_tokamak_from_plasma(
-        radial_build=[
-            (paramak.LayerType.GAP, 10),
-            (paramak.LayerType.SOLID, 60),
-            (paramak.LayerType.SOLID, 20),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.PLASMA, 300),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.SOLID, 20),
-            (paramak.LayerType.SOLID, 120),
-            (paramak.LayerType.SOLID, 10),
-        ],
-        elongation=2,
-        triangularity=0.55,
-        rotation_angle=90,
-    ).toCompound()
-
-
-.. code-block:: python
-
-    import paramak
-    result = paramak.spherical_tokamak_from_plasma(
-        radial_build=[
-            (paramak.LayerType.GAP, 10),
-            (paramak.LayerType.SOLID, 60),
-            (paramak.LayerType.SOLID, 20),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.PLASMA, 300),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.SOLID, 20),
-            (paramak.LayerType.SOLID, 120),
-            (paramak.LayerType.SOLID, 10),
-        ],
-        elongation=2,
-        triangularity=0.55,
-        rotation_angle=90,
-    )
-    result.save('reactor.step')
-
-
-
+- This is characterized by a blanket that only goes around the outboard sides of the plasma and the center column has no inboard breeding.
 
 Spherical tokamak
 -----------------
@@ -133,8 +79,64 @@ Spherical tokamak
 
     result.save(f"spherical_tokamak_minimal.step")
 
-Reactor with divertor(s)
-------------------------
+
+Spherical tokamak from plasma
+-----------------------------
+
+- The spherical_tokamak_from_plasma function provides a parametric tokamak shaped reactor.
+- This reactor requires minimal arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build.
+
+
+.. cadquery::
+    :gridsize: 0
+    :select: result
+    :color: #00cd00
+    :width: 100%
+    :height: 600px
+
+    import paramak
+    result = paramak.spherical_tokamak_from_plasma(
+        radial_build=[
+            (paramak.LayerType.GAP, 10),
+            (paramak.LayerType.SOLID, 60),
+            (paramak.LayerType.SOLID, 20),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.PLASMA, 300),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 20),
+            (paramak.LayerType.SOLID, 120),
+            (paramak.LayerType.SOLID, 10),
+        ],
+        elongation=2,
+        triangularity=0.55,
+        rotation_angle=90,
+    ).toCompound()
+
+
+.. code-block:: python
+
+    import paramak
+    result = paramak.spherical_tokamak_from_plasma(
+        radial_build=[
+            (paramak.LayerType.GAP, 10),
+            (paramak.LayerType.SOLID, 60),
+            (paramak.LayerType.SOLID, 20),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.PLASMA, 300),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 20),
+            (paramak.LayerType.SOLID, 120),
+            (paramak.LayerType.SOLID, 10),
+        ],
+        elongation=2,
+        triangularity=0.55,
+        rotation_angle=90,
+    )
+    result.save('reactor.step')
+
+
+Spherical tokamak with divertor(s)
+----------------------------------
 
 - ll reactors support adding additional radial builds for the lower_divertor and or the upper_divertor.
 - This example adds two divertors to a spherical_tokamak_from_plasma reactor but and other reactor would also work.
@@ -198,8 +200,8 @@ Reactor with divertor(s)
     )
     result.save('reactor.step')
 
-Reactor with poloidal field coils
----------------------------------
+Spherical tokamak with poloidal field coils
+-------------------------------------------
 
 - All reactors support adding a sequence of CadQuery shapes (e.g. workplanes) to the reactor using the extra_cut_shapes argument
 - This example adds PF coils to a spherical_tokamak_from_plasma reactor but and other reactor would also work.
@@ -298,8 +300,8 @@ Reactor with poloidal field coils
     result.save(f"spherical_tokamak_from_plasma_with_pf_magnets.step")
 
 
-Reactor with toroidal field coils
----------------------------------
+Spherical tokamak with toroidal field coils
+-------------------------------------------
 
 - In a similar way to adding poloidal field coils one can also add toroidal field coils by making use of the extra_cut_shapes argument.
 - All reactors support adding a sequence of CadQuery shapes (e.g. workplanes) to the reactor using the extra_cut_shapes argument
@@ -316,7 +318,7 @@ Reactor with toroidal field coils
 
     import paramak
 
-    tf = paramak.toroidal_field_coil_rectangle(
+    tf_style_1 = paramak.toroidal_field_coil_rectangle(
         horizontal_start_point = (10, 520),
         vertical_mid_point = (600, 0),
         thickness = 50,
@@ -325,7 +327,7 @@ Reactor with toroidal field coils
         azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
     )
 
-    result = paramak.spherical_tokamak_from_plasma(
+    result1 = paramak.spherical_tokamak_from_plasma(
         radial_build=[
             (paramak.LayerType.GAP, 70),
             (paramak.LayerType.SOLID, 10),
@@ -340,14 +342,47 @@ Reactor with toroidal field coils
         elongation=2.5,
         rotation_angle=180,
         triangularity=0.55,
-        extra_cut_shapes=[tf]
-    ).toCompound()
+        extra_cut_shapes=[tf_style_1]
+    ).toCompound().translate((700, 0, 0))
+
+    tf_style_2 = paramak.toroidal_field_coil_princeton_d(
+        r1=5,
+        r2=610,
+        azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
+        thickness = 50,
+        distance = 40
+    )
+
+    result2 = paramak.spherical_tokamak_from_plasma(
+        radial_build=[
+            (paramak.LayerType.GAP, 70),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.GAP, 50),
+            (paramak.LayerType.PLASMA, 300),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 60),
+            (paramak.LayerType.SOLID, 10),
+        ],
+        elongation=2.5,
+        rotation_angle=180,
+        triangularity=0.55,
+        extra_cut_shapes=[tf_style_2]
+    ).toCompound().translate((-700, 0, 0))
+
+    import cadquery as cq
+    result = cq.Assembly()
+    result.add(result1)
+    result.add(result2)
+    result = result.toCompound()
+
 
 .. code-block:: python
 
     import paramak
 
-    tf = paramak.toroidal_field_coil_rectangle(
+    tf_style_1 = paramak.toroidal_field_coil_rectangle(
         horizontal_start_point = (10, 520),
         vertical_mid_point = (600, 0),
         thickness = 50,
@@ -357,6 +392,34 @@ Reactor with toroidal field coils
     )
 
     result = paramak.spherical_tokamak_from_plasma(
+        radial_build=[
+            (paramak.LayerType.GAP, 70),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.GAP, 50),
+            (paramak.LayerType.PLASMA, 300),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 60),
+            (paramak.LayerType.SOLID, 10),
+        ],
+        elongation=2.5,
+        rotation_angle=180,
+        triangularity=0.55,
+        extra_cut_shapes=[tf_style_1]
+    )
+
+    result.save(f"spherical_tokamak_with_rectangular_tf.step")
+
+    tf_style_2 = paramak.toroidal_field_coil_princeton_d(
+        r1=5,
+        r2=610,
+        azimuthal_placement_angles = [120, 150, 180],
+        thickness = 50,
+        distance = 40
+    )
+
+    result2 = paramak.spherical_tokamak_from_plasma(
         radial_build=[
             (paramak.LayerType.GAP, 70),
             (paramak.LayerType.SOLID, 10),
@@ -374,7 +437,7 @@ Reactor with toroidal field coils
         extra_cut_shapes=[tf]
     )
 
-    result.save(f"spherical_tokamak_minimal.step")
+    result2.save(f"spherical_tokamak_with_princeton_tf.step")
 
 
 Spherical tokamak with negative triangularity
diff --git a/docs/usage_tokamak.rst b/docs/usage_tokamak.rst
index 6e05aa9b..c44987c2 100644
--- a/docs/usage_tokamak.rst
+++ b/docs/usage_tokamak.rst
@@ -1,10 +1,15 @@
 
-Tokamak from plasma
--------------------
+Tokamak
+=======
 
-- The tokamak_from_plasma function provides a parametric tokamak shaped reactor.
 - This is characterized by a continuous blanket that goes around the inboard and outboard sides of the plasma.
-- This reactor requires few arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build.
+
+Tokamak
+-------
+
+- The spherical_tokamak function provides a parametric tokamak shaped reactor.
+- This is characterized by a blanket that only goes around the outboard sides of the plasma.
+- This reactor allows for a separate vertical and radial build which allows different thickness layers in the blanket.
 
 .. cadquery::
     :gridsize: 0
@@ -14,7 +19,8 @@ Tokamak from plasma
     :height: 600px
 
     import paramak
-    result = paramak.tokamak_from_plasma(
+
+    result = paramak.tokamak(
         radial_build=[
             (paramak.LayerType.GAP, 10),
             (paramak.LayerType.SOLID, 30),
@@ -29,16 +35,26 @@ Tokamak from plasma
             (paramak.LayerType.SOLID, 120),
             (paramak.LayerType.SOLID, 10),
         ],
-        elongation=2,
+        vertical_build=[
+            (paramak.LayerType.SOLID, 15),
+            (paramak.LayerType.SOLID, 80),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.GAP, 50),
+            (paramak.LayerType.PLASMA, 700),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 40),
+            (paramak.LayerType.SOLID, 15),
+        ],
         triangularity=0.55,
-        rotation_angle=90,
+        rotation_angle=180,
     ).toCompound()
 
-
 .. code-block:: python
 
     import paramak
-    result = paramak.tokamak_from_plasma(
+
+    result = paramak.tokamak(
         radial_build=[
             (paramak.LayerType.GAP, 10),
             (paramak.LayerType.SOLID, 30),
@@ -53,18 +69,29 @@ Tokamak from plasma
             (paramak.LayerType.SOLID, 120),
             (paramak.LayerType.SOLID, 10),
         ],
-        elongation=2,
+        vertical_build=[
+            (paramak.LayerType.SOLID, 15),
+            (paramak.LayerType.SOLID, 80),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.GAP, 50),
+            (paramak.LayerType.PLASMA, 700),
+            (paramak.LayerType.GAP, 60),
+            (paramak.LayerType.SOLID, 10),
+            (paramak.LayerType.SOLID, 40),
+            (paramak.LayerType.SOLID, 15),
+        ],
         triangularity=0.55,
-        rotation_angle=90,
+        rotation_angle=180,
     )
 
+    result.save(f"tokamak_minimal.step")
 
-Tokamak
--------
 
-- The spherical_tokamak function provides a parametric tokamak shaped reactor.
-- This is characterized by a blanket that only goes around the outboard sides of the plasma.
-- This reactor allows for a separate vertical and radial build which allows different thickness layers in the blanket.
+Tokamak from plasma
+-------------------
+
+- The tokamak_from_plasma function provides a parametric tokamak shaped reactor.
+- This reactor requires few arguments to create as it keeps the vertical build of the blanket layers the same thickness as the radial build.
 
 .. cadquery::
     :gridsize: 0
@@ -74,8 +101,7 @@ Tokamak
     :height: 600px
 
     import paramak
-
-    result = paramak.tokamak(
+    result = paramak.tokamak_from_plasma(
         radial_build=[
             (paramak.LayerType.GAP, 10),
             (paramak.LayerType.SOLID, 30),
@@ -90,26 +116,16 @@ Tokamak
             (paramak.LayerType.SOLID, 120),
             (paramak.LayerType.SOLID, 10),
         ],
-        vertical_build=[
-            (paramak.LayerType.SOLID, 15),
-            (paramak.LayerType.SOLID, 80),
-            (paramak.LayerType.SOLID, 10),
-            (paramak.LayerType.GAP, 50),
-            (paramak.LayerType.PLASMA, 700),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.SOLID, 10),
-            (paramak.LayerType.SOLID, 40),
-            (paramak.LayerType.SOLID, 15),
-        ],
+        elongation=2,
         triangularity=0.55,
-        rotation_angle=180,
+        rotation_angle=90,
     ).toCompound()
 
+
 .. code-block:: python
 
     import paramak
-
-    result = paramak.tokamak(
+    result = paramak.tokamak_from_plasma(
         radial_build=[
             (paramak.LayerType.GAP, 10),
             (paramak.LayerType.SOLID, 30),
@@ -124,23 +140,22 @@ Tokamak
             (paramak.LayerType.SOLID, 120),
             (paramak.LayerType.SOLID, 10),
         ],
-        vertical_build=[
-            (paramak.LayerType.SOLID, 15),
-            (paramak.LayerType.SOLID, 80),
-            (paramak.LayerType.SOLID, 10),
-            (paramak.LayerType.GAP, 50),
-            (paramak.LayerType.PLASMA, 700),
-            (paramak.LayerType.GAP, 60),
-            (paramak.LayerType.SOLID, 10),
-            (paramak.LayerType.SOLID, 40),
-            (paramak.LayerType.SOLID, 15),
-        ],
+        elongation=2,
         triangularity=0.55,
-        rotation_angle=180,
+        rotation_angle=90,
     )
 
-    result.save(f"tokamak_minimal.step")
+Tokamak with divertor(s)
+------------------------
+TODO
+
+Tokamak with poloidal field coils
+---------------------------------
+TODO
 
+Tokamak with toroidal field coils
+---------------------------------
+TODO
 
 
 Tokamak with negative triangularity
diff --git a/docs/usage_vis.rst b/docs/usage_vis.rst
new file mode 100644
index 00000000..c80d7373
--- /dev/null
+++ b/docs/usage_vis.rst
@@ -0,0 +1,24 @@
+
+Visualization options
+=====================
+
+The reactors (cadquery.Assembly objects) and components (cadquery.Workplane objects) can be visualized in a number of ways.
+First option is to export them to STEP, BREP or STL file then open the file with your favorite CAD software for example `FreeCAD <https://www.freecad.org/>`_.
+See the `CadQuery documentation on saving <https://cadquery.readthedocs.io/en/latest/importexport.html#exporting-step>`_ for more information.
+Other options include the built in visualization tools in CadQuery
+
+.. code-block:: python
+
+    from cadquery.vis import show
+    show(result)  # where result is the returned reactor or component object
+
+Jupyter_cadquery package can be installed (pip install jupyter-cadquery) which allows for interactive 3D visualization in a web browser.
+This package can also be used to export portable html files with the 3D visualization.
+
+.. code-block:: python
+
+    # pip install jupyter_cadquery
+    # might needed to downgrade pip with ... python -m pip  install pip==24.0
+    from jupyter_cadquery import show
+    view = show(result)  # where result is the returned reactor or component object
+    view.export_html("3d.html")
diff --git a/examples/spherical_tokamak_from_plasma_with_tf_magnets.py b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py
index d458fa83..fdc5e4fa 100644
--- a/examples/spherical_tokamak_from_plasma_with_tf_magnets.py
+++ b/examples/spherical_tokamak_from_plasma_with_tf_magnets.py
@@ -4,13 +4,21 @@
 
 import paramak
 
-tf = paramak.toroidal_field_coil_rectangle(
+tf_style_1 = paramak.toroidal_field_coil_rectangle(
     horizontal_start_point = (10, 520),
     vertical_mid_point = (600, 0),
     thickness = 50,
     distance = 40,
     with_inner_leg = True,
-    azimuthal_placement_angles = [0, 30, 60, 90, 120, 150, 180],
+    azimuthal_placement_angles = [0, 30, 60, 90],
+)
+
+tf_style_2 = paramak.toroidal_field_coil_princeton_d(
+    r1=5,
+    r2=610,
+    azimuthal_placement_angles = [120, 150, 180],
+    thickness = 50,
+    distance = 40
 )
 
 result = paramak.spherical_tokamak_from_plasma(
@@ -28,7 +36,7 @@
     elongation=2.5,
     rotation_angle=180,
     triangularity=0.55,
-    extra_cut_shapes=[tf]
+    extra_cut_shapes=[tf_style_1, tf_style_2]
 )
 
 result.save(f"spherical_tokamak_minimal.step")
diff --git a/src/paramak/__init__.py b/src/paramak/__init__.py
index f793fab8..376795d8 100644
--- a/src/paramak/__init__.py
+++ b/src/paramak/__init__.py
@@ -14,6 +14,7 @@
 from .workplanes.poloidal_field_coil_case import poloidal_field_coil_case
 from .workplanes.toroidal_field_coil_rectangle import toroidal_field_coil_rectangle
 from .workplanes.u_shaped_dome import u_shaped_dome
+from .workplanes.toroidal_field_coil_princeton_d import toroidal_field_coil_princeton_d
 
 from .utils import LayerType
 
diff --git a/src/paramak/workplanes/toroidal_field_coil_princeton_d.py b/src/paramak/workplanes/toroidal_field_coil_princeton_d.py
new file mode 100644
index 00000000..28870aa5
--- /dev/null
+++ b/src/paramak/workplanes/toroidal_field_coil_princeton_d.py
@@ -0,0 +1,206 @@
+import typing
+import numpy as np
+from ..utils import create_wire_workplane_from_points, rotate_solid
+from scipy import integrate
+from scipy.optimize import minimize
+from typing import List, Tuple
+
+def _compute_inner_points(R1, R2):
+    """Computes the inner curve points
+
+    Args:
+        R1 (float): smallest radius (cm)
+        R2 (float): largest radius (cm)
+
+    Returns:
+        (list, list, list): R, Z and derivative lists for outer curve
+        points
+    """
+
+    def error(z_0, R0, R2):
+        segment = get_segment(R0, R2, z_0)
+        return abs(segment[1][-1])
+
+    def get_segment(a, b, z_0,num=5):
+        a_R = np.linspace(a, b, num=num, endpoint=True)
+        asol = integrate.odeint(solvr, [z_0[0], 0], a_R)
+        return a_R, asol[:, 0], asol[:, 1]
+
+    def solvr(Y, R):
+        return [Y[1], -1 / (k * R) * (1 + Y[1] ** 2) ** (3 / 2)]
+
+    R0 = (R1 * R2) ** 0.5
+    k = 0.5 * np.log(R2 / R1)
+
+    # computing of z_0
+    # z_0 is computed by ensuring outer segment end is zero
+    z_0 = 10  # initial guess for z_0
+    res = minimize(error, z_0, args=(R0, R2))
+    z_0 = res.x
+
+    # compute inner and outer segments
+    segment1 = get_segment(R0, R1, z_0)
+    segment2 = get_segment(R0, R2, z_0)
+
+    r_values = np.concatenate(
+        [
+            np.flip(segment1[0]),
+            segment2[0][1:],
+            np.flip(segment2[0])[1:],
+            segment1[0][1:],
+        ]
+    )
+    z_values = np.concatenate(
+        [
+            np.flip(segment1[1]),
+            segment2[1][1:],
+            -np.flip(segment2[1])[1:],
+            -segment1[1][1:],
+        ]
+    )
+    return r_values, z_values
+
+
+def add_thickness(x: List[float], y: List[float], thickness: float, dy_dx: List[float] = None) -> Tuple[list, list]:
+    """Computes outer curve points based on thickness
+
+    Args:
+        x (list): list of floats containing x values
+        y (list): list of floats containing y values
+        thickness (float): thickness of the magnet
+        dy_dx (list): list of floats containing the first order
+            derivatives
+
+    Returns:
+        R and Z lists for outer curve points
+    """
+
+    if dy_dx is None:
+        dy_dx = np.diff(y) / np.diff(x)
+
+    x_outer, y_outer = [], []
+    for i in range(len(dy_dx)):
+        if dy_dx[i] == float("-inf"):
+            nx, ny = -1, 0
+        elif dy_dx[i] == float("inf"):
+            nx, ny = 1, 0
+        else:
+            nx = -dy_dx[i]
+            ny = 1
+        if i != len(dy_dx) - 1:
+            if x[i] < x[i + 1]:
+                convex = False
+            else:
+                convex = True
+
+        if convex:
+            nx *= -1
+            ny *= -1
+        # normalize normal vector
+        normal_vector_norm = (nx**2 + ny**2) ** 0.5
+        nx /= normal_vector_norm
+        ny /= normal_vector_norm
+        # calculate outer points
+        val_x_outer = x[i] + thickness * nx
+        val_y_outer = y[i] + thickness * ny
+        x_outer.append(val_x_outer)
+        y_outer.append(val_y_outer)
+
+    return x_outer, y_outer
+
+
+def find_points(R1, R2, thickness,vertical_displacement):
+    """Finds the XZ points joined by connections that describe the 2D
+    profile of the toroidal field coil shape."""
+    # compute inner points
+    r_inner, z_inner = _compute_inner_points(R1 + thickness, R2)
+
+    # compute outer points
+    dz_dr = np.diff(z_inner) / np.diff(r_inner)
+    dz_dr[0] = float("-inf")
+    dz_dr = np.append(dz_dr, float("inf"))
+    r_outer, z_outer = add_thickness(r_inner, z_inner, thickness, dy_dx=dz_dr)
+    r_outer, z_outer = np.flip(r_outer), np.flip(z_outer)
+
+    # add vertical displacement
+    z_outer += vertical_displacement
+    z_inner += vertical_displacement
+
+    # extract helping points for inner leg
+    inner_leg_connection_points = [
+        (r_inner[0], z_inner[0]),
+        (r_inner[-1], z_inner[-1]),
+        (r_outer[0], z_outer[0]),
+        (r_outer[-1], z_outer[-1]),
+    ]
+
+    # add connections
+    inner_points = [[r, z, "spline"] for r, z in zip(r_inner, z_inner)]
+    outer_points = [[r, z, "spline"] for r, z in zip(r_outer, z_outer)]
+
+    inner_points[-1][2] = "straight"
+    outer_points[-1][2] = "straight"
+
+    points = inner_points + outer_points
+    outer_points = np.vstack((r_outer, z_outer)).T
+    inner_points = np.vstack((r_inner, z_inner)).T
+
+    return points, inner_leg_connection_points, inner_points, outer_points
+
+def toroidal_field_coil_princeton_d(
+    r1: float = 100,
+    r2: float = 300,
+    thickness: float = 30,
+    distance: float = 20,
+    name: str = "toroidal_field_coil",
+    with_inner_leg: bool = True,
+    azimuthal_placement_angles: typing.Sequence[float] = [0],
+    vertical_displacement: float = 0.0,
+    color: typing.Tuple[float, float, float, typing.Optional[float]] = (0.0, 0.0, 1.0),
+    plane: str = "XZ",
+    origin: typing.Tuple[float, float, float] = (0.0, 0.0, 0.0),
+    obj=None,
+):
+    """
+    Creates a toroidal field coil with a Princeton-D shape.
+
+    Args:
+        r1 (float, optional): Inner radius of the coil. Defaults to 100.
+        r2 (float, optional): Outer radius of the coil. Defaults to 300.
+        thickness (float, optional): Thickness of the coil. Defaults to 30.
+        distance (float, optional): Distance to extrude the coil. Defaults to 20.
+        name (str, optional): Name of the coil. Defaults to "toroidal_field_coil".
+        with_inner_leg (bool, optional): Whether to include the inner leg of the coil. Defaults to True.
+        azimuthal_placement_angles (typing.Sequence[float], optional): Angles for azimuthal placement. Defaults to [0].
+        vertical_displacement (float, optional): Vertical displacement of the coil. Defaults to 0.0.
+        color (typing.Tuple[float, float, float, typing.Optional[float]], optional): Color of the coil. Defaults to (0.0, 0.0, 1.0).
+        plane (str, optional): Plane in which to create the coil. Defaults to "XZ".
+        origin (typing.Tuple[float, float, float], optional): Origin point for the coil. Defaults to (0.0, 0.0, 0.0).
+        obj (optional): Existing object to modify. Defaults to None.
+
+    Returns:
+        solid: The created toroidal field coil solid.
+    """
+    points, inner_leg_connection_points, inner_points, outer_points = find_points(r1, r2, thickness, vertical_displacement)
+    # need to get square end, it appears to miss the last point in the solid, TODO fix so this append is not needed
+    points.append(points[-1])
+    wire = create_wire_workplane_from_points(points=points, plane=plane, origin=origin, obj=obj)
+    solid = wire.extrude(until=distance / 2, both=True)
+    solid = rotate_solid(angles=azimuthal_placement_angles, solid=solid)
+
+    if with_inner_leg:
+        inner_leg_connection_points=[(x,z,'straight') for x,z in inner_leg_connection_points]
+        # need to get square end, it appears to miss the last point in the solid, TODO fix so this append is not needed
+        inner_leg_connection_points.append(inner_leg_connection_points[-1])
+        for i in inner_leg_connection_points:
+            print('inner_leg_connection_points' , i)
+        inner_wire = create_wire_workplane_from_points(
+            points=inner_leg_connection_points, plane=plane, origin=origin, obj=obj
+        )
+        inner_solid = inner_wire.extrude(until=distance / 2, both=True)
+        inner_solid = rotate_solid(angles=azimuthal_placement_angles, solid=inner_solid)
+        solid = solid.union(inner_solid)
+
+    solid.name = name
+    solid.color = color
+    return solid
\ No newline at end of file
diff --git a/src/paramak/workplanes/toroidal_field_coil_rectangle.py b/src/paramak/workplanes/toroidal_field_coil_rectangle.py
index 61c46977..f8ad3c74 100644
--- a/src/paramak/workplanes/toroidal_field_coil_rectangle.py
+++ b/src/paramak/workplanes/toroidal_field_coil_rectangle.py
@@ -10,7 +10,7 @@ def toroidal_field_coil_rectangle(
     distance: float = 20,
     name: str = "toroidal_field_coil",
     with_inner_leg: bool = True,
-    azimuthal_placement_angles: typing.Sequence[float] = [0, 90, 180],
+    azimuthal_placement_angles: typing.Sequence[float] = [0],
     vertical_displacement: float = 0.0,
     color: typing.Tuple[float, float, float, typing.Optional[float]] = (0.0, 0.0, 1.0),
     plane: str = "XZ",
diff --git a/tests/test_workplanes/test_toroidal_field_coil_princeton_d.py b/tests/test_workplanes/test_toroidal_field_coil_princeton_d.py
new file mode 100644
index 00000000..d1f6cf0d
--- /dev/null
+++ b/tests/test_workplanes/test_toroidal_field_coil_princeton_d.py
@@ -0,0 +1,8 @@
+
+import paramak
+
+def test_creation_of_inner_leg():
+    solid_without_inner_leg = paramak.toroidal_field_coil_princeton_d(with_inner_leg=False)
+    solid_with_inner_leg = paramak.toroidal_field_coil_princeton_d(with_inner_leg=True)
+
+    assert solid_without_inner_leg.val().Volume() < solid_with_inner_leg.val().Volume()
\ No newline at end of file