diff --git a/README.md b/README.md
index e524c3f..3432c3f 100644
--- a/README.md
+++ b/README.md
@@ -12,6 +12,7 @@
     - standard [ISCWSA] method
     - new mesh based method using the [Flexible Collision Library]
 ## New Features!
+  - **Well Path Creation:** the addition of the `connector` module enables drilling well paths to be created simply by providing start and end locations (with some vector data like inclination and azimuth). No need to indicate *how* to connect the points, the module will figure that out itself.
   - **Fast visualization of well trajectory meshes:** addition of the `visual` module for quick and simple viewing and QAQC of well meshes.
   - **Mesh Based Collision Detection:** the current method for determining the Separation Factor between wells is constrained by the frequency and location of survey stations or necessitates interpolation of survey stations in order to determine if Anti-Collision Rules have been violated. Meshing the well bore interpolates between survey stations and as such is a more reliable method for identifying potential well bore collisions, especially wth more sparse data sets.
   - More coming soon!
@@ -36,7 +37,7 @@ pip install -e .
 ```
 Make sure you include that `.` in the final line (it's not a typo) as this ensures that any changes to your development version are immediately implemented on save.
 ## Quick Start
-Here's an example using `welleng` to construct a couple of simple well trajectories with `numpy`, creating survey listings for the wells with well bore uncertainty data, using these surveys to create well bore meshes and finally printing the results and plotting the meshes with the closest lines and SF data.
+Here's an example using `welleng` to construct a couple of simple well trajectories with the `connector` module, creating survey listings for the wells with well bore uncertainty data, using these surveys to create well bore meshes and finally printing the results and plotting the meshes with the closest lines and SF data.
 
 ```
 import welleng as we
@@ -45,29 +46,37 @@ from tabulate import tabulate
 
 # construct simple well paths
 print("Constructing wells...")
-md = np.linspace(0,3000,100) # 30 meter intervals to 3000 mTD
-inc = np.concatenate((
-    np.zeros(30), # vertical section
-    np.linspace(0,90,60), # build section to 60 degrees
-    np.full(10,90) # hold section at 60 degrees
-))
-azi1 = np.full(100,60) # constant azimuth at 60 degrees
-azi2 = np.full(100,225) # constant azimuth at 225 degrees
-
-# make a survey object and calculate the uncertainty covariances
+connector_reference = we.connector.Connector(
+    pos1=[0,0,0],
+    inc1=0,
+    azi1=0,
+    pos2=[-100,0,2000.],
+    inc2=90,
+    azi2=60,
+).survey(step=50)
+
+connector_offset = we.connector.Connector(
+    pos1=[0,0,0],
+    inc1=0,
+    azi1=225,
+    pos2=[-280,-600,2000],
+    inc2=90,
+    azi2=270,
+).survey(step=50)
+
+# make a survey objects and calculate the uncertainty covariances
 print("Making surveys...")
 survey_reference = we.survey.Survey(
-    md,
-    inc,
-    azi1,
+    md=connector_reference.md,
+    inc=connector_reference.inc_deg,
+    azi=connector_reference.azi_deg,
     error_model='ISCWSA_MWD'
 )
 
-# make another survey with offset surface location and along another azimuth
 survey_offset = we.survey.Survey(
-    md,
-    inc,
-    azi2,
+    md=connector_offset.md,
+    inc=connector_offset.inc_deg,
+    azi=connector_offset.azi_deg,
     start_nev=[100,200,0],
     error_model='ISCWSA_MWD'
 )
@@ -117,21 +126,21 @@ we.visual.plot(
 
 print("Done!")
 ```
-This results in a quick, interactive visualization of the well meshes that's great for QAQC.
+This results in a quick, interactive visualization of the well meshes that's great for QAQC. What's interesting about these results is that the ISCWSA method does not explicitly detect a collision in this scenario wheras the mesh method does.
 
-![image](https://user-images.githubusercontent.com/41046859/100718537-c3b12700-33bb-11eb-856e-cf1bd77d3cbf.png)
+![image](https://user-images.githubusercontent.com/41046859/102106351-c0dd1a00-3e30-11eb-82f0-a0454dfce1c6.png)
 
 For more examples, check out the [examples].
 
 ## Todos
- - Add a Target class to see what you're aiming for!
- - Export to Landmark's .wbp format so survey listings can be modified in COMPASS
+ - Add a Target class to see what you're aiming for - **in progress**
+ - Export to Landmark's .wbp format so survey listings can be modified in COMPASS - **in progress**
  - Documentation
  - Generate a scene of offset wells to enable fast screening of collision risks (e.g. hundreds of wells in seconds)
- - Well trajectory planning - construct your own trajectories using a range of methods (and of course, including some novel ones)
+ - Well trajectory planning - construct your own trajectories using a range of methods (and of course, including some novel ones) **- DONE!**
  - More error models
  - WebApp for those that just want answers
- - Viewer - a 3D viewer to quickly visualize the data and calculated results - **DONE!**
+ - Viewer - a 3D viewer to quickly visualize the data and calculated results **- DONE!**
 
 It's possible to generate data for visualizing well trajectories with [welleng], as can be seen with the rendered scenes below.
 ![image](https://user-images.githubusercontent.com/41046859/97724026-b78c2e00-1acc-11eb-845d-1220219843a5.png)
diff --git a/examples/connect_two_random_points.py b/examples/connect_two_random_points.py
new file mode 100644
index 0000000..12be7b9
--- /dev/null
+++ b/examples/connect_two_random_points.py
@@ -0,0 +1,152 @@
+import welleng as we
+import numpy as np
+from vedo import Arrows, Lines
+import random
+
+# Some code for testing the connector module.
+
+# Generate some random pairs of points
+pos1 = [0,0,0]
+md1 = 0
+
+pos2 = np.random.random(3) * 1000
+
+vec1 = np.random.random(3)
+vec1 /= np.linalg.norm(vec1)
+inc1, azi1 = np.degrees(we.utils.get_angles(vec1, nev=True)).reshape(2).T
+
+vec2 = np.random.random(3)
+vec2 /= np.linalg.norm(vec2)
+
+inc2, azi2 = np.degrees(we.utils.get_angles(vec2, nev=True)).reshape(2).T
+
+md2 = 100 + random.random() * 1000
+
+# Define some random input permutations
+
+number = 7
+rand = random.random()
+
+# 1: test only md2 (hold)
+if rand < 1 * 1 / number:
+    option = 1
+    expected_method = 'hold'
+    vec2, pos2, inc2, azi2 = None, None, None, None
+
+# 2: test md2 and an inc2
+elif rand < 1 * 2 / number:
+    option = 2
+    expected_method = 'min_curve'
+    pos2, vec2, azi2 = None, None, None
+
+# 3: test md2 and azi2
+elif rand < 1 * 3 / number:
+    option = 3
+    expected_method = 'min_curve'
+    pos2, vec2, inc2 = None, None, None
+
+# 4: test md2, inc2 and azi2
+elif rand < 1 * 4 / number:
+    option = 4
+    expected_method = 'min_curve'
+    pos2, vec2 = None, None
+
+# 5 test pos2
+elif rand < 1 * 5 / number:
+    option = 5
+    expected_method = 'min_dist_to_target'
+    vec2, inc2, azi2, md2 = None, None, None, None
+
+# 6 test pos2 vec2
+elif rand < 1 * 6 / number:
+    option = 6
+    expected_method = 'curve_hold_curve'
+    md2, inc2, azi2, = None, None, None
+
+# 7 test pos2, inc2 and azi2
+else:
+    option = 7
+    expected_method = 'curve_hold_curve'
+    md2, vec2 = None, None
+
+# Print the input parameters
+
+print(
+    f"Option: {option}\tExpected Method: {expected_method}\n"
+    f"md1: {md1}\tpos1: {pos1}\tvec1: {vec1}\tinc1: {inc1}\tazi1: {azi1}\n"
+    f"md2: {md2}\tpos2: {pos2}\tvec2: {vec2}\tinc2: {inc2}\tazi2: {azi2}\n"
+)
+
+# Initialize a connector object and connect the inputs
+section = we.connector.Connector(
+    pos1=[0.,0.,0],
+    vec1=vec1,
+    md2=md2,
+    pos2=pos2,
+    vec2=vec2,
+    inc2=inc2,
+    azi2=azi2,
+    degrees=True
+)
+
+# Print some pertinent calculation data
+
+print(
+    f"Method: {section.method}\n",
+    f"radius_design: {section.radius_design}\t",
+    f"radius_critical: {section.radius_critical}\n"
+    f"radius_design2: {section.radius_design2}\t",
+    f"radius_critical2: {section.radius_critical2}\n"
+    f"iterations: {section.iterations}"
+)
+
+# Create a survey object of the section with interpolated points and coords
+survey = section.survey(radius=5, step=30)
+
+# As a QAQC step, check that the wellpath hits the defined turn points
+start_points = np.array([section.pos1])
+end_points = np.array([section.pos_target])
+if section.pos2 is not None:
+    start_points = np.concatenate((start_points, [section.pos2]))
+    end_points = np.concatenate(([section.pos2], [section.pos_target]))
+if section.pos3 is not None:
+    start_points = np.concatenate((start_points, [section.pos3]))
+    end_points = np.concatenate(([section.pos2], [section.pos3], [section.pos_target]))
+lines = Lines(
+    startPoints=start_points,
+    endPoints=end_points,
+    c='green',
+    lw=5
+)
+
+# Add some arrows to represent the vectors at the start and end positions
+scalar=150
+arrows = Arrows(
+    startPoints=np.array([
+        section.pos1,
+        section.pos_target
+    ]),
+    endPoints=np.array([
+        section.pos1 + scalar * section.vec1,
+        section.pos_target + scalar * section.vec_target
+    ]),
+    s=0.5,
+    res=24
+)
+
+# generate a mesh of the generated section from the survey
+# use the 'circle' method to construct a cylinder with constant radius
+mesh = we.mesh.WellMesh(
+    survey=survey,
+    method='circle',
+    n_verts=24,
+)
+
+# plot the results
+we.visual.plot(
+    [mesh.mesh],
+    lines=lines,
+    arrows=arrows,
+)
+
+print("Done")
\ No newline at end of file
diff --git a/examples/simple_example.py b/examples/simple_example.py
index 0956a90..4748d38 100644
--- a/examples/simple_example.py
+++ b/examples/simple_example.py
@@ -4,29 +4,37 @@
 
 # construct simple well paths
 print("Constructing wells...")
-md = np.linspace(0,3000,100) # 30 meter intervals to 3000 mTD
-inc = np.concatenate((
-    np.zeros(30), # vertical section
-    np.linspace(0,90,60), # build section to 60 degrees
-    np.full(10,90) # hold section at 60 degrees
-))
-azi1 = np.full(100,60) # constant azimuth at 60 degrees
-azi2 = np.full(100,225) # constant azimuth at 225 degrees
+connector_reference = we.connector.Connector(
+    pos1=[0,0,0],
+    inc1=0,
+    azi1=0,
+    pos2=[-100,0,2000.],
+    inc2=90,
+    azi2=60,
+).survey(step=50)
 
-# make a survey object and calculate the uncertainty covariances
+connector_offset = we.connector.Connector(
+    pos1=[0,0,0],
+    inc1=0,
+    azi1=225,
+    pos2=[-280,-600,2000],
+    inc2=90,
+    azi2=270,
+).survey(step=50)
+
+# make a survey objects and calculate the uncertainty covariances
 print("Making surveys...")
 survey_reference = we.survey.Survey(
-    md,
-    inc,
-    azi1,
+    md=connector_reference.md,
+    inc=connector_reference.inc_deg,
+    azi=connector_reference.azi_deg,
     error_model='ISCWSA_MWD'
 )
 
-# make another survey with offset surface location and along another azimuth
 survey_offset = we.survey.Survey(
-    md,
-    inc,
-    azi2,
+    md=connector_offset.md,
+    inc=connector_offset.inc_deg,
+    azi=connector_offset.azi_deg,
     start_nev=[100,200,0],
     error_model='ISCWSA_MWD'
 )
diff --git a/requirements.txt b/requirements.txt
index abad2f5..c5c7f33 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,20 @@
-numpy
-scipy
-trimesh
-vedo
\ No newline at end of file
+cycler==0.10.0
+Cython==0.29.21
+decorator==4.4.2
+et-xmlfile==1.0.1
+jdcal==1.4.1
+kiwisolver==1.3.1
+matplotlib==3.3.3
+networkx==2.5
+numpy==1.19.4
+openpyxl==3.0.5
+Pillow==8.0.1
+pyparsing==2.4.7
+python-dateutil==2.8.1
+python-fcl==0.0.12
+scipy==1.5.4
+six==1.15.0
+tabulate==0.8.7
+trimesh==3.8.18
+vedo==2020.4.2
+vtk==9.0.1
diff --git a/welleng/__init__.py b/welleng/__init__.py
index 5fcc8ef..6c1308a 100644
--- a/welleng/__init__.py
+++ b/welleng/__init__.py
@@ -6,4 +6,7 @@
 import welleng.mesh
 import welleng.visual
 import welleng.version
-import welleng.errors.iscwsa_mwd
\ No newline at end of file
+import welleng.errors.iscwsa_mwd
+import welleng.exchange.wbp
+import welleng.target
+import welleng.connector
\ No newline at end of file
diff --git a/welleng/__pycache__/__init__.cpython-37.pyc b/welleng/__pycache__/__init__.cpython-37.pyc
index 6d3b6a5..6166454 100644
Binary files a/welleng/__pycache__/__init__.cpython-37.pyc and b/welleng/__pycache__/__init__.cpython-37.pyc differ
diff --git a/welleng/__pycache__/__init__.cpython-38.pyc b/welleng/__pycache__/__init__.cpython-38.pyc
index 3632c40..ec297e4 100644
Binary files a/welleng/__pycache__/__init__.cpython-38.pyc and b/welleng/__pycache__/__init__.cpython-38.pyc differ
diff --git a/welleng/__pycache__/clearance.cpython-38.pyc b/welleng/__pycache__/clearance.cpython-38.pyc
index 82aaad8..563c6cb 100644
Binary files a/welleng/__pycache__/clearance.cpython-38.pyc and b/welleng/__pycache__/clearance.cpython-38.pyc differ
diff --git a/welleng/__pycache__/connector.cpython-37.pyc b/welleng/__pycache__/connector.cpython-37.pyc
new file mode 100644
index 0000000..3f8bd4b
Binary files /dev/null and b/welleng/__pycache__/connector.cpython-37.pyc differ
diff --git a/welleng/__pycache__/connector.cpython-38.pyc b/welleng/__pycache__/connector.cpython-38.pyc
new file mode 100644
index 0000000..3c8c0f0
Binary files /dev/null and b/welleng/__pycache__/connector.cpython-38.pyc differ
diff --git a/welleng/__pycache__/error.cpython-38.pyc b/welleng/__pycache__/error.cpython-38.pyc
index 99083b9..861f651 100644
Binary files a/welleng/__pycache__/error.cpython-38.pyc and b/welleng/__pycache__/error.cpython-38.pyc differ
diff --git a/welleng/__pycache__/mesh.cpython-37.pyc b/welleng/__pycache__/mesh.cpython-37.pyc
index 164b40c..4dede42 100644
Binary files a/welleng/__pycache__/mesh.cpython-37.pyc and b/welleng/__pycache__/mesh.cpython-37.pyc differ
diff --git a/welleng/__pycache__/mesh.cpython-38.pyc b/welleng/__pycache__/mesh.cpython-38.pyc
index bbc8003..76d54d2 100644
Binary files a/welleng/__pycache__/mesh.cpython-38.pyc and b/welleng/__pycache__/mesh.cpython-38.pyc differ
diff --git a/welleng/__pycache__/survey.cpython-37.pyc b/welleng/__pycache__/survey.cpython-37.pyc
index a7610c0..17d60c5 100644
Binary files a/welleng/__pycache__/survey.cpython-37.pyc and b/welleng/__pycache__/survey.cpython-37.pyc differ
diff --git a/welleng/__pycache__/survey.cpython-38.pyc b/welleng/__pycache__/survey.cpython-38.pyc
index 13b1c26..1c218cf 100644
Binary files a/welleng/__pycache__/survey.cpython-38.pyc and b/welleng/__pycache__/survey.cpython-38.pyc differ
diff --git a/welleng/__pycache__/target.cpython-37.pyc b/welleng/__pycache__/target.cpython-37.pyc
new file mode 100644
index 0000000..84e0939
Binary files /dev/null and b/welleng/__pycache__/target.cpython-37.pyc differ
diff --git a/welleng/__pycache__/target.cpython-38.pyc b/welleng/__pycache__/target.cpython-38.pyc
new file mode 100644
index 0000000..4505827
Binary files /dev/null and b/welleng/__pycache__/target.cpython-38.pyc differ
diff --git a/welleng/__pycache__/utils.cpython-37.pyc b/welleng/__pycache__/utils.cpython-37.pyc
index fa0d8a0..6f0cbfc 100644
Binary files a/welleng/__pycache__/utils.cpython-37.pyc and b/welleng/__pycache__/utils.cpython-37.pyc differ
diff --git a/welleng/__pycache__/utils.cpython-38.pyc b/welleng/__pycache__/utils.cpython-38.pyc
index e90c3f1..1e2ca19 100644
Binary files a/welleng/__pycache__/utils.cpython-38.pyc and b/welleng/__pycache__/utils.cpython-38.pyc differ
diff --git a/welleng/__pycache__/version.cpython-37.pyc b/welleng/__pycache__/version.cpython-37.pyc
index f9477cb..3cec3b0 100644
Binary files a/welleng/__pycache__/version.cpython-37.pyc and b/welleng/__pycache__/version.cpython-37.pyc differ
diff --git a/welleng/__pycache__/version.cpython-38.pyc b/welleng/__pycache__/version.cpython-38.pyc
new file mode 100644
index 0000000..d9b6cf5
Binary files /dev/null and b/welleng/__pycache__/version.cpython-38.pyc differ
diff --git a/welleng/__pycache__/visual.cpython-37.pyc b/welleng/__pycache__/visual.cpython-37.pyc
index b1c7f51..2bebda7 100644
Binary files a/welleng/__pycache__/visual.cpython-37.pyc and b/welleng/__pycache__/visual.cpython-37.pyc differ
diff --git a/welleng/__pycache__/visual.cpython-38.pyc b/welleng/__pycache__/visual.cpython-38.pyc
new file mode 100644
index 0000000..d491faa
Binary files /dev/null and b/welleng/__pycache__/visual.cpython-38.pyc differ
diff --git a/welleng/connector.py b/welleng/connector.py
new file mode 100644
index 0000000..1b4fd44
--- /dev/null
+++ b/welleng/connector.py
@@ -0,0 +1,1150 @@
+import numpy as np
+import math as m
+from .utils import get_vec, get_angles
+from .survey import Survey
+
+
+class Connector:
+    def __init__(
+        self,
+        pos1=[0,0,0],
+        vec1=None,
+        inc1=None,
+        azi1=None,
+        md1=0,
+        dls_design=3.0,
+        dls_design2=None,
+        md2=None,
+        pos2=None,
+        vec2=None,
+        inc2=None,
+        azi2=None,
+        degrees=True,
+        unit='meters',
+        min_error=1e-5,
+        delta_radius=10,
+        min_tangent=10,
+        max_iterations=1000,
+    ):
+        
+        """
+        A class to provide a fast, efficient method for determining well bore 
+        section geometry using the appropriate method based upon the provided
+        parameters, with the intent of assisting machine learning fitness
+        fitness functions. Interpolation between the returned control points
+        can be performed postumously - attempts are made to keep processing to
+        a minimum in the event that the Connector is being used for machine
+        learning.
+
+        Only specific combinations of input data are permitted, e.g. if you
+        input both a start vector AND a start inc and azi you'll get an error
+        (which one is correct after all?). Think about what you're trying to
+        achieve and provide the data that will facilitate that outcome.
+
+        Parameters
+        ----------
+            pos1: (3) list or array of floats (default: [0,0,0])
+                Start position in NEV coordinates.
+            vec1: (3) list or array of floats or None (default: None)
+                Start position unit vector in NEV coordinates.
+            inc1: float or None (default: None)
+                Start position inclination.
+            azi2: float or None (default: None)
+                Start position azimuth.
+            md1: float or None (default: None)
+                Start position measured depth.
+            dls_design: float (default: 3.0)
+                The desired Dog Leg Severity (DLS) for the (first) curved
+                section in degrees per 30 meters or 100 feet.
+            dls_design2: float or None (default: None)
+                The desired DLS for the second curve section in degrees per
+                30 meters or 100 feet. If set to None then `dls_design` will
+                be the default value.
+            md2: float or None (default: None)
+                The measured depth of the target position.
+            pos2: (3) list or array of floats or None (default: None)
+                The position of the target in NEV coordinates.
+            vec2: (3) list or array of floats or None (default: None)
+                The target unit vector in NEV coordinates.
+            inc1: float or None (default: None)
+                The inclination at the target position.
+            azi2: float or None (default: None)
+                The azimuth at the target position.
+            degrees: boolean (default: True)
+                Indicates whether the input angles (inc, azi) are in degrees
+                (True) or radians (False).
+            unit: string (default: 'meters')
+                Indicates the distance unit, either 'meters' or 'feet'.
+            min_error: float (default: 1e-5):
+                Infers the error tolerance of the results and is used to set
+                iteration stops when the desired error tolerance is met. Value
+                must be less than 1. Use with caution as the code may
+                become unstable if this value is changed.
+            delta_radius: float (default: 10)
+                The delta radius (first curve and second curve sections) used
+                as an iteration stop when balancing radii. If the resulting
+                delta radius yielded from `dls_design` and `dls_design2` is
+                larger than `delta_radius`, then `delta_radius` defaults to
+                the former.
+            min_tangent: float (default: 10)
+                The minimum tangent length in the `curve_hold_curve` method
+                used to mitigate instability during iterations (where the
+                tangent section approaches or equals 0).
+            max_iterations: int (default: 1000)
+                The maximum number of iterations before giving up trying to
+                fit a `curve_hold_curve`. This number is limited by Python's
+                depth of recursion, but if you're hitting the default stop
+                then consider changing `delta_radius` and `min_tangent` as
+                your expectations may be unrealistic (this is open source
+                software after all!)     
+
+        Results
+        -------
+            connector: welleng.connector.Connector object
+        """
+
+        # TODO: remove self.step
+
+        # Set up a lookup dictionary to use with the logic to determine
+        # what connector method to deploy. Use a binary string to
+        # represent the inputs in the order:
+        # (md2, inc2, azi2, pos2, vec2)
+        # Initially I used boolean logic, but it quickly became non-
+        # transparent and difficult to debug.
+        self._get_initial_methods()
+
+        # METHODS = [
+        #     'hold',
+        #     'curve_hold', # 
+        #     'min_dist_to_target',
+        #     'min_curve_to_target',
+        #     'curve_hold_curve',
+        #     'min_curve'
+        # ]
+
+        # quick check that inputs are workable and if not some steer to 
+        # the user.
+        assert vec1 is not None or (inc1 is not None and azi1 is not None), "Require either vec1 or (inc1 and azi1)"
+        if vec1 is not None:
+            assert inc1 is None and azi1 is None, "Either vec1 or (inc1 and azi1)"
+        if (inc1 is not None or azi1 is not None):
+            assert vec1 is None, "Either vec1 or (inc1 and azi1)"
+
+        assert (
+            md2 is not None
+            or pos2 is not None
+            or vec2 is not None
+            or inc2 is not None
+            or azi2 is not None
+        ), "Missing target parameters"
+
+        if vec2 is not None:
+            assert not (inc2 or azi2), "Either vec2 or (inc2 and azi2)"
+        if (inc2 or azi2):
+            assert not vec2, "Either vec2 or (inc2 and azi2)"
+        if md2:
+            assert not pos2, "Either md2 or pos2"
+
+        assert dls_design > 0, "dls_design must be greater than zero"
+        assert min_error < 1, "min_error must be less than 1.0"
+
+        # figure out what method is required to connect the points
+        target_input = convert_target_input_to_booleans(md2, inc2, azi2, pos2, vec2)
+        self.initial_method = self.initial_methods[target_input]
+        
+        # do some more initialization stuff
+        self.min_error = min_error
+        self.min_tangent = min_tangent
+        self.iterations = 0
+        self.max_iterations = max_iterations
+        self.errors = []
+        self.dogleg_old, self.dogleg2_old = 0, 0
+        self.dist_curve2 = 0
+        self.pos1 = np.array(pos1)
+
+        self.pos2, self.vec2, self.inc2, self.azi2, self.md2 = (
+            None, None, None, None, None
+        )
+
+        # fill in the input data gaps
+        if (inc1 is not None and azi1 is not None):
+            if degrees:
+                self.inc1 = np.radians(inc1)
+                self.azi1 = np.radians(azi1)
+            else:
+                self.inc1 = inc1
+                self.azi1 = azi1
+            self.vec1 = np.array(get_vec(self.inc1, self.azi1, nev=True, deg=False))
+        else:
+            self.vec1 = np.array(vec1)
+            self.inc1, self.azi1 = get_angles(self.vec1, nev=True).reshape(2)
+
+        self.md1 = md1
+        self.pos_target = None if pos2 is None else np.array(pos2)
+        self.md_target = md2
+
+        if vec2 is not None:
+            self.vec_target = np.array(vec2)
+            self.inc_target, self.azi_target = get_angles(
+                self.vec_target,
+                nev=True
+            ).reshape(2)
+        elif (inc2 is not None and azi2 is not None):
+            if degrees:
+                self.inc_target = np.radians(inc2)
+                self.azi_target = np.radians(azi2)
+            else:
+                self.inc_target = inc2
+                self.azi_target = azi2
+            self.vec_target = get_vec(
+                self.inc_target, self.azi_target, nev=True, deg=False
+            ).reshape(3)
+        elif inc2 is None and azi2 is None:
+            self.inc_target, self.azi_target, self.vec_target = (
+                self.inc1, self.azi1, self.vec1
+            )
+        elif inc2 is None:
+            self.inc_target = self.inc1
+            if degrees:
+                self.azi_target = np.radians(azi2)
+            else:
+                self.azi_target = azi2
+            self.vec_target = get_vec(self.inc_target, self.azi_target, nev=True, deg=False)
+        elif azi2 is None:
+            self.azi_target = self.azi1
+            if degrees:
+                self.inc_target = np.radians(inc2)
+            else:
+                self.inc_target = inc2
+            self.vec_target = get_vec(self.inc_target, self.azi_target, nev=True, deg=False)
+        else:
+            self.vec_target = vec2
+            self.inc_target = inc2
+            self.azi_target = azi2
+
+        self.unit = unit
+        if self.unit == 'meters':
+            self.denom = 30
+        else:
+            self.denom = 100
+
+        if degrees:
+            self.dls_design = np.radians(dls_design)
+            if dls_design2:
+                self.dls_design2 = np.radians(dls_design2)
+        else:
+            self.dls_design = dls_design
+            if dls_design2:
+                self.dls_design2 = dls_design2
+        if not dls_design2:
+            self.dls_design2 = self.dls_design
+
+        self.radius_design = self.denom / self.dls_design
+        self.radius_design2 = self.denom / self.dls_design2
+
+        # check that the `delta_radius` actually makes sense
+        self.delta_radius = max(
+            delta_radius,
+            abs(self.radius_design - self.radius_design2)
+        )
+
+        # some more initialization stuff
+        self.tangent_length = None
+        self.dogleg2 = None
+
+        self.pos3, self.vec3, self.inc3, self.azi3, self.md3 = (
+            None, None, None, None, None
+        )
+        self.radius_critical, self.radius_critical2 = np.inf, np.inf
+
+        # Things fall apart if the start and end vectors exactly equal
+        # one another, so need to check for this and if this is the 
+        # case, modify the end vector slightly. This is a lazy way of
+        # doing this, but it's fast. Probably a more precise way would
+        # be to split the dogelg in two, but that's more hassle than
+        # it's worth.
+        if (
+            self.vec_target is not None
+            and np.array_equal(self.vec_target, self.vec1 * -1)
+        ):
+            (
+                self.vec_target,
+                self.inc_target,
+                self.azi_target
+            ) = mod_vec(self.vec_target, self.min_error)
+
+        # properly figure out the method
+        self._get_method()
+
+        # and finally, actually do something...
+        self._use_method()
+    
+    def _min_dist_to_target(self):
+        (
+            self.tangent_length,
+            self.dogleg
+        ) = min_dist_to_target(self.radius_design, self.distances)
+        self.dogleg = check_dogleg(self.dogleg)
+        self.dist_curve, self.func_dogleg = get_curve_hold_data(
+            self.radius_design, self.dogleg
+        )
+        self.vec_target = get_vec_target(
+            self.pos1,
+            self.vec1,
+            self.pos_target,
+            self.tangent_length,
+            self.dist_curve,
+            self.func_dogleg
+        )
+        self._get_angles_target()
+        self._get_md_target()
+        self.pos2 = (
+            self.pos_target - (
+                self.tangent_length * self.vec_target
+            )
+        )
+        self.md2 = self.md1 + abs(self.dist_curve)
+        self.md_target = self.md2 + self.tangent_length
+        self.vec2 = self.vec_target
+
+    def _min_curve_to_target(self):
+        (
+            self.tangent_length,
+            self.radius_critical,
+            self.dogleg
+        ) = min_curve_to_target(self.distances)
+        self.dogleg = check_dogleg(self.dogleg)
+        self.dist_curve, self.func_dogleg = get_curve_hold_data(
+            min(self.radius_design, self.radius_critical), self.dogleg
+        )
+        self.vec_target = get_vec_target(
+            self.pos1,
+            self.vec1,
+            self.pos_target,
+            self.tangent_length,
+            self.dist_curve,
+            self.func_dogleg
+        )
+        self._get_angles_target()
+        self._get_md_target()
+    
+    def _use_method(self):
+        if self.method == 'hold':
+            self._hold()
+        elif self.method == 'min_curve':
+            self._min_curve()
+        elif self.method == 'curve_hold_curve':
+            self.pos2_list, self.pos3_list = [], [self.pos_target]
+            self.vec23 = [np.array([0,0,0])]
+            self._target_pos_and_vec_defined(self.pos_target)
+        else:
+            self.distances = self._get_distances(
+                self.pos1, self.vec1, self.pos_target
+            )
+            if self.radius_design <= get_radius_critical(
+                self.radius_design, self.distances, self.min_error
+            ):
+                self.method = 'min_dist_to_target'
+                self._min_dist_to_target()
+            else:
+                self.method = 'min_curve_to_target'
+                self._min_curve_to_target()
+
+    
+    def _get_method(self):
+        assert self.initial_method not in [
+            'no_input',
+            'vec_and_inc_azi',
+            'md_and_pos'
+        ], f"{self.initial_method}"
+        if self.initial_method == 'hold':
+            self.method = 'hold'
+        elif self.initial_method[-8:] == '_or_hold':
+            if np.array_equal(self.vec_target, self.vec1):
+                self.method = 'hold'
+            else:
+                self.method = self.initial_method[:-8]
+        else:
+            self.method = self.initial_method
+    
+    def _get_initial_methods(self):
+        # TODO: probably better to load this in from a yaml file
+        self.initial_methods = {
+            '00000': 'no_input',
+            '00001': 'min_curve_or_hold',
+            '00010': 'curve_hold',
+            '00011': 'curve_hold_curve',
+            '00100': 'min_curve_or_hold',
+            '00101': 'vec_and_inc_azi',
+            '00110': 'curve_hold',
+            '00111': 'vec_and_inc_azi',
+            '01000': 'min_curve_or_hold',
+            '01001': 'vec_and_inc_azi',
+            '01010': 'curve_hold_or_hold',
+            '01011': 'vec_and_inc_azi',
+            '01100': 'min_curve_or_hold',
+            '01101': 'vec_and_inc_azi',
+            '01110': 'curve_hold_curve',
+            '01111': 'vec_and_inc_azi',
+            '10000': 'hold',
+            '10001': 'min_curve_or_hold',
+            '10010': 'md_and_pos',
+            '10011': 'md_and_pos',
+            '10100': 'min_curve_or_hold',
+            '10101': 'vec_and_inc_azi',
+            '10110': 'md_and_pos',
+            '10111': 'md_and_pos',
+            '11000': 'min_curve_or_hold',
+            '11001': 'vec_and_inc_azi',
+            '11010': 'md_and_pos',
+            '11011': 'md_and_pos',
+            '11100': 'min_curve',
+            '11101': 'vec_and_inc_azi',
+            '11110': 'md_and_pos',
+            '11111': 'md_and_pos'
+        }
+    
+    def _min_curve(self):
+        self.dogleg = get_dogleg(
+            self.inc1, self.azi1, self.inc_target, self.azi_target
+        )
+        self.dogleg = check_dogleg(self.dogleg)
+        if self.md_target is None:
+            self.md2 = None
+            self.dist_curve, self.func_dogleg = get_curve_hold_data(
+                        self.radius_design, self.dogleg
+                    )
+            self.md_target = self.md1 + abs(self.dist_curve)
+            self.pos_target = get_pos(
+                    self.pos1,
+                    self.vec1,
+                    self.vec_target,
+                    self.dist_curve,
+                    self.func_dogleg
+                ).reshape(3)
+        else:
+            self.radius_critical = abs((self.md_target - self.md1) / self.dogleg)
+            if (
+                self.radius_critical > self.radius_design
+                or (
+                    np.around(self.dogleg, decimals=5)
+                    == np.around(np.pi, decimals=5)
+                )
+            ):
+                self.md2 = (
+                    self.md1
+                    + min(self.radius_design, self.radius_critical) * self.dogleg
+                )
+                (
+                    self.inc2, self.azi2, self.vec2
+                ) = self.inc_target, self.azi_target, self.vec_target
+                self.dist_curve, self.func_dogleg = get_curve_hold_data(
+                        min(self.radius_design, self.radius_critical),
+                        self.dogleg
+                )
+                self.pos2 = get_pos(
+                    self.pos1,
+                    self.vec1,
+                    self.vec2,
+                    self.dist_curve,
+                    self.func_dogleg
+                ).reshape(3)
+                self.pos_target = self.pos2 + (
+                    self.vec2 * (self.md_target - self.md2)
+                )
+            else:    
+                self.dist_curve, self.func_dogleg = get_curve_hold_data(
+                        self.radius_critical, self.dogleg
+                    )
+                self.md2 = None
+                self.pos_target = get_pos(
+                    self.pos1,
+                    self.vec1,
+                    self.vec_target,
+                    self.dist_curve,
+                    self.func_dogleg
+                ).reshape(3)
+    
+    def _hold(self):        
+        self.pos_target = (
+            self.pos1 + self.vec1 * (self.md_target - self.md1)
+        )
+
+    def _get_angles_target(self):
+        self.inc_target, self.azi_target = get_angles(self.vec_target, nev=True).reshape(2)
+
+    def _get_md_target(self):
+        self.md_target = (
+            self.dist_curve
+            + self.tangent_length
+            + self.dist_curve2
+            + self.md1
+        )
+
+    def _target_pos_and_vec_defined(self, pos3, vec_old=[0,0,0]):
+        self.pos3 = pos3
+        self.distances1 = self._get_distances(self.pos1, self.vec1, self.pos3)
+
+        radius_temp1 = get_radius_critical(
+            self.radius_design,
+            self.distances1,
+            self.min_error
+        )
+        if radius_temp1 < self.radius_critical:
+            self.radius_critical = radius_temp1
+
+        (
+            self.tangent_length,
+            self.dogleg
+        ) = min_dist_to_target(
+            min(self.radius_design, self.radius_critical),
+            self.distances1
+        )
+
+        self.dogleg = check_dogleg(self.dogleg)
+
+        self.dist_curve, self.func_dogleg = get_curve_hold_data(
+                    min(self.radius_critical, self.radius_design),
+                    self.dogleg
+                )
+        self.vec3 = get_vec_target(
+            self.pos1,
+            self.vec1,
+            self.pos3,
+            self.tangent_length,
+            self.dist_curve,
+            self.func_dogleg
+        )
+
+        tangent_temp1 = self._get_tangent_temp(self.tangent_length)
+
+        self.pos2 = (
+            self.pos3 - (
+                tangent_temp1 * self.vec3
+            )
+        )
+
+        self.distances2 = self._get_distances(
+            self.pos_target,
+            self.vec_target * -1,
+            self.pos2
+        )
+
+        radius_temp2 = get_radius_critical(
+            self.radius_design2,
+            self.distances2,
+            self.min_error
+        )
+        if radius_temp2 < self.radius_critical2:
+            self.radius_critical2 = radius_temp2
+
+        (
+            self.tangent_length2,
+            self.dogleg2
+        ) = min_dist_to_target(
+            min(self.radius_design2, self.radius_critical2),
+            self.distances2
+        )
+
+        self.dogleg2 = check_dogleg(self.dogleg2)
+
+        self.dist_curve2, self.func_dogleg2 = get_curve_hold_data(
+                    min(self.radius_critical2, self.radius_design2),
+                    self.dogleg2
+                )
+        self.vec2 = get_vec_target(
+            self.pos_target,
+            self.vec_target * -1,
+            self.pos2,
+            self.tangent_length2,
+            self.dist_curve2,
+            self.func_dogleg2
+        )
+
+        tangent_temp2 = self._get_tangent_temp(self.tangent_length2)
+
+        self.pos3 = (
+            self.pos2 - (
+                tangent_temp2 * self.vec2
+            )
+        )
+
+        self.vec23.append((self.pos3 - self.pos2) / np.linalg.norm(self.pos3 - self.pos2))
+
+        self.error = np.allclose(
+            self.vec23[-1],
+            vec_old,
+            equal_nan=True,
+            rtol=self.min_error * 10,
+            atol=self.min_error * 0.1
+        )
+
+        self.errors.append(self.error)
+        self.pos3_list.append(self.pos3)
+        self.pos2_list.append(self.pos2)
+        self.md_target = self.dist_curve + tangent_temp2 + self.dist_curve2
+
+        # a bit of recursive magic - proceed with caution
+        if (
+            not self.error
+            or
+            (
+                (
+                    self.radius_critical2 < self.radius_design2
+                    or self.radius_critical < self.radius_design
+                )
+                and
+                abs(self.radius_critical - self.radius_critical2) > self.delta_radius
+            )
+        ):
+            # A solution will typically be found within a few iterations,
+            # however it will likely be unbalanced in terms of dls between the
+            # two curve sections. The code below will re-initiate iterations
+            # to balance the curvatures until the delta_radius parameter is met.
+            if self.error:
+                self.radius_critical = self.md_target / (
+                    abs(self.dogleg) + abs(self.dogleg2)
+                )
+                self.radius_critical2 = self.radius_critical
+            self.iterations += 1
+            # prevent a loop ad infinitum
+            assert self.iterations < self.max_iterations, "Out of iteration ammo"
+            self._target_pos_and_vec_defined(self.pos3, self.vec23[-1])
+        else:
+            # get pos1 to pos2 curve data
+            self.dist_curve, self.func_dogleg = get_curve_hold_data(
+                min(self.radius_critical, self.radius_design),
+                self.dogleg
+            )
+
+            self.vec3 = get_vec_target(
+                self.pos1,
+                self.vec1,
+                self.pos3,
+                self.tangent_length,
+                self.dist_curve,
+                self.func_dogleg
+            )
+
+            self.vec2 = self.vec3
+
+            self.pos2 = get_pos(
+                self.pos1,
+                self.vec1,
+                self.vec3,
+                self.dist_curve,
+                self.func_dogleg
+            ).reshape(3)
+
+            self.md2 = self.md1 + abs(self.dist_curve)
+
+            self.dist_curve2, self.func_dogleg2 = get_curve_hold_data(
+                min(self.radius_critical2, self.radius_design2),
+                self.dogleg2
+            )
+
+            self.pos3 = get_pos(
+                self.pos_target,
+                self.vec_target * -1,
+                self.vec3 * -1,
+                self.dist_curve2,
+                self.func_dogleg2
+            ).reshape(3)
+
+            tangent_length = np.linalg.norm(
+                self.pos3 - self.pos2
+            )
+
+            self.md3 = self.md2 + tangent_length
+
+            self.md_target = self.md3 + abs(self.dist_curve2)
+    
+    def interpolate(self, step=30):
+        return interpolate_well([self], step)
+
+    def survey(self, radius=10, step=30):
+        interpolation = self.interpolate(step)
+
+        survey = get_survey(interpolation, start_nev=self.pos1, radius=10)
+
+        return survey
+
+    def _get_tangent_temp(self, tangent_length):
+        if np.isnan(tangent_length):
+            tangent_temp = self.min_tangent
+        else:
+            tangent_temp = max(tangent_length, self.min_tangent)
+
+        return tangent_temp
+
+    def _mod_pos(self, pos):
+        pos_rand = np.random.random(3) * self.delta_radius
+        pos += pos_rand
+
+    def _get_distances(self, pos1, vec1, pos_target):
+        # When initializing a `curve_hold_curve` and pos_target is directly
+        # ahead it can cause issues (it's a hold and not a curve). So this code
+        # checks for that condition and if it's the case, will move the
+        # target_pos a sufficient amount to prevent issues.
+        vec_temp = np.array(pos_target - pos1)
+        vec_temp = vec_temp / np.linalg.norm(vec_temp)
+        if np.array_equal(vec1, vec_temp):
+            self._mod_pos(pos_target)
+        if np.allclose(pos1, pos_target):
+            return (0, 0, 0)
+        
+        else:
+            dist_to_target = (
+                np.linalg.norm((pos_target - pos1))
+            )
+
+            dist_perp_to_target = (
+                np.dot((pos_target - pos1), vec1)
+            )
+
+            dist_norm_to_target = (
+                (
+                    dist_to_target ** 2
+                    - dist_perp_to_target ** 2
+                ) ** 0.5
+            )
+
+            return (
+                dist_to_target,
+                dist_perp_to_target,
+                dist_norm_to_target
+            )
+
+
+def check_dogleg(dogleg):
+    # the code assumes angles are positive and clockwise
+    if dogleg < 0:
+        dogleg_new = dogleg + 2 * np.pi
+        return dogleg_new
+    else:
+        return dogleg
+
+def mod_vec(vec, error=1e-5):
+    # if it's not working then twat it with a hammer
+    vec_mod = vec * np.array([1, 1, 1 - error])
+    vec_mod /= np.linalg.norm(vec_mod)
+    inc_mod, azi_mod = get_angles(vec_mod, nev=True).T
+
+    return vec_mod, inc_mod, azi_mod
+
+def get_pos(pos1, vec1, vec2, dist_curve, func_dogleg):
+    # this function is redundent, but I'd like to reinstate it one day
+    # TODO: this seems to be a bit dodgy if you send it an array for
+    # vec2 - investigate.
+    inc1, azi1 = get_angles(vec1, nev=True).T
+    inc2, azi2 = get_angles(vec2, nev=True).T
+
+    pos2 = pos1 + (
+        (
+            dist_curve * func_dogleg / 2
+        )
+        *
+        np.array([
+            np.sin(inc1) * np.cos(azi1) + np.sin(inc2) * np.cos(azi2),
+            np.sin(inc1) * np.sin(azi1) + np.sin(inc2) * np.sin(azi2),
+            np.cos(inc1) + np.cos(inc2)
+        ]).T
+    )
+
+    return pos2
+
+def get_vec_target(
+    pos1,
+    vec1,
+    pos_target,
+    tangent_length,
+    dist_curve,
+    func_dogleg
+):
+    vec_target = (
+        (
+            pos_target - pos1 - (
+                    dist_curve
+                    * func_dogleg
+                ) / 2 * vec1
+        )
+        /
+        (
+            (
+                dist_curve * func_dogleg / 2
+            ) + tangent_length
+        )
+    )
+    
+    vec_target /= np.linalg.norm(vec_target)
+
+    return vec_target
+
+def get_curve_hold_data(radius, dogleg):
+    dist_curve = radius * dogleg
+    func_dogleg = shape_factor(dogleg)
+
+    return (
+        dist_curve,
+        func_dogleg
+    )
+
+def shape_factor(dogleg):
+    """
+    Function to determine the shape factor of a dogleg.
+
+    Parameters
+    ----------
+        dogleg: float
+            The dogleg angle in radians of a curve section.
+    """
+    return np.tan(dogleg / 2) / (dogleg / 2)
+
+def min_dist_to_target(radius, distances):
+    """
+    Calculates the control points for a curve and hold section from the
+    start position and vector to the target position.
+    """
+    (
+        dist_to_target,
+        dist_perp_to_target,
+        dist_norm_to_target
+    ) = distances
+
+    tangent_length = (
+        dist_to_target ** 2
+        - 2 * radius * dist_norm_to_target
+    ) ** 0.5
+
+    # determine the dogleg angle of the curve section
+    dogleg = 2 * np.arctan2(
+        (dist_perp_to_target - tangent_length)
+        ,
+        (
+            2 * radius - dist_norm_to_target
+        )
+    )
+
+    return tangent_length, dogleg
+
+def min_curve_to_target(distances):
+    """
+    Calculates the control points for a curve section from the start 
+    position and vector to the target position which is not achievable with
+    the provided dls_design.
+    """
+    (
+        dist_to_target,
+        dist_perp_to_target,
+        dist_norm_to_target
+    ) = distances
+
+    radius_critical = (
+        dist_to_target ** 2 / (
+            2 * dist_norm_to_target
+        )
+    )
+
+    dogleg = (
+        2 * np.arctan2(
+            dist_norm_to_target,
+            dist_perp_to_target
+        )
+    )
+
+    tangent_length = 0
+
+    return (
+        tangent_length,
+        radius_critical,
+        dogleg
+    )
+
+def get_radius_critical(radius, distances, min_error):
+    (
+        dist_to_target,
+        dist_perp_to_target,
+        dist_norm_to_target
+    ) = distances
+
+    radius_critical = (
+        dist_to_target ** 2 / (
+            2 * dist_norm_to_target
+        )
+    ) * (1 - min_error)
+
+    return radius_critical
+
+def angle(vec1, vec2, acute=True):
+    angle = np.arccos(
+        np.dot(vec1, vec2)
+        /
+        (
+            np.linalg.norm(vec1) * np.linalg.norm(vec2)
+        )
+    )
+
+    if acute:
+        return angle
+
+    else:
+        return 2 * np.pi - angle
+
+def get_dogleg(inc1, azi1, inc2, azi2):
+    dogleg = (
+        2 * np.arcsin(
+            (
+                np.sin((inc2 - inc1) / 2) ** 2
+                + np.sin(inc1) * np.sin(inc2)
+                * np.sin((azi2 - azi1) / 2) ** 2
+            ) ** 0.5
+        )
+    )
+
+    return dogleg
+
+def interpolate_well(sections, step=30):
+    """
+    Constructs a well survey from a list of sections of control points.
+
+    Parameters
+    ----------
+        sections: list of welleng.connector.Connector objects
+        step: float (default: 30)
+            The desired delta measured depth between survey points, in
+            addition to the control points.
+
+    Results
+    -------
+        survey: `welleng.survey.Survey` object
+    """
+    method = {
+        'hold': get_interpolate_hold,
+        'min_dist_to_target': get_interpolate_min_dist_to_target,
+        'min_curve_to_target': get_interpolate_min_curve_to_target,
+        'curve_hold_curve': get_interpololate_curve_hold_curve,
+        'min_curve': get_min_curve
+    }
+    
+    data = []
+    for s in sections:
+        data.extend(method[s.method](s, step))
+
+    return data
+
+def get_survey(section_data, start_nev=[0,0,0], radius=10):
+    """
+    Constructs a well survey from a list of sections of control points.
+
+    Parameters
+    ----------
+        section_data: list of dicts with section data
+        start_nev: (3) array of floats (default: [0,0,0])
+            The starting position in NEV coordinates.
+        radius: float (default: 10)
+            The radius is passed to the `welleng.survey.Survey` object
+            and represents the radius of the wellbore. It is also used
+            when visualizing the results, so can be used to make the
+            wellbore *thicker* in the plot.
+
+    Results
+    -------
+        survey: `welleng.survey.Survey` object
+    """
+
+    # generate lists for survey
+    md, inc, azi = np.vstack([np.array(list(zip(
+            s['md'].tolist(),
+            s['inc'].tolist(),
+            s['azi'].tolist()
+        )))
+        for s in section_data
+    ]).T
+
+    survey = Survey(
+        md=md,
+        inc=inc,
+        azi=azi,
+        start_nev=start_nev,
+        deg=False,
+        radius=radius,
+    )
+
+    return survey
+
+def interpolate_curve(md1, pos1, vec1, vec2, dist_curve, dogleg, func_dogleg, step, endpoint=False):
+    start_md = step - (md1 % step)
+    end_md = abs(dist_curve)
+    md = np.arange(start_md, end_md, step)
+    md = np.concatenate(([0.], md))
+    if endpoint:
+        md = np.concatenate((md, [end_md]))
+    dogleg_interp = (dogleg / dist_curve * md).reshape(-1,1)
+
+    vec = (
+        (
+            np.sin(dogleg - dogleg_interp) / np.sin(dogleg) * vec1
+        )
+        +
+        (
+            np.sin(dogleg_interp) / np.sin(dogleg) * vec2
+        )
+    )
+    vec = vec / np.linalg.norm(vec, axis=1).reshape(-1,1)
+    inc, azi = get_angles(vec, nev=True).T
+
+    data = dict(
+        md = md + md1,
+        vec = vec,
+        inc = inc,
+        azi = azi,
+        dogleg = np.concatenate((
+            np.array([0.]), np.diff(dogleg_interp.reshape(-1))
+        )),
+    )
+
+    return data
+
+def interpolate_hold(md1, pos1, vec1, md2, step, endpoint=False):
+    start_md = step - (md1 % step)
+    end_md = md2 - md1
+    md = np.arange(start_md, end_md, step)
+    md = np.concatenate(([0.], md)) 
+    if endpoint:
+        md = np.concatenate((md, [end_md]))
+    vec = np.full((len(md),3), vec1)
+    inc, azi = get_angles(vec, nev=True).T
+    dogleg = np.full_like(md, 0.)
+
+    data = dict(
+        md = md + md1,
+        vec = vec,
+        inc = inc,
+        azi = azi,
+        dogleg = dogleg,
+    )
+
+    return data
+
+def get_min_curve(section, step=30, data=None):
+    if section.md2 is None:
+        result = (
+            get_interpolate_min_curve_to_target(
+                section, step, data
+            )
+        )
+    else:
+        result = (
+            get_interpolate_min_dist_to_target(
+                section, step, data
+            )
+        )
+    return result
+
+def get_interpolate_hold(section, step=30, data=None):
+    if data is None:
+        data = []
+
+    data.append(interpolate_hold(
+        md1=section.md1,
+        pos1=section.pos1,
+        vec1=section.vec1,
+        md2=section.md_target,
+        step=step,
+        endpoint=True
+    ))
+
+    return data
+
+def get_interpolate_min_curve_to_target(section, step=30, data=None):
+    if data is None:
+        data = []
+
+    data.append(interpolate_curve(
+        md1=section.md1,
+        pos1=section.pos1,
+        vec1=section.vec1,
+        vec2=section.vec_target,
+        dist_curve=section.dist_curve,
+        dogleg=section.dogleg,
+        func_dogleg=section.func_dogleg,
+        step=step,
+        endpoint=True
+    ))
+
+    return data  
+
+def get_interpolate_min_dist_to_target(section, step=30, data=None):
+    if data is None:
+        data = []
+
+    # the first curve section
+    data.append(interpolate_curve(
+        md1=section.md1,
+        pos1=section.pos1,
+        vec1=section.vec1,
+        vec2=section.vec2,
+        dist_curve=section.dist_curve,
+        dogleg=section.dogleg,
+        func_dogleg=section.func_dogleg,
+        step=step
+    ))
+
+    # the hold section
+    data.append(interpolate_hold(
+        md1=section.md2,
+        pos1=section.pos2,
+        vec1=section.vec2,
+        md2=section.md_target,
+        step=step,
+        endpoint=True
+    ))
+
+    return data    
+
+def get_interpololate_curve_hold_curve(section, step=30, data=None):
+    if data is None:
+        data = []
+
+    # the first curve section
+    data.append(interpolate_curve(
+        md1=section.md1,
+        pos1=section.pos1,
+        vec1=section.vec1,
+        vec2=section.vec2,
+        dist_curve=section.dist_curve,
+        dogleg=section.dogleg,
+        func_dogleg=section.func_dogleg,
+        step=step
+    ))
+
+    # the hold section
+    data.append(interpolate_hold(
+        md1=section.md2,
+        pos1=section.pos2,
+        vec1=section.vec2,
+        md2=section.md3,
+        step=step
+    ))
+
+    # the second curve section
+    data.append(interpolate_curve(
+        md1=section.md3,
+        pos1=section.pos3,
+        vec1=section.vec3,
+        vec2=section.vec_target,
+        dist_curve=section.dist_curve2,
+        dogleg=section.dogleg2,
+        func_dogleg=section.func_dogleg2,
+        step=step,
+        endpoint=True
+    ))
+
+    return data
+
+def convert_target_input_to_booleans(*inputs):
+    input = [
+        "0" if i is None else "1" for i in inputs
+    ]
+
+    return ''.join(input)
\ No newline at end of file
diff --git a/welleng/errors/__pycache__/__init__.cpython-38.pyc b/welleng/errors/__pycache__/__init__.cpython-38.pyc
index 5d74a52..e97aa24 100644
Binary files a/welleng/errors/__pycache__/__init__.cpython-38.pyc and b/welleng/errors/__pycache__/__init__.cpython-38.pyc differ
diff --git a/welleng/errors/__pycache__/iscwsa_mwd.cpython-38.pyc b/welleng/errors/__pycache__/iscwsa_mwd.cpython-38.pyc
index c7e8350..5c7e681 100644
Binary files a/welleng/errors/__pycache__/iscwsa_mwd.cpython-38.pyc and b/welleng/errors/__pycache__/iscwsa_mwd.cpython-38.pyc differ
diff --git a/welleng/exchange/__init__.py b/welleng/exchange/__init__.py
new file mode 100644
index 0000000..99585fc
--- /dev/null
+++ b/welleng/exchange/__init__.py
@@ -0,0 +1,6 @@
+ 
+"""
+welleng/exchange
+----------------
+Contains the importers and exporters for various survey formats.
+"""
\ No newline at end of file
diff --git a/welleng/exchange/__pycache__/__init__.cpython-37.pyc b/welleng/exchange/__pycache__/__init__.cpython-37.pyc
new file mode 100644
index 0000000..3f87ba3
Binary files /dev/null and b/welleng/exchange/__pycache__/__init__.cpython-37.pyc differ
diff --git a/welleng/exchange/__pycache__/__init__.cpython-38.pyc b/welleng/exchange/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000..c3a2057
Binary files /dev/null and b/welleng/exchange/__pycache__/__init__.cpython-38.pyc differ
diff --git a/welleng/exchange/__pycache__/wbp.cpython-37.pyc b/welleng/exchange/__pycache__/wbp.cpython-37.pyc
new file mode 100644
index 0000000..683ae25
Binary files /dev/null and b/welleng/exchange/__pycache__/wbp.cpython-37.pyc differ
diff --git a/welleng/exchange/__pycache__/wbp.cpython-38.pyc b/welleng/exchange/__pycache__/wbp.cpython-38.pyc
new file mode 100644
index 0000000..fc9ecc0
Binary files /dev/null and b/welleng/exchange/__pycache__/wbp.cpython-38.pyc differ
diff --git a/welleng/exchange/wbp.py b/welleng/exchange/wbp.py
new file mode 100644
index 0000000..303b4a0
--- /dev/null
+++ b/welleng/exchange/wbp.py
@@ -0,0 +1,165 @@
+import numpy as np
+from ..survey import get_sections
+
+class WellPlan:
+    def __init__(
+        self,
+        survey,
+        location_type,
+        plan_method,
+        plan_name,
+        parent_name,
+        depth_unit="meters",
+        surface_unit="meters",
+        dirty_flag=3,
+        sidetrack_id=0,
+        extension="0.00000",
+        **targets
+    ):
+
+        LOCATIONS = [
+            "unknown",
+            "surface",
+            "ow_sidetrack",
+            "wp_sidetrack",
+            "lookahead",
+            "ow_well",
+            "complex_extension",
+            "site",
+            "site_based_plan",
+            "compass_well"
+        ]
+
+        METHODS = [
+            "curve_only",
+            "curve_hold",
+            "opt_align"
+        ]
+
+        UNITS = [
+            "feet",
+            "meters"
+        ]
+
+        assert location_type in LOCATIONS, "location_type not in LOCATIONS"
+        assert plan_method in METHODS, "plan_method not in METHODS"
+        assert len(plan_name) < 21, "plan_name too long (max length 20)"
+        assert len(parent_name) < 21, "parent_name too long (max length 20)"
+        assert surface_unit in UNITS, "surface_units must be feet or meters"
+        assert depth_unit in UNITS, "depth_units must be feet or meters"
+
+        self.survey = survey
+        self.location = str(LOCATIONS.index(location_type))
+        self.surface_unit = surface_unit
+        self.depth_unit = depth_unit
+        self.method = str(METHODS.index(plan_method))
+        self.flag = str(dirty_flag)
+        self.st_id = str(sidetrack_id).rjust(8)
+        self.plan_name = str(plan_name).ljust(21)
+        self.parent_name = str(parent_name).ljust(20)
+        self.dls = f"{str(np.around(self.survey.dls[1]))[:5]}".rjust(5)
+        self.kickoff_dls = f"{survey.dls[1]:.2f}".rjust(7)
+        self.extension = str(extension).rjust(8)
+        self.targets = targets
+
+        self._get_unit()
+
+        self.doc = []
+
+        self._make_header()
+        # self._make_tie_on(target=0)
+        self._make_wellplan()
+
+    def _get_unit(self):
+        if self.depth_unit == "meters":
+            if self.surface_unit == "meters":
+                self.unit = 2
+            else:
+                self.unit = 3
+        elif self.depth_unit == "feet":
+            if self.surface_unit == "feet":
+                self.unit = 1
+            else:
+                self.unit = 4
+
+    def _add_plan(
+        self,
+        section
+    ):
+        self.doc.append(
+            f"P:"
+            f"{str(section.md)[:8].rjust(8)}"
+            f" {str(section.azi)[:7].rjust(7)}"
+            f" {str(section.inc)[:7].rjust(7)}"
+            f" {str(section.build_rate)[:7].rjust(7)}"
+            f" {str(section.turn_rate)[:7].rjust(7)}"
+            f" {str(np.around(section.dls, decimals=1))[:7].rjust(7)}"
+            f" {str(np.degrees(section.toolface))[:7].rjust(7)}"
+            f" {str(section.method).rjust(4)}"
+            f" {str(section.target).rjust(10)}"
+        )
+        self.doc.append(
+            f"L:"
+            f"{str(section.x)[:13].rjust(13)}"
+            f"{str(section.y)[:13].rjust(13)}"
+            f"{str(section.z)[:11].rjust(11)}"
+        )
+
+
+    def _make_wellplan(
+        self,
+    ):
+        sections = get_sections(self.survey)
+
+        for s in sections:
+            self._add_plan(s)
+    
+    def _make_header(self):
+        self.doc.append(
+            f"DEPTH {self.unit}"
+        )
+
+        if self.targets:
+            self.doc.append(
+                "TARGETS:"
+            )
+            # TODO: add function to retrieve and list target data and implement
+            # a Target class.
+
+        self.doc.append(
+            "WELLPLANS:",
+        )
+        self.doc.append(
+            (
+                f"W:{self.location}{self.unit}{self.method}{self.flag}   0   0"
+                f"{self.st_id} {self.plan_name} {self.parent_name} {self.dls}"
+                f"{self.extension}{self.kickoff_dls}"
+            )
+        )
+
+    # I think this is redundent as it's the same line as the first section
+    def _make_tie_on(self, target):
+        if self.location == "unknown":
+            self.doc.append(
+                f"P:0.000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000    0"
+            )
+        else:
+            self.doc.append((
+                f"P:"
+                f"{str(self.survey.md[0])[:8].rjust(8)}"
+                f" {str(self.survey.azi_deg[0])[:7].rjust(7)}"
+                f" {str(self.survey.inc_deg[0])[:7].rjust(7)}"
+                f" {str(0.00000).rjust(7)}" # build rate
+                f" {str(0.00000).rjust(7)}" # turn rate
+                f" {str(0.00000).rjust(7)}" # DLS
+                f" {str(np.degrees(self.survey.toolface[0]))[:7].rjust(7)}"
+                f" {self.method.rjust(4)}"
+                f" {str(target).rjust(10)}"
+            ))
+
+            self.doc.append(
+                f"L:"
+                f"{str(self.survey.x)[:13].rjust(13)}"
+                f"{str(self.survey.y)[:13].rjust(13)}"
+                f"{str(self.survey.z)[:11].rjust(11)}"
+            )
diff --git a/welleng/mesh.py b/welleng/mesh.py
index f7fc9df..a5d83b0 100644
--- a/welleng/mesh.py
+++ b/welleng/mesh.py
@@ -14,7 +14,7 @@ def __init__(
         sigma=3.0,
         sigma_pa=0.5,
         Sm=0,
-        method="mesh_ellipse",
+        method="ellipse",
     ):
         """
         Create a WellMesh object from a welleng Survey object and a
@@ -50,11 +50,12 @@ def __init__(
         self.Sm = Sm
         self.sigma_pa = sigma_pa
 
-        assert method in ["mesh_ellipse", "mesh_pedal_curve"], \
+        assert method in ["ellipse", "pedal_curve", "circle"], \
             "Invalid method (ellipse or pedal_curve)"
         self.method = method
 
-        self.sigmaH, self.sigmaL, self.sigmaA = get_sigmas(self.s.cov_hla)
+        if self.method != 'circle':
+            self.sigmaH, self.sigmaL, self.sigmaA = get_sigmas(self.s.cov_hla)
         self.nevs = np.array([self.s.n, self.s.e, self.s.tvd]).T
         self._get_vertices()
         self._align_verts()
@@ -135,22 +136,26 @@ def _get_vertices(
         Determine the positions of the vertices on the desired uncertainty
         circumference.
         '''
+        if self.method == "circle":
+            h = self.s.radius.reshape(-1,1)
+            l = h
 
-        h = (
-            np.array(self.sigmaH) * self.sigma
-            + self.radius + self.Sm
-            + self.sigma_pa / 2
-        ).reshape(-1,1)
-
-        l = (
-            np.array(self.sigmaL) * self.sigma
-            + self.radius
-            + self.Sm
-            + self.sigma_pa / 2
-        ).reshape(-1,1)
-        # a = self.s.sigmaA * self.c.k + self.s.radius + self.c.Sm
-
-        if self.method == "mesh_ellipse":
+        else:
+            h = (
+                np.array(self.sigmaH) * self.sigma
+                + self.radius + self.Sm
+                + self.sigma_pa / 2
+            ).reshape(-1,1)
+
+            l = (
+                np.array(self.sigmaL) * self.sigma
+                + self.radius
+                + self.Sm
+                + self.sigma_pa / 2
+            ).reshape(-1,1)
+            # a = self.s.sigmaA * self.c.k + self.s.radius + self.c.Sm
+
+        if self.method in ["ellipse", "circle"]:
             lam = np.linspace(0, 2 * pi, self.n_verts, endpoint=False)
             # theta = np.zeros_like(lam)
 
diff --git a/welleng/survey.py b/welleng/survey.py
index eb968aa..a55118e 100644
--- a/welleng/survey.py
+++ b/welleng/survey.py
@@ -138,6 +138,7 @@ def __init__(
         self.vec = vec
 
         self._min_curve()
+        self._get_toolface_and_rates()
 
         # initialize errors
         # TODO: read this from a yaml file in errors
@@ -172,6 +173,8 @@ def _min_curve(self):
         self.rf = mc.rf
         self.delta_md = mc.delta_md
         self.dls = mc.dls
+        self.pos = mc.poss
+
         if self.x is None:
             # self.x, self.y, self.z = (mc.poss + self.start_xyz).T
             self.x, self.y, self.z = (mc.poss).T
@@ -242,6 +245,71 @@ def _get_errors(self):
             self.cov_nev += self.start_cov_nev
             self.cov_hla = NEV_to_HLA(self.survey_rad, self.cov_nev.T).T
 
+    def _curvature_to_rate(self, curvature):
+        with np.errstate(divide='ignore', invalid='ignore'):
+            radius = 1 / curvature
+        circumference = 2 * np.pi * radius
+        if self.unit == 'meters':
+            x = 30
+        else:
+            x = 100
+        rate = np.absolute(np.degrees(2 * np.pi / circumference) * x)
+
+        return rate
+    
+    def _get_toolface_and_rates(self):
+        """
+        Reference SPE-84246.
+        theta is inc, phi is azi
+        """
+        # split the survey
+        s = SplitSurvey(self)
+
+        if self.unit == 'meters':
+            x = 30
+        else:
+            x = 100
+
+        # this is lazy I know, but I'm using this mostly for flags
+        with np.errstate(divide='ignore', invalid='ignore'):
+            t1 = np.arctan(
+                np.sin(s.inc2) * np.sin(s.delta_azi) /
+                (
+                    np.sin(s.inc2) * np.cos(s.inc1) * np.cos(s.delta_azi)
+                    - np.sin(s.inc1) * np.cos(s.inc2)
+                )
+            )
+            t1 = np.nan_to_num(
+                np.where(t1 < 0, t1 + 2 * np.pi, t1),
+                nan=np.nan
+            )
+            t2 = np.arctan(
+                np.sin(s.inc1) * np.sin(s.delta_azi) /
+                (
+                    np.sin(s.inc2) * np.cos(s.inc1)
+                    - np.sin(s.inc1) * np.cos(s.inc2) * np.cos(s.delta_azi)
+                )
+            )
+            t2 = np.nan_to_num(
+                np.where(t2 < 0, t2 + 2 * np.pi, t2),
+                nan=np.nan
+            )
+            self.curve_radius = (360 / self.dls * x) / (2 * np.pi)
+        
+            curvature_dls = 1 / self.curve_radius
+
+            self.toolface = np.concatenate((t1, np.array([t2[-1]])))
+
+            curvature_turn = curvature_dls * (np.sin(self.toolface) / np.sin(self.inc_rad))
+            self.turn_rate = self._curvature_to_rate(curvature_turn)
+
+            curvature_build = curvature_dls * np.cos(self.toolface)
+            self.build_rate = self._curvature_to_rate(curvature_build)
+
+        # calculate plan normals
+        n12 = np.cross(s.vec1, s.vec2)
+        with np.errstate(divide='ignore', invalid='ignore'):
+            self.normals = n12 / np.linalg.norm(n12, axis=1).reshape(-1,1)
 
 def interpolate_survey(survey, x=0, index=0):
     """
@@ -398,4 +466,158 @@ def make_long_cov(arr):
         [ac, bc, cc]
     ]).T
 
-    return cov
\ No newline at end of file
+    return cov
+
+class SplitSurvey:
+    def __init__(
+        self,
+        survey,
+    ):
+        self.md1, self.inc1, self.azi1 = survey.survey_rad[:-1].T
+        self.md2, self.inc2, self.azi2 = survey.survey_rad[1:].T
+        self.delta_azi = self.azi2 - self.azi1
+        self.delta_inc = self.inc2 - self.inc1
+
+        self.vec1 = survey.vec[:-1]
+        self.vec2 = survey.vec[1:]
+        self.dogleg = survey.dogleg[1:]   
+
+def get_circle_radius(survey, **targets):
+    # TODO: add target data to sections
+    ss = SplitSurvey(survey)
+
+    x1, y1, z1 = np.cross(ss.vec1, survey.normals).T
+    x2, y2, z2 = np.cross(ss.vec2, survey.normals).T
+
+    b1 = np.array([y1, x1, z1]).T
+    b2 = np.array([y2, x2, z2]).T
+    nev = np.array([survey.n, survey.e, survey.tvd]).T
+
+    cc1 = (
+        nev[:-1] - b1
+        / np.linalg.norm(b1, axis=1).reshape(-1,1)
+        * survey.curve_radius[:-1].reshape(-1,1)
+    )
+    cc2 = (
+        nev[1:] - b2
+        / np.linalg.norm(b2, axis=1).reshape(-1,1)
+        * survey.curve_radius[1:].reshape(-1,1)
+    )
+
+    starts = np.vstack((cc1, cc2))
+    ends = np.vstack((nev[:-1], nev[1:]))
+
+    n = 1
+
+
+    return (starts, ends)
+
+def get_sections(survey, **targets):
+    # TODO: add target data to sections
+    ss = SplitSurvey(survey)
+    
+    continuous = np.all(
+        np.isclose(
+            survey.normals[:-1],
+            survey.normals[1:],
+            rtol=1e-01, atol=1e-02,
+            equal_nan=True
+        ), axis=-1
+    )
+
+    ends = np.concatenate(
+        (np.where(continuous == False)[0] + 1, np.array([len(continuous)-1]))
+    )
+    starts = np.concatenate(([0], ends[:-1]))
+    
+    actions = [
+        "hold" if a else "curve"
+        for a in np.isnan(survey.toolface[starts])
+    ]
+
+    sections = []
+    for s, e, a in zip(starts, ends, actions):
+        md = survey.md[s]
+        inc = survey.inc_deg[s]
+        azi = survey.azi_deg[s]
+        x = survey.e[s]
+        y = survey.n[s]
+        z = -survey.tvd[s]
+
+        target = ""
+        if survey.unit == 'meters':
+            denominator = 30
+        else:
+            denominator = 100
+
+        if a == "hold":
+            dls = 0.0
+            toolface = 0.0
+            build_rate = 0.0
+            turn_rate = 0.0
+            method = ""
+        else:
+            dls = survey.dls[e]
+            toolface = survey.toolface[s]
+            delta_md = survey.md[e] - md
+
+            # TODO: should sum this line by line to avoid issues with long sections
+            build_rate = abs(
+                (survey.inc_deg[e] - survey.inc_deg[s])
+                / delta_md * denominator
+            )
+
+            # TODO: should sum this line by line to avoid issues with long sections
+            # need to be careful with azimuth straddling north
+            delta_azi_1 = abs(survey.azi_deg[e] - survey.azi_deg[s])
+            delta_azi_2 = abs(360 - delta_azi_1)
+            delta_azi = min(delta_azi_1, delta_azi_2)
+            turn_rate = delta_azi / delta_md * denominator
+
+        section = WellSection(
+            md,
+            inc,
+            azi,
+            build_rate,
+            turn_rate,
+            dls,
+            toolface,
+            method=0,
+            target="",
+            x=x,
+            y=y,
+            z=z
+        )
+
+        sections.append(section)
+    
+    return sections
+
+class WellSection:
+    def __init__(
+        self,
+        md,
+        inc,
+        azi,
+        build_rate,
+        turn_rate,
+        dls,
+        toolface,
+        x,
+        y,
+        z,
+        method=0,
+        target=None
+    ):
+        self.md = md
+        self.inc = inc
+        self.azi = azi
+        self.build_rate = build_rate
+        self.turn_rate = turn_rate
+        self.dls = dls
+        self.toolface = toolface
+        self.method = method
+        self.target = target
+        self.x = np.around(x, decimals=2)
+        self.y = np.around(y, decimals=2)
+        self.z = np.around(z, decimals=2)
\ No newline at end of file
diff --git a/welleng/target.py b/welleng/target.py
new file mode 100644
index 0000000..bf49df8
--- /dev/null
+++ b/welleng/target.py
@@ -0,0 +1,70 @@
+import numpy as np
+from vedo import Circle
+
+class Target:
+    def __init__(
+        self,
+        name,
+        n,
+        e,
+        tvd,
+        shape,
+        locked=0,
+        orientation=0,
+        dip=0,
+        color='green',
+        alpha=0.5,
+        **geometry
+    ):
+        """
+        Parameters
+        ----------
+            geometry: 
+        """
+
+        SHAPES = [
+            'circle',
+            'ellipse',
+            'rectangle',
+            'polygon',
+        ]
+
+        GEOMETRY = dict(
+            rectangle = ['pos1', 'pos2'],
+            circle = ['radius'],
+            ellipse = {'radius_1': 0, 'radius_2': 0, 'res': 120},
+        )
+
+        assert shape in SHAPES, "shape not in SHAPES"
+        assert set(geometry.keys()) == set(GEOMETRY[shape]), 'wrong geometry'
+
+        self.name = name
+        self.n = n
+        self.e = e
+        self.tvd = tvd
+        self.shape = shape
+        self.locked = locked
+        self.orientation = orientation
+        self.dip = dip
+        self.color = color
+        self.alpha = alpha
+        self.geometry = geometry
+
+    def plot_data(self):
+        pos = [self.n, self.e, self.tvd]
+        if self.shape == "circle":
+            g = Circle(
+                pos=pos,
+                r=self.geometry['radius'],
+                c=self.color,
+                alpha=self.alpha,
+                # res=self.geometry['res']
+            )
+            g.flag = self.name
+            g.pos=[self.n, self.e, 0]
+            g.rotate(self.dip, point=pos, axis=(0,1,0))
+            g.rotate(self.orientation, point=pos, axis=(1,0,0))
+
+        return g
+
+
diff --git a/welleng/utils.py b/welleng/utils.py
index 461ee76..f5c9af8 100644
--- a/welleng/utils.py
+++ b/welleng/utils.py
@@ -121,7 +121,7 @@ def __init__(
             ) + self.start_xyz
         )
 
-def get_vec(inc, azi, r=1, deg=True):
+def get_vec(inc, azi, nev=False, r=1, deg=True):
     """
     Convert inc and azi into a vector.
 
@@ -144,8 +144,13 @@ def get_vec(inc, azi, r=1, deg=True):
     y = r * np.sin(inc_rad) * np.cos(azi_rad)
     x = r * np.sin(inc_rad) * np.sin(azi_rad)
     z = r * np.cos(inc_rad)
-    
-    return np.array([x,y,z]).T
+
+    if nev:
+        vec = np.array([y, x, z]).T
+    else:    
+        vec = np.array([x,y,z]).T
+
+    return vec / np.linalg.norm(vec)
 
 def get_nev(poss, start_xyz=[0,0,0], start_nev=[0,0,0]):
     """
@@ -168,22 +173,29 @@ def get_nev(poss, start_xyz=[0,0,0], start_nev=[0,0,0]):
     
     return (np.array([n, e, v]).T + np.array([start_nev]))
 
-def get_angles(vec):
+def get_angles(vec, nev=False):
     '''
     Determines the inclination and azimuth from a vector.
 
     Params:
         vec: (n,3) array of floats
+        nev: boolean (default: False)
+            Indicates if the vector is in (x,y,z) or (n,e,v) coordinates.
 
     Returns:
-        [inc, azi]: (2,n) array of floats
+        [inc, azi]: (n,2) array of floats
             A numpy array of incs and azis in radians
 
     '''
     # make sure it's a unit vector
-    vec = vec / np.linalg.norm(vec, axis=-1)
+    vec = vec / np.linalg.norm(vec, axis=-1).reshape(-1,1)
     vec = vec.reshape(-1,3)
 
+    # if it's nev then need to do the shuffle
+    if nev:
+        y, x, z = vec.T
+        vec = np.array([x, y, z]).T
+
     xy = vec[:,0] ** 2 + vec[:,1] ** 2
     inc = np.arctan2(np.sqrt(xy), vec[:,2]) # for elevation angle defined from Z-axis down
     azi = (np.arctan2(vec[:,0], vec[:,1]) + (2 * np.pi)) % (2 * np.pi)
diff --git a/welleng/version.py b/welleng/version.py
index 63eb0cb..0876fc4 100644
--- a/welleng/version.py
+++ b/welleng/version.py
@@ -1 +1 @@
-__version__ = '0.1.6'
\ No newline at end of file
+__version__ = '0.1.7'
\ No newline at end of file
diff --git a/welleng/visual.py b/welleng/visual.py
index beefaaa..05fa5dc 100644
--- a/welleng/visual.py
+++ b/welleng/visual.py
@@ -1,5 +1,5 @@
 import trimesh
-from vedo import show, Box, Axes, trimesh2vedo, Lines, colorMap
+from vedo import show, Box, Axes, trimesh2vedo, Lines, colorMap, Arrows, Text
 import numpy as np
 
 from .version import __version__ as VERSION
@@ -23,7 +23,15 @@ def __init__(
             height
         ).wireframe()
 
-def plot(data, names=None, colors=None, lines=None):
+def plot(
+    data,
+    names=None,
+    colors=None,
+    lines=None,
+    targets=None,
+    arrows=None,
+    text=None
+):
     """
     A vedo wrapper for quicly visualizing well trajectories for QAQC purposes.
 
@@ -54,6 +62,7 @@ def plot(data, names=None, colors=None, lines=None):
                 "Colors must be length of meshes list, 1 else None"
 
     meshes_vedo = []
+    # TODO: if inly a single mesh is provided then drop the need for a list
     for i, mesh in enumerate(meshes):
         if i == 0:
             vertices = np.array(mesh.vertices)
@@ -89,6 +98,8 @@ def plot(data, names=None, colors=None, lines=None):
         meshes_vedo,
         w.world,
         lines,
+        targets,
+        arrows,
         axes,
         bg='lightgrey',
         bg2='lavender',